Matlab Code for the Gauss Legendre Quadrature

564 次查看(过去 30 天)
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

采纳的回答

Tommy
Tommy 2020-4-16
编辑:Tommy 2020-4-16
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.
  4 个评论
John D'Errico
John D'Errico 2022-3-28
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
Lainie Suarez 2022-3-29
I see what happened now, thank you. Now I am now getting the right answer.

请先登录,再进行评论。

更多回答(2 个)

Manikanta Balaji
Manikanta Balaji 2022-6-9
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

Mushimiyimana
Mushimiyimana 2024-1-15
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

类别

Help CenterFile Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by