复化梯形公式和复化Simpson公式.docx
- 文档编号:5697446
- 上传时间:2022-12-31
- 格式:DOCX
- 页数:9
- 大小:19.80KB
复化梯形公式和复化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 % (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 二、三点数值微分
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 梯形 公式 Simpson
![提示](https://static.bdocx.com/images/bang_tan.gif)