I want to extract w and x from this code, but unfortunately i cant and it gives me only 1 answer.

1 view (last 30 days)
function [sol]=numintla(fun,n)
[xx,w]=gaussla(6); % Call the gaussla.m file
sx=size(xx,1);
sol = 0;
for i=1:sx,
x=xx(i);
fx=eval(fun);
sol=sol+w(i)*fx;
end
end
function [xx,ww]=gaussla(n)
format long;
b(1)=1;
for i=1:n, % Calculation of the binomial coefficients
b(i+1)=(factorial(n))/(factorial(i)*(factorial(n+1-i)));
end
for i=1:n+1, % The polynomial coefficients
poly(i)=((-1)^(n+1-i))*b(i)*(factorial(n)/(factorial(n+1-i)));
end
xx=roots(poly);% The polynomial roots
for i=1:n, % Coefficients of the first derivative of the polynomial
polycd(i)=poly(i)*(n+1-i);
end
for i=1:n, % Evaluation
x=xx(i);
solde=0;
for k=1:n,
solde=solde+polycd(k)*(x^(n-k));
end
ww(i,1)=(factorial(n)^2)/(xx(i)*(solde^2));
end
end

Answers (2)

Walter Roberson
Walter Roberson on 11 Apr 2022
function [sol, w, x] = numintla(fun,n)
I would recommend that you reconsider your requirement. Your x variable is whatever was last written into x by your for i loop. It would make more sense if you were to return xx instead of x.

Riccardo Scorretti
Riccardo Scorretti on 11 Apr 2022
Hi Alireza,
I'm afraid that the fact your code returns only one value is the last of your problems. More specifically, you provides two functions: one returns two values and the other returns two. So, it depends on what you want to do: if your purpose is to get the quadrature formula, you can invoke guass1a in this way (the number 3 is just an example: you can use any positive integer):
[x, w] = gaussla(3)
In my opinion (but I may be wrong) the true problems are more fundamental. You seems to be willing to compute (gaussla) the nodes and the weights of a Gauss-Legendre quadrature formula, which is usually defined in the range [-1 1]. The function numintla integrates numerically the function fun over a not precised domain (which I suppose is [-1 1]).
I'm afraid the function gaussla returns the wrong result. For instance, you can check the first coefficients here: https://en.wikipedia.org/wiki/Gaussian_quadrature. For n = 1 and n = 2 you get respectively:
x =
1
w =
1
and
x =
1.000000000000000 + 1.000000000000000i
1.000000000000000 - 1.000000000000000i
w =
-0.500000000000000 + 0.500000000000001i
-0.500000000000000 - 0.500000000000001i
One observes that coefficients are complex (this is wrong). The error is likely ot be in the computation of the coefficients of Legendre polynomials (see https://en.wikipedia.org/wiki/Legendre_polynomials).

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!