高斯求积公式 matlab
1. 分别用三点和四点Gauss-Chebyshev公式计算积分
并与准确积分值2arctan4比较误差。若用同样的三点和四点Gauss-Legendre公式计算，也给出误差比较结果。
2*atan(4)
ans =
2.6516
Gauss-Chebyshev：
function I = gausscheby(f,a,b,n)
syms t;
t= findsym(sym(f));
ta = (b-a)/2;
tb = (a+b)/2;
switch n
case 3
I=pi/n*ta*(subs(sym(f),t,ta*cos(pi/(2*n))+tb)*sqrt(1-cos(pi/(2*n))^2)+...
subs(sym(f),t,ta*cos(3*pi/(2*n))+tb)*sqrt(1-cos(3*pi/(2*n))^2)+...
subs(sym(f),t,ta*cos(5*pi/(2*n))+tb)*sqrt(1-cos(5*pi/(2*n))^2));
case 4
I=pi/n*ta*(subs(sym(f),t,ta*cos(pi/(2*n))+tb)*sqrt(1-cos(pi/(2*n))^2)+...
subs(sym(f),t,ta*cos(3*pi/(2*n))+tb)*sqrt(1-cos(3*pi/(2*n))^2)+...
subs(sym(f),t,ta*cos(5*pi/(2*n))+tb)*sqrt(1-cos(5*pi/(2*n))^2)+...
subs(sym(f),t,ta*cos(7*pi/(2*n))+tb)*sqrt(1-cos(7*pi/(2*n))^2));
end
I=simplify(I);
I=vpa(I,6);
syms x
f=1/(1+x^2);
a=-4;b=4;
n=3;
y=gausscheby(f,a,b,n)
y =
4.511
N=4:
y =
1.90041
Gauss-Legendre：
function I = IntGaussLegen(f,a,b,n)
syms t;
t= findsym(sym(f));
ta = (b-a)/2;
tb = (a+b)/2;
switch n
case 0,
I=2*ta*subs(sym(f),t,tb);
case 1,
I=ta*(subs(sym(f),t,ta*0.5773503+tb)+...
subs(sym(f),t,-ta*0.5773503+tb));
case 2,
I=ta*(0.55555556*subs(sym(f),t,ta*0.7745967+tb)+...
0.55555556*subs(sym(f),t,-ta*0.7745967+tb)+...
0.88888889*subs(sym(f),t,tb));
case 3,
I=ta*(0.3478548*subs(sym(f),t,ta*0.8611363+tb)+...
0.3478548*subs(sym(f),t,-ta*0.8611363+tb)+...
0.6521452*subs(sym(f),t,ta*0.3398810+tb) +...
0.6521452*subs(sym(f),t,-ta*0.3398810+tb));
case 4,
I=ta*(0.2369269*subs(sym(f),t,ta*0.9061793+tb)+...
0.2369269*subs(sym(f),t,-ta*0.9061793+tb)+...
0.4786287*subs(sym(f),t,ta*0.5384693+tb) +...
0.4786287*subs(sym(f),t,-ta*0.5384693+tb)+...
0.5688889*subs(sym(f),t,tb));
case 5,
I=ta*(0.1713245*subs(sym(f),t,ta*0.9324695+tb)+...
0.1713245*subs(sym(f),t,-ta*0.9324695+tb)+...
0.3607616*subs(sym(f),t,ta*0.6612094+tb)+...
0.3607616*subs(sym(f),t,-ta*0.6612094+tb)+...
0.4679139*subs(sym(f),t,ta*0.2386292+tb)+...
0.4679139*subs(sym(f),t,-ta*0.2386292+tb));
end
I=simplify(I);
I=vpa(I,6);
y =
2.04798
N=4:
y =
3.08862
2. 分别用三点和四点Gauss-Lagurre公式计算积分
function I = GaussLagurre(f,n)
syms t;
t= findsym(sym(f));
switch n
case 2
I=0.7110930*subs(sym(f),t,0.4157746)+...
0.2785177*subs(sym(f),t,2.2942804)+...
0.0103893*subs(sym(f),t,6.2899451);
case 3
I=0.6031541*subs(sym(f),t,0.3225477)+...
0.3574187*subs(sym(f),t,1.7457611)+...
0.0388879*subs(sym(f),t,4.5366203) +...
0.0005393*subs(sym(f),t,9.3950710);
end
I=simplify(I);
I=vpa(I,6);
syms x
f=exp(-10*x)*sin(x);
f=f./exp(-x);
a=0;b=inf;
n=2;
y= GaussLagurre(f,n)
y =
0.00680897
N=4:
y =
0.0104892
3. 设，分别取，，用以下三个公式计算，
列表比较三个公式的计算误差，从误差可以得出什么结论？
function [df1,df2,df3,w1,w2,w3]=MidPoint(func,a)
if (nargin == 3 && h == 0.0)
disp('h不能为0');
return;
end
for k=1:6
h=1/10^k;
y0=subs(sym(func), findsym(sym(func)),a);
y1 = subs(sym(func), findsym(sym(func)),a+h);
y2 = subs(sym(func), findsym(sym(func)),a-h);
df1(k) = (y1-y0)/h;
df2(k) = (y1-y2)/(2*h);
y3=subs(sym(func), findsym(sym(func)),a+2*h);
y4=subs(sym(func), findsym(sym(func)),a-2*h);
df3(k)=(y4-8*y2+8*y1-y3)/(12*h);
w1(k)=1/a-df1(k);
w2(k)=1/a-df2(k);
w3(k)=1/a-df3(k);
end
df1=simplify(df1); df1=vpa(df1,6);
df2=simplify(df2); df2=vpa(df2,6);
df3=simplify(df3); df3=vpa(df3,6);
w1=simplify(w1); w1=vpa(w1,6);
w2=simplify(w2); w2=vpa(w2,6);
w3=simplify(w3); w3=vpa(w3,6);
syms x
f=log(x);
a=0.7;
[y1,y2,y3,w1,w2,w3]=MidPoint(f,a)
y1 =
[ 1.33531, 1.41846, 1.42755, 1.42847, 1.42856, 1.42857]
y2 =
[ 1.43841, 1.42867, 1.42857, 1.42857, 1.42857, 1.42857]
y3 =
[ 1.42806, 1.42857, 1.42857, 1.42857, 1.42857, 1.42857]
w1 =
[ 0.0932575, 0.0101079, 0.00101944, 0.000102031, 0.0000102043, 0.00000101738]
w2 =
[ -0.00983893, -0.0000971936, -9.71814e-7, -9.72193e-9, 1.92131e-10, -3.02526e-9]
w3 =
[ 0.000513166, 4.76342e-8, 8.04334e-12, -1.10191e-10, 3.37991e-10, -4.5859e-9]
