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

1、复化梯形公式和复化Simpson公式数值计算方法上机题目3一、计算定积分的近似值：要求：（1）若用复化梯形公式和复化Simpson公式计算，要求误差限，分别利用他们的余项估计对每种算法做出步长的事前估计；（2）分别利用复化梯形公式和复化Simpson公式计算定积分；（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：clcclear

2、syms 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*10(-8) %精度要求值 a=1 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值for n=2:100000

3、0 %求等分数n Rn=-(b-a)/12*(b-a)/n)2*f2(x1) %计算余项 if abs(Rn)e %用余项进行判断 break % 符合要求时结束 endendh=(b-a)/n %求hTn1=0 for k=1:n-1 %求连加和 xk=a+k*h Tn1=Tn1+f(xk)endTn=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)输出结果显示：用复化

4、梯形算法计算的结果 Tn= 7.3891等分数 n=7019已知值与计算值的误差 R= 2.8300e-0082. Simpson公式程序：程序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) + x*exp(x),

5、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值for n=2:1000000 %求等分数n Rn=-(b-a)/180*(b-a)/(2*n)4*f2(x1) %计算余项 if abs(Rn)=8 error(为了保证NewtonCotes积分的稳定性，最多只能有9个等距节

6、点！) elseif nn=2 error(fun构成应为：第一列为x，第二列为y，并且个数为小于10的等距节点！) end xk=fun(1,:); fk=fun(2,:); a=min(xk); b=max(xk); n=mm-1; elseif nargin=4 xk=linspace(a,b,n+1); if isa(fun,function_handle) fx=fun(xk); else error(fun积分函数的句柄，且必须能够接受矢量输入！) end else error(输入参数错误，请参考函数帮助！) end Ck=cotescoeff(n); Ak=(b-a)*Ck;

7、y=Ak*fx; （2）function Ck=cotescoeff(n) for i=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）function f=intfun(t,n,k) f=1; for i=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) % 牛顿-科特斯数值积分公式 % 参数说明： % fu

8、n，积分表达式，这里有两种选择 %(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，牛顿-科特斯数公式的阶数，必须满足1n=8时不能保证公式的稳定性 % (1)n=1，即梯形公式 % (2)n=2，即辛普森公式% (3)n=4，即科特斯公式 % y，数值积分结果 % Ck，科特斯系数

9、 % 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) if nargin=1 mm,nn=size(fun); if mm=8 error(为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！) elseif nn=2 error(fun构成应为：第一列为x，第二列为y，并且个数为小于10的等距节点！) end xk=fun(1,:);

10、 fk=fun(2,:); a=min(xk); b=max(xk); n=mm-1; elseif nargin=4 % 计算积分节点xk和节点函数值fx xk=linspace(a,b,n+1); if isa(fun,function_handle) fx=fun(xk); else error(fun积分函数的句柄，且必须能够接受矢量输入！) end else error(输入参数错误，请参考函数帮助！) end % 计算科特斯系数 Ck=cotescoeff(n); % 计算求积系数 Ak=(b-a)*Ck; % 求和算积分 y=Ak*fx; function Ck=cotescoe

11、ff(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阶科特斯系数 for i=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 function f=intfun(t,n,k) % 科特斯系数中的积分表达式 f=1; for i=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二、三点数值微分