Matlab Code for the Gauss Legendre Quadrature

I need help contsructing the code for the gauss legendre quadrature using n =2 , 4, and 6. i was able to get the value for n =2 but after that im kind of lost.

%% parameters
a = -1; % lower bound
b = 1; % upper bound 
n = [-1/sqrt(3) 1/sqrt(3)]; %location values for n=2
n1 = [-0.86113 -0.33998 0.33998 0.86113]; %location values for n=4
n2 = [-0.93246 -0.66120 -0.23861 0.23861 0.66120 0.993246]; % location values for n=6
w = 1; %weight for n=2
w1 = [0.34785 0.625214]; %weights for n=4
w2 = [0.17132 0.36076 0.6791]; % weights for n=6
f = @(x) (1)/(x.^2 * (sqrt(x.^2 +1)));%function 
%% Gauss_Legendre Quadrature 
% n=2
for i = 1:length(n)
    h = n(i);
    y(i) =  w*f(h);   
end
gaussleg1 = sum(y); % sum of wi*fi for n=2
%n=4
%n=6

Answer 1

[Tommy]

Have your vectors of weights be the same size as your vectors of locations, ordered so that corresponding weights and locations are in the same position:

%% parameters
a = -1; % lower bound
b = 1; % upper bound 
n = [-1/sqrt(3) 1/sqrt(3)]; %location values for n=2
n1 = [-0.86113 -0.33998 0.33998 0.86113]; %location values for n=4
n2 = [-0.93246 -0.66120 -0.23861 0.23861 0.66120 0.993246]; % location values for n=6
w = [1 1]; %weight for n=2
w1 = [0.34785 0.625214 0.625214 0.34785]; %weights for n=4
w2 = [0.17132 0.36076 0.6791 0.6791 0.36076 0.17132]; % weights for n=6
f = @(x) (1)./(x.^2 .* (sqrt(x.^2 +1)));%function 
%% Gauss_Legendre Quadrature 
% n=2
gaussleg = sum(w.*f(n));
%n=4
gaussleg1 = sum(w1.*f(n1));
%n=6
gaussleg2 = sum(w2.*f(n2));

Higher precision in your weights and locations will give more accurate results.

(edit) I also changed f slightly so that you could pass in the entire vector of locations at once.

[Jose Puentes]

thanks this worked good !!

[Lainie Suarez]

When I run this in matlab, it says,"

Array indices

must be positive

integers or

logical values.

Error in gausslegone (line 13)

gaussleg = sum(w.*f(n));"

