1、1,5.3 Romberg算法,综合前几节的内容,我们知道,梯形公式,Simpson公式,Cotes公式的代数精度分别为,1次,3次和5次,复合梯形、复合Simpson、复合Cotes公式的收敛阶分别为,2阶、4阶和6阶,无论从代数精度还是收敛速度,复合梯形公式都是较差的,有没有办法改善梯形公式呢?,2,一、复合梯形公式的递推化,各节点为,复合梯形(Trapz)公式为,-(1),-(2),3,-(3),4,则由(1)(2)(3)式,有,5,因此(1)(2)(3)式可化为如下递推公式,(4)-,上式称为递推的梯形公式,6,三种公式之间的关系,Romberg算法求解步骤,（3，-3）（-1,2）（
2、0,1）（1，-1）（3，-4）,7,二、外推加速公式,由复合梯形公式的余项公式,可得,由(3)式,8,复合Simpson公式,-(5),-(6),9,因此由复合Simpson公式的余项,可得,即,当然,令,自己证明,-(6),-(7),10,-(8),即,当然,同样由复合Cotes公式的余项,得,令,-(9),11,外推加速公式,以上整个过程称为Romberg算法,将上述结果综合后,12,其中外推加速公式可简化为,-(9),Romberg算法求解步骤,13,romberg.m,例:,前侧矩形公式 z1 = 0.99212545660563 z11 = 0.99212545660563后侧矩形
3、公式 z2 = 1.00783341987358 z22 = 1.00783341987358梯形公式 z3 = 0.99997943823961Simpson公式 z4 = 1.000000002016138阶simpson公式 z5 =1.00000000000000自选步长梯形公式 z6 = 0.99999921563419自选步长Simpson公式 z7 =1.00000051668471Romberg公式 z8 = 0.99999999999802Mote-Carlo算法 z9 = 0.99821071589516,-0.00787454339437 -0.0078745433943
4、7 0.00783341987358 0.00783341987358 -0.00002056176039 0.00000000201613 -0.00000000000000 -0.00000078436581 0.00000051668471 -0.00000000000198 -0.00178928410484,Jifenbijiao.m,积分法,积分值,绝对误差,14,如何构造Romberg算法,3 龙贝格(Romberg)积分方法,我们已经知道,当被积函数f(x)在区间a,b上连续时,要使得复合梯形公式或复合抛物线公式比较精确地代替定积分 可将分点(即基点)加密,也就是将区间a,b细
5、分,然后利用复合梯形公式或复合抛物线公式求积。,若用Tm表示把a,b作m等分并按复合梯形公式求积的结果,将每一小段再对分,令新的小段的长h=h/2,则T2m与Tm之间有如下关系:,(528),其中,另外,若用Sm表示把a,b分成m(偶数)个小段按复合抛物线公式计算的结果,那么只要把Sm中的m改为2m,h改为h就有,从Tm的定义可得到关系式,(529),我们再举一个计算上半单位圆面积的例子(它的准确面积为/2)。现用内接正多边形的逼近方法来计算。 如图5.6,图(a) 、 图(b)是用同样的内接正多边形计算上半单位圆的面积。图(a)是用梯形方法计算其面积,图(b)是用三角形方法计算其面积。,图
6、5.6,设正多边形边数为n=2k,则由图(b)利用三角形公式算得面积为,同理,如果组合一下,就会得到更精确的结果,即,同理,再以类似方法组合得,这样继续下去,其值越来越接近上半单位圆面积/2。这种方法可以用到计算定积分,为了推广公式(529)和上述计算上半单位圆面积的组合方法,我们引进龙贝格求积算法。 龙贝格求积算法本来是利用所谓外推法构造出的一种计算积分的方法。为了避免从外推引入而带来理论上的麻烦,我们将直接从构造一个T数表开始。 首先将a,b依次作20,21,22,等分,记,按复合梯形公式(520)算得的值相应地记为T(k)0(k=0,1,2,);把按式(529)算得的S2m依次记为T(k
7、)1(k=0,1,2,崐),而这每一个S2m又理解为由T2m与Tm的线性组合得到的改进值,即,我们可按照类似的方法继续进行改进,也即由S2m与Sm的线性组合得到改进值,依次记为T(k)2(k=0,1,2,),即,这样就可构造出一个数表,(5-30),其中除第0列(即最左一列)的T(k)0是按复合梯形公式计算外,其余各列都按下述规则(对m),(531),递推地计算出来。箭头表示计算流程。其计算步骤为: (1)将区间a,b等分为20,用梯形公式计算T(0)0,即,(2)将区间a,b等分为21,用梯形公式算出T(1)0,即,再由T(0)0,T(1)0根据公式(531)算出T(0)1,即,若 T(0)
8、1-T(0)0, (为预给的精度) 则停止计算;否则继续往下计算;,(3)依次分别算出T(2)0,T(1)1,T(0)2,这一行地往下推算,每一行算完,就得验证T(0)m(m=1,2,)是否满足预给的精度,即若,则停止计算;否则继续进行下一行。为了便于在计算机上实现,可运用下列公式编制程序:,例 4 计算积分,精确到10-4。 解,于是,由于,实际上,简单的数值方法与基本概念,1. 简单欧拉法（Euler） 2后退的欧拉法 3梯形法 4改进Euler法,2、初值问题的数值解法单步法,1. 简单的欧拉(Euler)方法,考虑模型：,在精度要求不高时,通过欧拉方法的讨论,弄清常微方程初值问题数值解
9、法的一些基本概念和构造方法的思路.,欧拉方法,最简单而直观实用方法,2. 欧拉方法的导出,把区间a,b,分为n个小区间,步长为,要计算出解函数 y(x) 在一系列节点,节点,处的近似值,N等分,对微分方程(1.1)两端从,进行积分,右端积分用左矩形数值求积公式,得,亦称为欧拉折线法 /* Eulers polygonal arc method*/,每步计算,只用到,或用向前差商近似导数,依上述公式逐次计算可得：,故也称Euler为单步法。,公式右端只含有已知项,所以又称为显格式的单步法。,例1 用欧拉公式求解初值问题,解 取步长h=0.1，欧拉公式的具体形式为,其中xn=nh=0.1n (n=
10、0,1,10), 已知y0 =1, 由此式可得,依次计算下去，部分计算结果见下表.,与准确解 相比，可看出欧拉公式的计算结果精度很差.,欧拉公式具有明显的几何意义, 就是用折线近似代替方程的解曲线，因而常称公式(2.1)为欧拉折线法.,还可以通过几何直观来考察欧拉方法的精度.假设yn=y(xn),即顶点Pn落在积分曲线y=y(x)上，那么，,按欧拉方法做出的折线PnPn+1便是y=y(x)过点Pn的切线.从图形上看,这样定出的顶点Pn+1显著地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的.,Ri 的主项/* leading term */, 欧拉法的局部截断误差：,欧拉法具有 1 阶精度。,
11、5. 欧拉公式的改进：, 隐式欧拉法 /* implicit Euler method */,由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公式。,一般先用显式计算一个初值，再迭代求解。, 隐式欧拉法的局部截断误差：,即隐式欧拉公式具有 1 阶精度。,设用欧拉公式,给出迭代初值 ，用它代入(2.5)式的右端，使之转化为显式，直接计算得,然后再用 代入(2.5)式，又有,如此反复进行，得,6.梯形公式 /* trapezoid formula */, 显、隐式两种算法的平均,注：的确
12、有局部截断误差 ， 即梯形公式具有2 阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。, 中点欧拉公式 /* midpoint formula */,假设 ，则可以导出即中点公式具有 2 阶精度。,需要2个初值 y0和 y1来启动递推过程，这样的算法称为双步法 /* double-step method */，而前面的三种算法都是单步法 /* single-step method */。,简单,精度低,稳定性最好,精度低, 计算量大,精度提高,计算量大,精度提高, 显式,多一个初值, 可能影响精度,改进的欧拉公式,我们看到，梯形方法虽然
13、提高了精度，但其算法复杂，在应用迭代公式(2.9)进行实际计算时，每迭代一次，都要重新计算函数 f(x, y )的值，而迭代又要反复进行若干次，计算量很大，而且往往难以预测. 为了控制计算量，通常只迭代一两次就转入下一步的计算，这就简化了算法.,具体地说，我们先用欧拉公式求得一个初步的近似值 ，称之为预测值，此预测值 的精度可能很差，再用梯形公式(2.7)将它校正一次，即按(2.8)式迭代一次，这个结果称之为校正值., 改进欧拉法 /* modified Eulers method */,注：此法亦称为预测-校正法 /* predictor-corrector method */。可以证明该算法具有 2 阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。,
