数值计算方法上机题目3一、计算定积分的近似值：221x e xe dx =⎰ 要求：（1）若用复化梯形公式和复化Simpson 公式计算，要求误差限71021−⨯=ε，分别利用他们的余项估计对每种算法做出步长的事前估计；（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：clcclearsyms 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:1000000 %求等分数nRn=-(b-a)/12*((b-a)/n)^2*f2(x1) %计算余项if abs(Rn)<e %用余项进行判断break% 符合要求时结束endendh=(b-a)/n %求hTn1=0for k=1:n-1 %求连加和xk=a+k*hTn1=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)输出结果显示：用复化梯形算法计算的结果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))','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 %求等分数nRn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1) %计算余项if abs(Rn)<e %用余项进行判断break% 符合要求时结束endendh=(b-a)/n %求hSn1=0Sn2=0for k=0:n-1 %求两组连加和xk=a+k*hxk1=xk+h/2Sn1=Sn1+f(xk1)Sn2=Sn2+f(xk)endSn=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=243、柯斯特公式求积分：程序代码：（1）function [y,Ck,Ak]=NewtonCotes(fun,a,b,n)if nargin==1[mm,nn]=size(fun);if mm>=8error('为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！')elseif nn~=2error('fun构成应为：第一列为x，第二列为y，并且个数为小于10的等距节点！')endxk=fun(1,:);fk=fun(2,:);a=min(xk);b=max(xk);n=mm-1;elseif nargin==4xk=linspace(a,b,n+1);if isa(fun,'function_handle')fx=fun(xk);elseerror('fun积分函数的句柄，且必须能够接受矢量输入！') endelseerror('输入参数错误，请参考函数帮助！')endCk=cotescoeff(n);Ak=(b-a)*Ck;y=Ak*fx';（2）function Ck=cotescoeff(n)for i=1:n+1k=i-1;Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/n*quadl(@(t)int fun(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)% 牛顿-科特斯数值积分公式% 参数说明：% 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<n<7，因为n>=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)if nargin==1[mm,nn]=size(fun);if mm>=8error('为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！')elseif nn~=2error('fun构成应为：第一列为x，第二列为y，并且个数为小于10的等距节点！')endxk=fun(1,:);fk=fun(2,:);a=min(xk);b=max(xk);n=mm-1;elseif nargin==4% 计算积分节点xk和节点函数值fxxk=linspace(a,b,n+1);if isa(fun,'function_handle')fx=fun(xk);elseerror('fun积分函数的句柄，且必须能够接受矢量输入！')endelseerror('输入参数错误，请参考函数帮助！')end% 计算科特斯系数Ck=cotescoeff(n);% 计算求积系数Ak=(b-a)*Ck;% 求和算积分y=Ak*fx';function Ck=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阶科特斯系数for i=1:n+1k=i-1;Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/n*quadl(@(t)intfun( t,n,k),0,n);endfunction 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二、三点数值微分。