[John D'Errico]

@Lainie Suarez

Almost certainly you have created a variable with the name sum. DON'T DO THAT!!!!!

Look at the list of variables in your workspace. Think about what happens now when MATLAB tries to use the function sum. It gets confused, and then you get that equally confusing error message.

[Lainie Suarez]

I see what happened now, thank you. Now I am now getting the right answer.

Answer 2

[Manikanta Balaji]

Y=input('Enter function directly f(x) =','s');
f=inline(Y);
a=input('Enter Initial interval point ''a'' =');
b=input('Enter final interval point ''b'' =');
n=input('Enter Gauss point ''n'' =');
if (n>=2 && n<=9)
    %% data
    %abscissa values
    X=[-0.57735 -0.7746 -0.86114 -0.90618 -0.932470 -0.949108 -0.96029 -0.96816
        0.57735  0.0000 -0.33998 -0.53847 -0.661209 -0.741531 -0.79666 -0.83603
        0.00000  0.7746  0.33998  0.00000 -0.238619 -0.405845 -0.52553 -0.61337
        0.00000  0.0000  0.86114  0.53847  0.238619  0.000000 -0.18343 -0.32425
        0.00000  0.0000  0.00000  0.90618  0.661209  0.405845  0.18343  0.00000
        0.00000  0.0000  0.00000  0.00000  0.932470  0.741531  0.52553  0.32425
        0.00000  0.0000  0.00000  0.00000  0.000000  0.949108  0.79666  0.61337
        0.00000  0.0000  0.00000  0.00000  0.000000  0.000000  0.96029  0.83603
        0.00000  0.0000  0.00000  0.00000  0.000000  0.000000  0.00000  0.96816];
    %weight values
    w=[1.00000 0.5555 0.34785 0.23693 0.171324 0.129484 0.10123 0.08127
        1.00000 0.8888 0.65214 0.47863 0.360761 0.279705 0.22238 0.18064
        0.00000 0.5555 0.65214 0.56889 0.467914 0.381830 0.31370 0.26061
        0.00000 0.0000 0.34785 0.47863 0.467914 0.417959 0.36268 0.31234
        0.00000 0.0000 0.00000 0.23693 0.360761 0.381830 0.36268 0.33024
        0.00000 0.0000 0.00000 0.00000 0.171324 0.279705 0.31370 0.31234
        0.00000 0.0000 0.00000 0.00000 0.000000 0.129484 0.22238 0.26061
        0.00000 0.0000 0.00000 0.00000 0.000000 0.000000 0.10123 0.18064
        0.00000 0.0000 0.00000 0.00000 0.000000 0.000000 0.00000 0.08127];
    %% Gauss Quadrature
    G=0;
    F=@(t) f(((b-a)*t+(b+a))/2);
    for i=1:n
        G=G+w(i,n-1)*F(X(i,n-1))*(b-a)/2;
    end
    fprintf('Integration by %d Gauss Quadrature = %d',n,G);
else
    fprintf('Formula not available for %d point',n)
end

Answer 3

[Mushimiyimana]

Y=input('Enter function directly f(x) =','s');

f=inline(Y);

a=input('Enter Initial interval point ''a'' =');

b=input('Enter final interval point ''b'' =');

n=input('Enter Gauss point ''n'' =');

if (n>=2 && n<=9)

%% data

%abscissa values

X=[-0.57735 -0.7746 -0.86114 -0.90618 -0.932470 -0.949108 -0.96029 -0.96816

0.57735 0.0000 -0.33998 -0.53847 -0.661209 -0.741531 -0.79666 -0.83603

0.00000 0.7746 0.33998 0.00000 -0.238619 -0.405845 -0.52553 -0.61337

0.00000 0.0000 0.86114 0.53847 0.238619 0.000000 -0.18343 -0.32425

0.00000 0.0000 0.00000 0.90618 0.661209 0.405845 0.18343 0.00000

0.00000 0.0000 0.00000 0.00000 0.932470 0.741531 0.52553 0.32425

0.00000 0.0000 0.00000 0.00000 0.000000 0.949108 0.79666 0.61337

0.00000 0.0000 0.00000 0.00000 0.000000 0.000000 0.96029 0.83603

0.00000 0.0000 0.00000 0.00000 0.000000 0.000000 0.00000 0.96816];

%weight values

w=[1.00000 0.5555 0.34785 0.23693 0.171324 0.129484 0.10123 0.08127

1.00000 0.8888 0.65214 0.47863 0.360761 0.279705 0.22238 0.18064

0.00000 0.5555 0.65214 0.56889 0.467914 0.381830 0.31370 0.26061

0.00000 0.0000 0.34785 0.47863 0.467914 0.417959 0.36268 0.31234

0.00000 0.0000 0.00000 0.23693 0.360761 0.381830 0.36268 0.33024

0.00000 0.0000 0.00000 0.00000 0.171324 0.279705 0.31370 0.31234

0.00000 0.0000 0.00000 0.00000 0.000000 0.129484 0.22238 0.26061

0.00000 0.0000 0.00000 0.00000 0.000000 0.000000 0.10123 0.18064

0.00000 0.0000 0.00000 0.00000 0.000000 0.000000 0.00000 0.08127];

%% Gauss Quadrature

G=0;

F=@(t) f(((b-a)*t+(b+a))/2);

for i=1:n

G=G+w(i,n-1)*F(X(i,n-1))*(b-a)/2;

end

fprintf('Integration by %d Gauss Quadrature = %d',n,G);

else

fprintf('Formula not available for %d point',n)

end