复化梯形公式和复化Simpson公式
吉林大学工程数学计算方法(第三章习题答案)
第三章习题答案1.分别用梯形公式、Simpson公式、Cotes公式计算积分1,I=⎰并估计误差。
解：1）用梯形公式有：()()110.51[10.5]10.42678242f f⎛-≈+=+≈⎝⎭⎰()()()333333220.512.6042107.36571012124Tb aE f fηηη-----⎛⎫''=-=--=⨯≤⨯⎪⎝⎭事实上，()()()()()()110.430964410.50.510.4267767210.50.510.00418772Tf x II f fE f f f===-≈+=⎡⎤⎣⎦-∴=-+=⎡⎤⎣⎦⎰⎰2）Simpson公式()110.53111410.43093 642122f f f⎛-⎡⎤⎛⎫⎛⎫≈++=+=⎪ ⎪⎢⎥⎝⎭⎝⎭⎣⎦⎝⎭⎰[]()()44744211111522 1.1837710180218028Sb a b aE f fηη--⎛⎫--⎪⎛⎫--⎛⎫=-=--≤⨯⎪ ⎪⎪⎝⎭⎝⎭⎪⎝⎭3122()''()48T f fb aE事实上，()()()10.510.50.510.5410.000030462SE f f f f-⎡+⎤⎛⎫=-++=⎪⎢⎥⎝⎭⎣⎦⎰3）由Cotes公式有：()() ()111537270.5321232719084814.9497525.2982210.3923029.9332670.43096180f f f f f-⎡⎤⎛⎫⎛⎫⎛⎫≈++++⎪ ⎪ ⎪⎢⎥⎝⎭⎝⎭⎝⎭⎣⎦=++++=⎰15732127)18088()6116211294522 2.697410945464C E f η--⎛⎫⨯ ⎪⎛⎫=-⨯-≤⨯ ⎪ ⎪⎝⎭⎪⎝⎭7(6)945*42()()82Cf b aEf事实上，()0.0000003C E f =2．证明Simpson 公式()2.8具有三次代数精度。
复化梯形公式,复化辛普森公式,复化柯特斯公式
复化梯形公式,复化辛普森公式,复化柯特斯公式
复化梯形公式、复化辛普森公式和复化柯特斯公式都是用来计算定积分的近似值的方法。
1. 复化梯形公式：将积分区间分成若干个小区间，在每个小区间上用梯形面积近似代替该小区间的曲边梯形面积，然后将这些梯形面积相加，得到积分的近似值。
2. 复化辛普森公式：将积分区间分成若干个等分小区间，在每个小区间上用矩形面积近似代替该小区间的曲边梯形面积，然后将这些矩形面积相加，得到积分的近似值。
3. 复化柯特斯公式：将积分区间分成若干个等分小区间，在每个小区间上用切线段长度近似代替该小区间的曲边梯形面积，然后将这些切线段长度相加，得到积分的近似值。
这三种方法都是通过将积分区间分成若干个小区间，然后在每个小区间上用近似方法计算该小区间的曲边梯形面积，最后将这些近似值相加得到积分的近似值。
它们的精度和误差都与分区间的大小有关。
复化求积公式
h[ 2
f ( x0 )
n1
2
k 1
f ( xk )
f ( xn )]
复化梯形公式
Tn
h 2
[
f
(
x0
)
n1
2
k 1
f ( xk )
f ( xn )]
复化梯形公式
计算方法
2.复化辛浦生公式
计算方法
在每个小区间[xk1, xk ]上应用辛浦生公式得：
xk
xk 1
则
f
b
( x)dx h[ f 6
计算方法
在 每 个 小 区 间[ xk1, xk ]上 应 用 梯 形 公 式 得 ：
xk xk 1
f ( x)dx
h 2
[
f
(
xk1
)
f ( xk )]
则
b
n
f (x)dx =
xk f (x)dx
a
k 1 xk1
n k 1
h[ f 2
(xk1)
f
(xk )]
计算方法
x0 x1 x2 x3
2
三、区间逐次分半求积法
计算方法
复化求积公式可有效提高计算精度，但对给定 的误差限，如何确定节点的个数，即[a,b]应多少等 份？由截断误差可以估计步长的取值情况，但需要 给出各阶导数的最大值，这往往是比较困难的，且 估计值往往偏大.
接下来，我们将考虑步长的更为实用的选取方 法.
计算方法
若用Tn及T2n分别表示将[a, b]n等分及2n等分的复化 梯形公式，则
f(x) 1 0.997 0.9896 0.976 0.95 0.936 0.908 0.877 0.841 3978 158 7267 8851 1556 8516 1925 4709
复化梯形公式和复化Simpson公式
数值计算方法上机题目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。
3.1用复化梯形公式、复化Simpson公式、Romberg方法和复
(2) f (x) ex2 cos20x,0 x 2。
1
3.4 设计自适应的Simpson方法求积分0 x xdx( 0.4) 的
近似值，即对不同的子区间分别按精度标准确定各自适当的步长， 计算各子区间上的积分近似值，然后将各个近似值相加，要求近似 值的绝对误差限为0.5107 。
(2) 0 1 x2 dx 4
3.2 用外推方法计算下列积分值，并对计算结果进行比较。 如果所得结果不满意，对算法进行适当修改。
2x
x2ห้องสมุดไป่ตู้
(1)
1
( 1
x2
)dx 0.51324025, 2
(2) sin2 xdx
0
2
第三章 数值积分与数值微分
3.3 用样条函数方法和外推法求下列函数的一阶和二阶导数， 并结合函数的图形说明精度与步长h的关系。
第三章 数值积分与数值微分
数值试验题3
3.1 用复化梯形公式、复化Simpson公式、Romberg方法 和复化Gauss-Legendre公式计算下列积分的近似值，使绝对误 差限为 0.5107 ，并将计算结果与精确解作比较以及比较各种 算法的计算量。
21
(1)1
dx ln 2, x
11
复化梯形公式和复化辛普生公式
void main()
{
tixing ti;
ti.integration();
}
2）复化simpson
#include<iostream.h>
#include<math.h>
class simpson
{
private:
double result;
int n,k;//n是等分数，
double h;//h步长
cin>>b>>a;
cout<<"输入你要使用tixing法则的数目(即等分数)";
cin>>n;
h=(b-a)/n;
sum=0.0;
for(k=1;k<n;k++)
{
sum+=func(k*h);
}
integral=(h/2.0)*(2.0*sum+func(b)+1);
cout<<endl<<"积分值="<<integral<<endl;
public:
double sine(double);
void integration();//实现积分的函数
};
double simpson::sine(double a)
{
result=sin(a)/a;//被积函数
cout<<"sin("<<a<<")/"<<a<<"="<<result<<endl;
源程序：
牛顿柯特斯公式
I
ab
f
( x)dx
n1 xi1
i0 xi
f
( x)dx
n1h [
i02
f
(xi )
f
( xi 1 )]
h[ 2
f
(a)
n1
2 f
i1
(xi )
f
(b)].
记为
Tn
n1h [
i02
f
(xi )
复化梯形公式Tn的余项R
若f ( x)在[a, b]上连续，则
ba h
在每个小区间：[xk , xk1 ]上，共三个点： xk , xk1/ 2 , xk1 n
所以这里在[0, 1]上实际上共有5个分点。
精品课件！
精品课件！
例2的计算结果表明，为达到相同的精度，用复化
Simpson公式所需的计算量比复化梯形公式少，这也说明 了复化Simpson公式的精度较高，实际计算时多采用复化 Simpson公式。
6
2
I
b
f (x)dx
b x3dx b4 a4
a
a
4
而 S b a ( f (a) 4 f ( a b ) f (b)) b a (a3 4( a b )3 b3 )
6
2
6
2
1
(b 4
a4
a(a2
b2)
4
a
b
3
(b
(1)
0.9460832
C2
1 2 90
7
f
(0) 32[
数学实验题目2 Romberg积分法
数学实验题目2 Romberg 积分法摘要考虑积分()()b aI f f x dx =⎰欲求其近似值，可以采用如下公式：（复化）梯形公式 110[()()]2n i i i hT f x f x -+==+∑ 2()12b a E h f η-''=- [,]a b η∈ （复化）辛卜生公式 11102[()4()()]6n i i i i hS f x f x f x -++==++∑4(4)()1802b a h E f η-⎛⎫=- ⎪⎝⎭ [,]a b η∈ （复化）柯特斯公式 111042[7()32()12()90n i i i i hC f x f x f x -++==+++∑31432()7()]i i f xf x +++6(6)2()()9454b a h E f η-⎛⎫=- ⎪⎝⎭[,]a b η∈ 这里，梯形公式显得算法简单，具有如下递推关系121021()22n n n i i h T T f x -+==+∑因此，很容易实现从低阶的计算结果推算出高阶的近似值，而只需要花费较少的附加函数计算。
但是，由于梯形公式收敛阶较低，收敛速度缓慢。
所以，如何提高收敛速度，自然是人们极为关心的课题。
为此，记0,k T 为将区间[,]a b 进行2k等份的复化梯形积分结果，1,k T 为将区间[,]a b 进行2k等份的复化辛卜生积分结果，2,k T 为将区间[,]a b 进行2k等份的复化柯特斯积分结果。
根据李查逊（Richardson ）外推加速方法，可得到1,11,,0,1,2,40,1,2,41m m k m km k m k T T T m -+-=-⎛⎫=⎪=-⎝⎭可以证明，如果()f x 充分光滑，则有,lim ()m k k T I f →∞= （m 固定）,0lim ()m m T I f →∞=这是一个收敛速度更快的一个数值求积公式，我们称为龙贝格积分法。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值计算方法上机题目3一、计算定积分的近似值：e2= xe x dx要求：（1）若用复化梯形公式和复化Simpson 公式计算，要求误差限= 110-7，分别利用他们的余项估计对每种算法做出步长的事前估计；（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：clearsyms x %定义自变量xclcf=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*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)输出结果显示：用复化梯形算法计算的结果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 end end h=(b-a)/n% 用余项进行判断% 符合要求时结束%求 hSn1=0Sn2=0 for k=0:n-1xk=a+k*h%求两组连加和xk1=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==4 xk=linspace(a,b,n+1);if isa(fun,'function_handle')fx=fun(xk);elseerror('fun 积分函数的句柄，且必须能够接受矢量输入！') end elseerror（'输入参数错误，请参考函数帮助！'） 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)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)% 牛顿- 科特斯数值积分公式% 参数说明：% 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 个点，必须等距(--O T皿⅛、A 、X..UnmOEeu2uu4φs φ(--莊、番芝6皿遥 □(*嘛、⅛⅛s φ^SEeOUolMeNH ⅛p ^-^o h ω8U人IUlU⅛XUn4ΦN S丄 UirIUIU一T H H U -P e U 4一(CN u n 4)SECOUoIMeNH艮 ％(9Ln .o o 1: Un4) SEeOUoIMeNHIZ%A>lx)un4HX4 (^φ-puPlIIUO一I o Un4-⊂n4)es一七U+uqeφoedsu=H>lx X4Wi 協莊相 t ⅛>lx j π?相 0赵M ⅛%寸H H W P e U4φs φTIUIUHUAMX)XeIUHqVx)U-IUHe (Z)un4u>l4【3error('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二三点数值微分。