数值微分与数值积分

关于数值微分与数值积分第1页，讲稿共47页，2023年5月2日，星期三（1）向前差商由泰勒（Taylor）展式：向前差商的截断误差为：2第2页，讲稿共47页，2023年5月2日，星期三（2）向后差商向后差商的截断误差阶为：（3）中心差商3由泰勒（Taylor）展式：第3页，讲稿共47页，2023年5月2日，星期三中心差商的截断误差阶为：4第4页，讲稿共47页，2023年5月2日，星期三2.1.2插值型数值微分以为插值点的插值多项式为当时误差项：5第5页，讲稿共47页，2023年5月2日，星期三例2.2给定三点并且求解过三点的插值多项式为6第6页，讲稿共47页，2023年5月2日，星期三得三点公式：7第7页，讲稿共47页，2023年5月2日，星期三对定义在区间［a,b］上的定积分以上公式多称为牛顿-莱布尼兹公式，F(x)为f(x)的原函数.但有时原函数不能用初等函数表示，有时原函数又十分复杂，难于求出或计算.如被积函数为:等函数的积分都无法解决，当被积函数为一组数据时，更是无能为力.为解决定积分的近似计算,从定积分的定义：2.2.1求积公式8§2.2数值积分第8页，讲稿共47页，2023年5月2日，星期三结束9

这样就避开了求原函数的运算.上式就叫做求积公式，Ak(k=0,1,…,n)与函数f(x)无关，叫做求积系数,显然要确定一个求积公式，要确定求积结点xk和求积系数Ak，或者说不同的求积结点和求积系数将确定不同的求积公式.2.2.2求积公式的余项和代数精度一般情况下，上式两端并不相等.我们称:为求积公式

的余项，或截断误差.第9页，讲稿共47页，2023年5月2日，星期三结束10

为考查一个求积公式的误差，通常用代数精度来表示，如果一个求积公式对于不超过m次的多项式都能够精确成立(R[f]≡0)，而对m+1次以上的多项式不能精确成立，则称该求积公式的代数精度为m.例如求积公式：验证当f(x)=xm，m=0,1,2,3,4时，是否有R[xm]=0第10页，讲稿共47页，2023年5月2日，星期三结束11所以以上求积公式的代数精度为3.

任何一个求积公式的代数精度至少为零,即取f(x)=1时公式应精确成立，这是求积系数应满足的起码条件，可以用它检验一个求积公式的系数的正确性.求代数精度的另一种方法：利用余项公式求代数精度m.如：若取：则而有所以求积公式的代数精度为3.第11页，讲稿共47页，2023年5月2日，星期三结束122.2.3插值型数值积分

由插值可知，对任一函数f(x)(包括表格形式的函数)可用一n次多项式对其插值，即当Pn

(x)为拉格朗日插值多项式时，即第12页，讲稿共47页，2023年5月2日，星期三结束13其中:通常将公式(2.1)叫做插值型求积公式.（2.2）称为求积公式的余项或误差。第13页，讲稿共47页，2023年5月2日，星期三例：建立[0,2]上以x0=0,x1=0.5,x2=2为节点的数值积分公式。解

第14页，讲稿共47页，2023年5月2日，星期三

为便于上机计算，通常在求积公式中取等距节点，即将积分区间[a,b]n等分，即令h=(b-a)/n，且记x0=a,xn=b，则节点为xk=x0+kh(k=0,1,…n)，作变换:t=(x-x0)/h，代入求积系数公式：结束15这种等距节点的求积公式通常叫做牛顿-柯特斯公式,下面介绍几个常用的公式:2.2.4Newton-Cotes积分第15页，讲稿共47页，2023年5月2日，星期三取a=x0,b=x1,(即n=1)，代入(2.3)式得结束1梯形公式所以梯形公式为16从图2.1看到，这是用梯形面积近似代替曲边梯形的面积，梯形公式的误差估计为：第16页，讲稿共47页，2023年5月2日，星期三ya0xb图2.1梯形求积公式的余项：梯形求积公式的代数精度：取所以，代数精度为1。第17页，讲稿共47页，2023年5月2日，星期三结束18例

