复化梯形公式和复化Simpson公式.docx
《复化梯形公式和复化Simpson公式.docx》由会员分享，可在线阅读，更多相关《复化梯形公式和复化Simpson公式.docx（9页珍藏版）》请在冰豆网上搜索。
复化梯形公式和复化Simpson公式
数值计算方法上机题目3
一、计算定积分的近似值：
要求：
（1）若用复化梯形公式和复化Simpson公式计算，要求误差限
，分别利用他们的余项估计对每种算法做出步长的事前估计；
（2）分别利用复化梯形公式和复化Simpson公式计算定积分；
（3）将计算结果与精确解比较，并比较两种算法的计算量。
1.复化梯形公式
程序：
程序1（求f（x）的n阶导数：
symsx
f=x*exp（x）%定义函数f（x）
n=input（'输入所求导数阶数:
'）
f2=diff（f,x,n）%求f（x）的n阶导数
结果1
输入n=2
f2=
2*exp（x）+x*exp（x）
程序2：
clc
clear
symsx%定义自变量x
f=inline（'x*exp（x）','x'）%定义函数f（x）=x*exp（x），换函数时只需换该函数表达式即可
f2=inline（'（2*exp（x）+x*exp（x））','x'）%定义f（x）的二阶导数,输入程序1里求出的f2即可。
f3='-（2*exp（x）+x*exp（x））'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值
e=5*10^（-8）%精度要求值
a=1%积分下限
b=2%积分上限
x1=fminbnd（f3,1,2）%求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值
forn=2:
1000000%求等分数n
Rn=-（b-a）/12*（（b-a）/n）^2*f2（x1）%计算余项
ifabs（Rn）break%符合要求时结束
end
end
h=（b-a）/n%求h
Tn1=0
fork=1:
n-1%求连加和
xk=a+k*h
Tn1=Tn1+f（xk）
end
Tn=h/2*（（f（a）+2*Tn1+f（b）））
z=exp
（2）
R=Tn-z%求已知值与计算值的差
fprintf（'用复化梯形算法计算的结果Tn='）
disp（Tn）
fprintf（'等分数n='）
disp（n）%输出等分数
fprintf（'已知值与计算值的误差R='）
disp（R）
输出结果显示：
用复化梯形算法计算的结果Tn=7.3891
等分数n=7019
已知值与计算值的误差R=2.8300e-008
2.Simpson公式
程序：
程序1：
（求f（x）的n阶导数）：
symsx
f=x*exp（x）%定义函数f（x）
n=input（'输入所求导数阶数:
'）
f2=diff（f,x,n）%求f（x）的n阶导数
结果1
输入n=4
f2=
4*exp（x）+x*exp（x）
程序2：
clc
clear
symsx%定义自变量x
f=inline（'x*exp（x）','x'）%定义函数f（x）=x*exp（x），换函数时只需换该函数表达式即可
f2=inline（'（4*exp（x）+x*exp（x））','x'）%定义f（x）的四阶导数，输入程序1里求出的f2即可
f3='-（4*exp（x）+x*exp（x））'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，一边求最大值
e=5*10^（-8）%精度要求值
a=1%积分下限
b=2%积分上限
x1=fminbnd（f3,1,2）%求负的四阶导数的最小值点，也就是求四阶导数的最大值点对应的x值
forn=2:
1000000%求等分数n
Rn=-（b-a）/180*（（b-a）/（2*n））^4*f2（x1）%计算余项
ifabs（Rn）break%符合要求时结束
end
end
h=（b-a）/n%求h
Sn1=0
Sn2=0
fork=0:
n-1%求两组连加和
xk=a+k*h
xk1=xk+h/2
Sn1=Sn1+f（xk1）
Sn2=Sn2+f（xk）
end
Sn=h/6*（f（a）+4*Sn1+2*（Sn2-f（a））+f（b））%因Sn2多加了k=0时的值，故减去f（a）
z=exp
（2）
R=Sn-z%求已知值与计算值的差
fprintf（'用Simpson公式计算的结果Sn='）
disp（Sn）
fprintf（'等分数n='）
disp（n）
fprintf（'已知值与计算值的误差R='）
disp（R）
输出结果显示：
用Simpson公式计算的结果Sn=7.3891
等分数n=24
已知值与计算值的误差R=2.7284e-008
用复化梯形公式计算的结果为：
7.3891，与精确解的误差为：
2.8300e-008。
等分数n=7019
用复化Simpson公式计算的结果为：
7.3891，与精确解的误差为：
2.7284e-008。
等分数n=24
3、柯斯特公式求积分：
程序代码：
（1）function[y,Ck,Ak]=NewtonCotes（fun,a,b,n）
ifnargin==1
[mm,nn]=size（fun）;
ifmm>=8
error（'为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！
'）
elseifnn~=2
error（'fun构成应为：
第一列为x，第二列为y，并且个数为小于10的等距节点！
'）
end
xk=fun（1,:
）;
fk=fun（2,:
）;
a=min（xk）;
b=max（xk）;
n=mm-1;
elseifnargin==4
xk=linspace（a,b,n+1）;
ifisa（fun,'function_handle'）
fx=fun（xk）;
else
error（'fun积分函数的句柄，且必须能够接受矢量输入！
'）
end
else
error（'输入参数错误，请参考函数帮助！
'）
end
Ck=cotescoeff（n）;
Ak=（b-a）*Ck;
y=Ak*fx';
（2）functionCk=cotescoeff（n）
fori=1:
n+1
k=i-1;
Ck（i）=（-1）^（n-k）/factorial（k）/factorial（n-k）/n*quadl（@（t）intfun（t,n,k）,0,n）;
end
（3）functionf=intfun（t,n,k）
f=1;
fori=[0:
k-1,k+1:
n]
f=f.*（t-i）;
end
代码解释：
function[y,Ck,Ak]=NewtonCotes（fun,a,b,n）
%y=NewtonCotes（fun,a,b,n）
%牛顿-科特斯数值积分公式
%参数说明：
%fun，积分表达式，这里有两种选择
%
（1）积分函数句柄，必须能够接受矢量输入，比如fun=@（x）sin（x）.*cos（x）
%
（2）x,y坐标的离散点，第一列为x，第二列为y，必须等距，且节点的个数小于9，比如：
fun=[1:
8;sin（1:
8）]'
%如果fun的表采用第二种方式，那么只需要输入第一个参数即可，否则还要输入a,b,n三个参数
%a，积分下限
%b，积分上限
%n，牛顿-科特斯数公式的阶数，必须满足1=8时不能保证公式的稳定性
%
（1）n=1，即梯形公式
%
（2）n=2，即辛普森公式
%（3）n=4，即科特斯公式
%y，数值积分结果
%Ck，科特斯系数
%Ak，求积系数
%
%Example
%fun1=@（x）sin（x）;%必须可以接受矢量输入
%fun2=[0:
0.1:
0.5;sin（0:
0.1:
0.5）];%最多8个点，必须等距
%y1=NewtonCotes（fun1,0,0.5,6）
%y2==NewtonCotes（fun2）
ifnargin==1
[mm,nn]=size（fun）;
ifmm>=8
error（'为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！
'）
elseifnn~=2
error（'fun构成应为：
第一列为x，第二列为y，并且个数为小于10的等距节点！
'）
end
xk=fun（1,:
）;
fk=fun（2,:
）;
a=min（xk）;
b=max（xk）;
n=mm-1;
elseifnargin==4
%计算积分节点xk和节点函数值fx
xk=linspace（a,b,n+1）;
ifisa（fun,'function_handle'）
fx=fun（xk）;
else
error（'fun积分函数的句柄，且必须能够接受矢量输入！
'）
end
else
error（'输入参数错误，请参考函数帮助！
'）
end
%计算科特斯系数
Ck=cotescoeff（n）;
%计算求积系数
Ak=（b-a）*Ck;
%求和算积分
y=Ak*fx';
functionCk=cotescoeff（n）
%由于科特斯系数最多7阶，为了方便我们可以直接使用，省得每次都计算
%A1=[1,1]/2
%A2=[1,4,1]/6
%A3=[1,3,3,1]/8
%A4=[7,32,12,32,1]/90
%A5=[19,75,50,50,75,19]/288
%A6=[41,216,27,272,27,216,41]/840
%A7=[751,3577,1323,2989,2989,1323,3577,751]/17280
%当时为了体现公式，我们使用程序计算n阶科特斯系数
fori=1:
n+1
k=i-1;
Ck（i）=（-1）^（n-k）/factorial（k）/factorial（n-k）/n*quadl（@（t）intfun（t,n,k）,0,n）;
end
functionf=intfun（t,n,k）
%科特斯系数中的积分表达式
f=1;
fori=[0:
k-1,k+1:
n]
f=f.*（t-i）;
end
输出结果：
fun=@（x）exp（x）;
a=-1;
b=1;
n=4;
NewtonCotes（fun,a,b,n）
ans=2.3505
二、三点数值微分
升级会员 