复化梯形公式和复化Simpson公式
大学精品课【工程计算基础】工程计算5数值积分和数值微分
Rtn
(b a)h2 12
f ()
(a,b)
2020/9/29
18
5.2 梯形公式和Simpson求积公式
如果采用Simpson公式，则须在区间[a，b]内有 2n +1个节点
a =x0< x1<…<x2n=b xk= a + kh k = 0，1，…，2n 其中，h = (ba)/2n
[a，b]分成n个子区间，每个子区间[x2n，x2n+2]
2020/9/29
4
5.2 梯形公式和Simpson求积公式
• 5.2.1 梯形公式和Simpson公式
采用插值原理构造数值求积公式。
对于两点插值{a，f(a)；b，f(b)}，其拉格朗日插值
公式和余项分别为
xb
xa
L1(x)
ab
f (a) ba
f (b)
R1 ( x)
1 2
f
( )(x a)(x b)
2020/9/29
10
5.2 梯形公式和Simpson求积公式
• 5.2.2牛顿-考特斯 (Newton-Cotes)公式
将区间[a，b]划分为n等分，步长h = (ba)/n，
选取等距节点xk=a+kh 构造出的插值性求积公式
n
In ( f ) (b a)
C(n) k
f
(
xk
)
k 0
称为牛顿-考特斯 (Newton-Cotes)公式。
12 其中，h =b a。
(a,b)
2020/9/29
6
5.2 梯形公式和Simpson求积公式
如果在[a，b]内取三个插值节点，{a，f(a)；(a +b)/2， f((a +b)/2)； b，f(b)}，则插值函数为
复化梯形公式,复化辛普森公式,复化柯特斯公式
复化梯形公式,复化辛普森公式,复化柯特斯公式
复化梯形公式、复化辛普森公式和复化柯特斯公式都是用来计算定积分的近似值的方法。
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公式
一、计算定积分的近似值：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=等分数 n=7019已知值与计算值的误差 R=2. 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=等分数 n=24已知值与计算值的误差 R=用复化梯形公式计算的结果为：，与精确解的误差为：。
复化Simpson公式
2012-2013（1）专业课程实践论文复化Simpson公式肖云龙，0818180214，R数学08-2班Simpson 公式是一个数值积分公式，在计算一些多项式函数（三次或三次一下）的定积分时会得出精确值。
但容易验证它对于)(x f =4x 通常是不准确的，因此，Simpson 公式实际上具有三次代数精度。
提高阶的途径并不总能取得满意的效果。
为了改善求积的精度，通常采用复化求积法。
Simpson 公式)]()2(4)([6b f ba f a f ab S +++-= 将定积分⎰=ba dx x f I )(的积分区间],[b a 分隔为n 等分，各节点为n j jh a x j ,,1,0, =+= nab h -=复化梯形公式的形式])()(2)([211∑-=++=n j j n b f x f a f hT记子区间[1,+j j x x ] 的中点为21+j x则复化Simpson 公式：n S ])(4)(2)()((6102111∑∑-=+-=+++=n j j n j j x f x f b f a f h#include<iostream.h>#include<math.h>double function(double x){ double s;s=x/(4+x*x);return s;}double ReiterationOfSimpson(double a,double b,double n,double f(double x)) { double h,fa,fb,xk,xj;h=(b-a)/n;fa=f(a);fb=f(b);double s1=0.0;double s2=0.0;for(int k=1;k<n;k++){ xk=a+k*h;s1=s1+f(xk);}for(int j=0;j<n;j++){ xj=a+(j+0.5)*h;s2=s2+f(xj);}double sn;sn=h/6*(fa+fb+2*s1+4*s2);return sn;}main(){double a,b,Result,n;cout<<"请输入积分下限："<<endl;cin>>a;cout<<"请输入积分上限："<<endl;cin>>b;cout<<"请输入分割区间数n："<<endl;cin>>n;cout<<"复化Simpson公式计算结果：";Result=ReiterationOfSimpson(a,b,n,function);cout<<Result<<endl;return 0;}例1．用Simpson 公式求积分⎰+1024x x解：运行程序(1) 请输入积分下限：输入a 的值为0，回车。
复化梯形公式和复化辛普生公式
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[
复化梯形求积公式
复化梯形求积公式
复化梯形求积公式是一种常用的数学计算方法，通常用于计算定积分的近似值。
该公式基于将被积函数的图像划分为若干个梯形，然后计算这些梯形的面积之和来近似计算定积分的值。
具体而言，复化梯形求积公式可以表示为：
∫a^b f(x)dx ≈ h/2 [f(a) + 2f(x1) + 2f(x2) + … + 2f(xn-1) + f(b)]
其中，a和b分别是定积分的上限和下限，f(x)是被积函数，h=(b-a)/n 是每个梯形的宽度，n是将被积函数图像分成的梯形个数，x1、x2、…、xn-1是每个梯形的右端点。
复化梯形求积公式的精度与分割后梯形的数量有关，通常情况下，分割的数量越多，计算结果越精确。
但也应该注意到，过多的分割会导致计算量增大，从而影响计算效率。
除了复化梯形求积公式外，还存在其他的数值积分方法，如复化矩形求积公式、辛普森公式等。
每种数值积分方法都有其适用的场景和特点，需要根据实际情况选择合适的方法进行计算。
总之，复化梯形求积公式是一种简单易用的数值积分方法，可以在不求解原函数的情况下近似计算定积分的值，具有较广泛的应用价值。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值计算方法上机题目3一、计算定积分的近似值：要求：（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)ifnargin==1[mm,nn]=size(fun);ifmm>=8error('为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！')elseifnn~=2error('fun构成应为：第一列为x，第二列为y，并且个数为小于10的等距节点！')endxk=fun(1,:);fk=fun(2,:);a=min(xk);b=max(xk);n=mm-1;elseifnargin==4xk=linspace(a,b,n+1);ifisa(fun,'function_handle')fx=fun(xk);elseerror('fun积分函数的句柄，且必须能够接受矢量输入！')endelseerror('输入参数错误，请参考函数帮助！')endCk=cotescoeff(n);Ak=(b-a)*Ck;y=Ak*fx';（2）functionCk=cotescoeff(n)fori=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）functionf=intfun(t,n,k)f=1;fori=[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)ifnargin==1[mm,nn]=size(fun);ifmm>=8error('为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！')elseifnn~=2error('fun构成应为：第一列为x，第二列为y，并且个数为小于10的等距节点！')endxk=fun(1,:);fk=fun(2,:);a=min(xk);b=max(xk);n=mm-1;elseifnargin==4%计算积分节点xk和节点函数值fxxk=linspace(a,b,n+1);ifisa(fun,'function_handle')fx=fun(xk);elseerror('fun积分函数的句柄，且必须能够接受矢量输入！')endelseerror('输入参数错误，请参考函数帮助！')end%计算科特斯系数Ck=cotescoeff(n);%计算求积系数Ak=(b-a)*Ck;%求和算积分y=Ak*fx';functionCk=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阶科特斯系数fori=1:n+1k=i-1;Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/n*quadl(@(t)intfun(t,n,k ),0,n);endfunctionf=intfun(t,n,k)%科特斯系数中的积分表达式f=1;fori=[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二、三点数值微分。