利用梯形公式计算解:2抛物形（辛卜生）公式取a=x0,(a+b)/2=x1,b=x2,(即n=2)，代入(2.3)式得第18页，讲稿共47页，2023年5月2日，星期三结束19所以抛物形（Simpson）公式为其中h=(b-a)/2,上式也可写成:Simpson公式的余项：第19页，讲稿共47页，2023年5月2日，星期三结束20或例2

利用抛物形公式计算解:Simpson公式的代数精度：取所以，代数精度为3。第20页，讲稿共47页，2023年5月2日，星期三结束(2.3)式给出213牛顿-柯特斯公式其中:可以看出，C(n)k不依赖函数f(x)和积分区间[a,b]，可以事先计算出来，通常叫做牛顿-柯特斯系数，下面给出n从1～6的牛顿--柯特斯系数表：第21页，讲稿共47页，2023年5月2日，星期三结束22n11/21/221/64/61/631/83/83/81/847/9032/9012/9032/907/90519/28875/28850/28850/28875/28819/288641/840216/84027/840272/84027/840216/84041/840Cotes系数表第22页，讲稿共47页，2023年5月2日，星期三结束当n为奇数，f(x)∈Cn+1［a,b］时，求积公式的余项为:23代数精度为n.4牛顿-柯特斯公式的余项与代数精度当n为偶数，f(x)∈Cn+2［a,b］时，求积公式的余项为:代数精度为n+1.第23页，讲稿共47页，2023年5月2日，星期三结束§2.3复化求积公式24

高阶的Newton-Cotes公式不能保证数值积分的系列的收敛性；同时，高阶Newton-Cotes积分的计算是不稳定的。因而，在实际应用中往往采用

将积分区间划分成若干个小区间，在各小区间上采用低次的求积公式(梯形公式或抛物形公式)，然后再利用积分的可加性，把各区间上的积分加起来，便得到新的求积公式.

这就是复化求积公式的基本思想。第24页，讲稿共47页，2023年5月2日，星期三结束2.3.1.复化梯形公式用n+1个分点将区间[a，b]n等分。每个区间长

在[xk，xk+1]上用梯形公式，则25Tn叫做复化梯形求积公式，下标n表示将积分区间等分的份数.第25页，讲稿共47页，2023年5月2日，星期三结束公式的特点：节点xk(k=1,2,…n-1)作为小区间的端点参与前、后两个小区间的计算，因而系数为2，端点a与b只参与一次计算，系数为1.26如果在Tn的基础上，将各小区间对分，这时节点数为2n+1，分段数为2n.记新的分点的函数值的和为σn，则T2n应为原内节点与新增节点函数值的和的两倍，加上两端点a,b的函数值之和再乘上新区间长度的一半，即第26页，讲稿共47页，2023年5月2日，星期三结束定理

设f(x)∈C2[a,b]，复化梯形公式的截断误差复化梯形积分的余项其中：h=(b-a)/n27从这一公式可以看出，将区间对分后，原复化梯形公式的值Tn作为一个整体保留.只需计算出新分点的函数值，便可得出对分后的积分值，不需重复计算原节点的函数值，从而减少了计算量.

利用复化求积公式的余项，我们可以估计出在满足精度的要求下，应将积分区间等分多少份，即n取多少.这种误差估计方法称为事前误差估计.第27页，讲稿共47页，2023年5月2日，星期三例

利用复化梯形公式计算使其误差限为10-4，应将区间［0,1］几等分?结束解:因为被积函数28取n=48可满足要求.第28页，讲稿共47页，2023年5月2日，星期三结束29将积分区间[a,b]2m等分，n=2m，节点为xk=a+kh(k=0,1,2,…,2m)，h=(b-a)/2m.在每两个小区间［x2k,x2k+2］(k=0,1,2,…,m-1)上用抛物形公式，则有:2.3.2复化抛物形（Simpson）公式S2m叫做复化抛物形求积公式，下标2m表示积分区间等分的份数，2m强调为偶数份.第29页，讲稿共47页，2023年5月2日，星期三结束30定理

