复化梯形公式和复化Simpson公式

1、数值计算方法上机题目3一、计算定积分的近似值：22 xe 二 xe dxi要求：(1) 若用复化梯形公式和复化Simpson公式计算，要求误差限10，分别利用他们的余项估计对每种算法做出步长的事前估2计；(2) 分别利用复化梯形公式和复化 Simps on公式计算定积分；(3) 将计算结果与精确解比较，并比较两种算法的计算量。1.复化梯形公式程序：程序1 (求f (x)的n阶导数：syms xf=x*exp(x)淀义函数 f (x)n=input(输入所求导数阶数:)f2=diff(f,x,n)%f(x)的 n阶导数结果1输入n=2f2 = 2*exp(x) + x*exp(x)程序2 :cl

2、cclearsyms x%定义自变量xf=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*10A(-8)%精度要求值a=1%积分下限b=2%积分上限x1=fmi nbn d(f3,1,2)%求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值f

3、or n=2:10000009求等分数nRn=-(b-a)/12*(b-a)/n)A2*f2(x1) if abs(Rn)ve breakendendh=(b-a)/nTn仁0for k=1:n-1%求连加和xk=a+k*h%计算余项%用余项进行判断%符合要求时结束%求 hTn 1=T n1+f(xk)endTn=h/2*(f(a)+2*T n1+f(b)z=exp(2)R=T n-z%求已知值与计算值的差fprintf(用复化梯形算法计算的结果Tn=)disp(T n)fprintf(等分数 n=)disp( n)%俞出等分数fprintf(已知值与计算值的误差 R=)disp(R)输出结果

4、显示：用复化梯形算法计算的结果 Tn二7.3891等分数n=7019已知值与计算值的误差R= 2.8300e-0082. Simps on 公式程序：程序1 :(求f (x )的n阶导数)：syms xf=x*exp(x)%定义函数 f (x)n=input(输入所求导数阶数:)f2=diff(f,x,n)%求f(x)的 n 阶导数结果1输入n=4f2 = 4*exp(x) + x*exp(x)程序2 :clcclearsyms x%定义自变量xf=inline(x*exp(x) , x ) %定义函数f(x)=x*exp(x),换函数时只需换该函数表达式即可f2=inline(4*exp(x

5、) + x*exp(x) , x )%定义 f(x)的四阶导数，输入程序 1 里求出的f2即可f3= -(4*exp(x) + x*exp(x)%因fminbnd ()函数求的是表达式的最小值，且要求表达式带引号，故取负号，一边求最大值e=5*10A(-8)a=1%积分下限b=2%积分上限x1=fmi nbn d(f3,1,2)最大值点对应的x值for n=2:1000000%精度要求值%求负的四阶导数的最小值点，也就是求四阶导数的%求等分数nRn=-(b-a)/180*(b-a)/(2* n)F4*f2(x1)%计算余项if abs(Rn)ve%用余项进行判断break%符合要求时结束end

6、%求 hend h=(b-a)/nSn仁0Sn2=0for k=0:n-1%求两组连加和xk=a+k*hxk仁xk+h/2Sn 1=S n1+f(xk1)Sn 2=S n2+f(xk)endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a)+f(b)%因Sn2多加了 k=0时的值，故减去 f(a)z=exp(2)R=S n-z%求已知值与计算值的差fprintf(用Simpson公式计算的结果 Sn=)disp(S n)fprintf(等分数 n=)disp (n)fprintf(已知值与计算值的误差 R=)disp(R)输出结果显示：用Simpson公式计算的结果 Sn二7.389

7、1等分数n=24已知值与计算值的误差R= 2.7284e-008用复化梯形公式计算的结果为：7.3891，与精确解的误差为：2.8300e-008。等分数 n=7019用复化Simpson公式计算的结果为：7.3891，与精确解的误差为:2.7284e-008。等分数 n=243、柯斯特公式求积分:程序代码：(1) function y,Ck,Ak二NewtonCotes(fun,a,b,n)if n arg in=1mm, nn =size(fu n);if mm=8error(为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！)elseif nn=2error(fun构成

8、应为：第一列为x，第二列为y,并且个数为小于10的等距节点！)endxk=fu n(1,:);fk=fu n( 2,:);a=mi n(xk);b=max(xk);n=m m-1;elseif nargin=4xk=li nspace(a,b, n+1);if isa(fu n,fun cti on_han dle)fx=fu n(xk);elseerror（fun积分函数的句柄，且必须能够接受矢量输入！）endelseerror(输入参数错误，请参考函数帮助！)endCk=cotescoeff( n);Ak=(b-a)*Ck;y=Ak*fx：(2) function Ck=cotescoef

9、f(n)for i=1: n+1k=i-1;Ck(i)=(-1)A( n-k)/factorial(k)/factorial( n-k)/n*quadl(t)i ntfun (t, n,k),0,n);end(3) function f二intfun(t,n,k)f=1;for i=0:k-1,k+1: nf=f.*(t-i);end代码解释:fun ctio n y,Ck,Ak二Newto nCotes(fu n,a,b ,n) % y=Newt on Cotes(f un, a,b ,n)%牛顿-科特斯数值积分公式%参数说明：% fun,积分表达式，这里有两种选择%(1)积分函数句柄，必须

10、能够接受矢量输入，比如fun二 (x)si n( x).*cos(x)% (2)x,y坐标的离散点，第一列为x，第二列为y，必须等距，且节点的个数小于9， 比如：fun二1:8;sin(1:8)%如果fun的表采用第二种方式，那么只需要输入第一个参数即可， 否则还要输入a,b,n三个参数%a，积分下限% b，积分上限% n，牛顿-科特斯数公式的阶数，必须满足1n=8时不能保证公式的稳定性%(1) n=1，即梯形公式%(2) n=2，即辛普森公式%(3) n=4，即科特斯公式%y，数值积分结果% Ck，科特斯系数% Ak，求积系数% Example% fun仁(x)sin(x);%必须可以接受矢

11、量输入% fun2二0:0.1:0.5;sin(0:0.1:0.5);% 最多 8 个点，必须等距% y仁N ewto nCotes(fu n1,0,0.5,6)% y2=Newto nCotes(fu n2)if n arg in=1mm, nn =size(f un);if mm=8error(为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！)elseif nn=2error(fun构成应为：第一列为x,第二列为y，并且个数为小于10的等距节点！)endxk=fu n(1,:);fk=fu n( 2,:);a=mi n(xk);b=max(xk);n=m m-1;els

12、eif nargin=4%计算积分节点xk和节点函数值fx xk=li nspace(a,b ,n+1);if isa(fu n,fu ncti on_han die)fx=fu n(xk);elseerror(fun积分函数的句柄，且必须能够接受矢量输入！)endelseerror(输入参数错误，请参考函数帮助！)end%计算科特斯系数Ck=cotescoeff( n);%计算求积系数Ak=(b-a)*Ck;%求和算积分y=Ak*fx;fun cti on Ck=cotescoeff( n)%由于科特斯系数最多7阶，为了方便我们可以直接使用，省 得每次都计算% A仁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阶科特斯系数for i=1: n+1k=i-1;Ck(i)=(-1F( n-k)/factorial(k)/factor