数值计算方法上机题目3
一、计算定积分的近似值：
2
21x e xe dx =⎰ 要求：
（1）若用复化梯形公式和复化Simpson 公式计算，要求误差限7102
1
-⨯=ε，分别利用他们的余项估计对每种算法做出步长的事前估计；
（2）分别利用复化梯形公式和复化Simpson 公式计算定积分；
（3）将计算结果与精确解比较，并比较两种算法的计算量。

1.复化梯形公式
程序：
程序1（求f （x ）的n 阶导数：
syms x
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
syms x%定义自变量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值
for n=2:1000000 %求等分数n
Rn=-(b-a)/12*((b-a)/n)^2*f2(x1) %计算余项
if abs(Rn)<e %用余项进行判断
break% 符合要求时结束
end
end
h=(b-a)/n %求h
Tn1=0
for k=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阶导数）：
syms x
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
syms x%定义自变量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值
for n=2:1000000 %求等分数n
Rn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1) %计算余项
if abs(Rn)<e %用余项进行判断
break% 符合要求时结束
end
end
h=(b-a)/n %求h
Sn1=0
Sn2=0
for k=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)
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,:);
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。