设函数f(x)∈C4[a,b]，则Simpson的余项为公式的特点为节点x2k,(k=1,2,…,m-1)作为小区间[x2k,x2k+2]的端点，参与前后两次的辛普生公式的计算，因而系数为2，而奇数节点x2k+1,(k=0,1,…,m-1)因辛普生公式中间点的求积系数为4而保留4，前面的h/3为辛普生公式的公共求积系数.复化抛物形（Simpson）公式的余项和代数精度Simpson的代数精度为3第30页，讲稿共47页，2023年5月2日，星期三结束31解:利用例3的结果因此只需将区间［0,1］四等分，即取m=2(n=4).例4

利用复化抛物形公式计算使其误差限为10-4，应将区间［0,1］几等分?第31页，讲稿共47页，2023年5月2日，星期三2.3.3复化积分的自动控制误差算法在计算中常用来估计误差。（1）的计算公式对于取分点n=1,h1=(b-a),得：对分(取分点n=2),h2=(b-a)/2得：32第32页，讲稿共47页，2023年5月2日，星期三再对分(取分点n=4),h4=(b-a)/4得：33……，再对分(取分点n=2m),h2m=(b-a)/(2m),得：或：其中第33页，讲稿共47页，2023年5月2日，星期三34（2）计算公式的误差控制比较若要只需后者比前者便于控制第34页，讲稿共47页，2023年5月2日，星期三复化梯形积分的算法1.输入误差控制精度tol，初始分点数n，端点a,b,

h=(b-a)/n

2.建立函数文件f(x)

3.T2=(f(a)+f(b))*h/2;T1=T2+100;6.S=Sum(S1);T2=T1/2+h*S;ClearS1

end

5.Fori=1:n/2

计算S1(i)=f(a+(2*i-1)*h)end

4.While|T2-T1|>tolT1=T2h=h/2;n=2n;%区间分半357.输出T2

第35页，讲稿共47页，2023年5月2日，星期三结束362.3.4龙贝格(Romberg)积分公式我们已知的T2n与Tn的关系1.Romberg积分公式于是可以逐次对分形成一个序列{T1,T2,T4,T8,…},此序列收敛于积分真值I.当|T2n-Tn|<ε时,取T2n为I的近似值.以上算法称为复化梯形公式的逐次分半公式.但由于此序列收敛太慢,因此并不实用.现我们试图将它改造成为收敛快的序列.第36页，讲稿共47页，2023年5月2日，星期三第37页，讲稿共47页，2023年5月2日，星期三结束38如认为则有于是有:记这样我们从收敛较慢的{Tn}序列推出了收敛较快的{Sn}序列.截断误差由O(h2)提高到O(h4).如认为则有于是有:记这样我们从{Sn}序列又推出了收敛更快的{Cn}柯特斯序列.截断误差提高到O(h6),代数精度为5.第38页，讲稿共47页，2023年5月2日，星期三结束39如认为则有于是有:记这样我们从{Cn}序列又推出了收敛更快的{Rn}序列.{Rn}序列也称为龙贝格序列.截断误差提高到O(h8),代数精度为7.Romberg积分公式：对每一个k，j从2做到k，一直做到｜Tk,k-Tk,k-1｜<tol时停止计算。第39页，讲稿共47页，2023年5月2日，星期三结束40T1T2S1T4S2C1T8S4C2R1T16S8C4R2﹕﹕﹕﹕表7-2该表四个序列都是收敛的.第40页，讲稿共47页，2023年5月2日，星期三结束41例5

利用龙贝格方法计算解:计算结果列如下表:i2iT序列S序列C序列R序列013.00000123.100003.13333243.131183.141573.14212383.138993.141593.141593.141594