常微分方程的数值解及实验

1、1 引言 在常微分方程中，我们已经掌握了一些典型方程的解法。但许多形式的方程只能用数值方法求近似解，也就是求在某些点上满足一定精度的近似解。现以求一阶常微分方程初值问题00( , )()dyf x ydxy xy(41) 在区间a，b上的解为例，介绍数值方法的基本思想。 设f(x，y)在带形区域 R：axb，-y+ 上为x，y的连续函数，且对任意的y满足李普希茨(Libusize)条件 f(x，y1)-f(x，y2)Ly1-y2 (42) 其中(x，y1)、(x，y2)R，L为正常数。在求初值问题(41) 的数值解时，我们通常采用离散化方法(数值微分、数值积分、泰勒展式等)，求在自变量x的离散

2、点 a=x0 x1x2xn=b 图 4.1 上的准确解y(x)的近似值 y0，y1，y2，yn 常取离散点x0，x1，x2，xn为等距，即 x i+1-xi=h，i=0，1，2，n-1 h称为步长。图41表示为初值问题(41)在n+1个离散点上的准确解y(x)的近似值。2 欧拉法和改进的欧拉法 2.1 欧拉法(折线法) 若将函数y(x)在点xi处的导数y(xi)用差商来表示，即111()()()()()iiiiiiiy xy xy xxxy xy xh 再用yi近似地代替y(xi)，则初值问题(41) 就化为100( ,)()0,1,2,iiiiyyhf x yyy xi(43) 式(43)就

3、是所求的欧拉公式。 欧拉公式有很明显的几何意义。我们知道初值问题(41)中的微分方程的解是xoy平面上的一簇积分曲线，这簇积分曲线上任意点(x，y)?的斜率为f(x，y)，而初值问题(41)的解是过点(x0，y0)的一条特定的积分曲线。 例1 用欧拉法求初值问题2(0)0yxyy的数值解(取h=0.1)。 解 因为2100.1 ()0,0,1,2,iiiiyyxyyi2( , )(0)0,0.1f x yxyyh故由欧拉计算公式(43)得 (44) 表 41 图 4.2 图 4.3 2.2 改进的欧拉法 欧拉法虽然形式简单，计算方便，但比较粗糙，精度也低。特别当y=y(x)?的曲线曲率较大时，

4、欧拉法的效果更差。为了达到较高精度的计算公?式，对欧拉法进行改进，将在一点(xi，yi)的切线斜率f(xi?，yi)用两点的平均斜率来代替，即 111( ,) ( ,)(,)2iiiiiif x yf x yf xy 代入(43)式得 111 ( ,)(,)20,1,2,1iiiiiihyyf x yf xyin(45) 这样得到的点列仍为一折线，只是用平均斜率来代替原来一点处的斜率。式(45)称为改进的欧拉公式。 不难发现，欧拉公式(43)是关于yi+1的显式?，只要已知yi，经一次计算可立即得到yi+1的值;而改进的欧?拉公式(45)中的yi+1以隐式给出，且yi+1含在函数f(xi+1

5、， yi+1)中，因此?，通常用迭代法求解。具体做法是：先用欧拉公式(43)?求出一个y(0)i+1作为初始近似，然后再用改进的欧拉公式(45)进行迭代，即(0)1(1)( )111( ,) ( ,)(,)20,1,2,iiiikkiiiiiiyyhf x yhyyf x yf xyk 直到满足 (1)( )11(1)11kkiikiiyyyy(为预给精度)取 再转到下一步计算。 这里必须特别说明，因为初值问题(41)满足李普 希茨条件(1)( )( )(1)111111( )(1)11(,)22kkkkiiiiiikkiihyyf xxfyhLyy 当h足够小时，可使得 12hLq 于是有(

6、1)( )( )(1)11112(1)(2)11(1)(0)11kkkkiiiikkiikiiyyq yyqyyqyy当k时，有qk0，故公式(46)收敛。 2.3 预估校正法 改进的欧拉公式在实际计算时要进行多次迭代，因而计算量较大。所谓预估校正法，就是先用(43)式算出yi + 1的预估值y( p )i + 1，然后再用(45)式进行一次迭代便得到校正值y(c)i+1，即()1( )()111( ,) ( ,)(,2piiiicpiiiiiiyyhf x yhyyf x yf xy(47) 预估： 校正： 并取 ( )11ciiyy 虽然式(47)仅迭代一次，但因进行了预先估计，故精度却有

7、较大的提高。 在实际计算时，还常常将式(47)写成下列形式： 121112( ,)(,)()20,1,2,iiiiiikf x ykf xh yhkhyykki(48) 2.4误差估计 初值问题(41)的等价积分方程为11()()( , ( )iixiixy xy xf x y x dx(49) 若对式(49)右端的积分采用各种不同的近似计算方法，就可以得到初值问题(41)的各种不同的数值解法。 例 如积分采用左矩形公式111( , ( )( ,)()()()( ,)iixiiiixiiiif x y x dxf x yxxy xy xhf x y 图 4.4 用yi、yi+1分别代替y(xi

8、)、y(xi+1)便得到欧拉公式(43)。 若积分采用梯形公式11111( , ( ) ( ,)(,)()2iixiiiiiixf x y x dxf x yf xyxx 在进行误差分析时，我们假设yi=y(xi)，考虑用yi+1代替y(x i+1)而产生慕囟衔蟛睿康氖俏伺卸吓拉公式和改进的欧拉公式的精确度。 设初值问题(41)的准确解为y=y(x)，则利用泰勒公式 123()()()()()2!()3!iiiiiiy xy xhhy xhy xyxhyx 1. 欧拉公式的截断误差 由式(43)知 1( ,)()( ,)iiiiiiiyyhf x yy xhf x y(411) 比较式(410

9、)和(411)得2112()()2!()iiihy xyyxO h(412) 2. 改进的欧拉公式的截断误差 由式(45)知11111 ( ,)(,)2() ( , ()(,)2iiiiiiiiiiihyyf x yf xyhy xf x y xf xy(413) 对(49)式右端的积分采用梯形公式并根据梯形公式的误差可得到1113()() ( , ()(, ()2( )12iiiiiihy xy xf x y xf xy xhf(414)其中(xi，xi+1)，比较式(413)和(414)得1111113() (, ()(,)2( )12iiiiiihy xyf xy xf xyhf(415

10、) 因此 31111311311()()( )212(1)()( )12()()iiiiiiiihLhy xyy xyfhq y xyfy xyO h 所以，改进的欧拉公式的截断误差为O(h3)，也即改进的欧拉法为二阶的。 可以验证，预估校正公式(47)与改进的欧拉公式的截断误差相同，均为O(h3)。这里略去证明。 例 2求解初值问题 2(0)1,01,0.1xyyyyxh 解 现分别用欧拉公式和改进的欧拉公式进行计算。 这里欧拉公式的具体形式为121111222()()2iiiiiiiixkyyxhkyhkyhkhyykk 其解析解为 21yx 表 42 3 龙格库塔法 3.1 泰勒级数展开

11、法 我们还是假设yi=y(xi)，利用泰勒级数展开求y(xi+1)。式(410)就是y(xi+1)的泰勒展开式，若取右端前有限项作为y(xi+1)的近似值，就可得到计算y(xi+1)的各种不同截断误差的数值公式。 例 如，取前两项可得到 y(xi+1)y(xi)+hy(xi) =y(xi)+hf(xi，y(xi) =yi+hf(xi，yi) 若取前三项，可得到截断误差为O(h3)的公式 212()()()()2()( , ()( , ()( , ()( , ()2iiiiiiixiiiiyiihy xy xhy xyxy xhf x y xhfx y xf x y xfx y x 这里 y(x

12、i)=f(xi，y(xi) y(xi)=fx(xi，y(xi)+fy(xi，y(xi)y(xi) =fx(xi，y(xi)+f(xi，y(xi)fy(xi，y(xi) 类似地，若取前k项作为y(xi+1)的近似值，便得到截断误差为O(hk)的数值计算公式。这些公式的计算必须依赖于求y(xi)的k阶导数，除非f(x，y)足够简单，否则直接用泰勒展开法求解较为复杂。但是泰勒级数展开法的基本思想是许多数值方法的基础。 3.2 龙格库塔法 前面已经知道，初值问题(41)等价于11()()( , ( )()(, ()01iixiixiiiy xy xf x y x dxy xhf xh y xh 龙格库

13、塔法的基本思想是：用f(x，y)在几个不同点的数值加权平均来代替f(xi+h，y(xi+h)的值，而使截断误差的阶尽可能高。 1. 二阶龙格库塔公式 将预估校正公式(48)改写成更一般的形式11 122121()( ,)(,)iiiiiiyyhkkkf x ykf xh yhk(416) 适当选取%、1、2%的值，使截断误差y(x i+1)-y i+1的阶数尽可能高。这里仍假定yi=y(xi)，显然 1()iky x 2. 四阶龙格库塔公式 二阶龙格库塔公式是由使用在两个不同点上的函数值的线性组合而得到的。同样，我们用四个不同点上的函数值的线性组合就可得到四阶龙格库塔公式。设 yi+1=yi+

14、h(1k1+2k2+3k3+4k4) (420) 这里k1、k2、k3、k4为四个不同点上的函数值，分别设其为 k1=f(xi，yi) k2=f(xi+1h，yi+11k1h) k3=f(xi+2h，yi+21k1h+22k2h) k4=f(xi+3h，yi+31k1h+32k2h+33k3h) 其中1、2、3、4、1、2、3、1 1、2 1、2 2、3 1、32、33均为待定系数。 类似于前面的讨论，把k2、k3、k4分别在xi点展成h的幂级数，代入线性组合式(420)中，将得到的公式与y(xi+1)在xi点上的泰勒展开式比较，使其两式右端直到h4的系数相等，经过较复杂的运算便可得到关于i，

15、i，ij的一组特解 1=2=11=22=1/2 21=31=32=0 3=33=1 1=4=1/6 2=3=1/3 (422) 从而得到常用的标准四阶龙格库塔公式：121334311234( ,)(,)22(,)22(,)(22)6iiiiiiiiiikf x yhhkf xykhhkf xykkf xh yhkhyykkkk(423)图 4.5 表 434 阿达姆斯方法 我们已经知道，初值问题(41)等价于积分方程(49)，即11()()( , ( )nnxnnxy xy xf x y x dx 对积分式分别采用矩形公式和梯形公式可得到欧拉公式和改进的欧拉公式，截断误差分别为O(h2)和O(

16、h3)。为此，我们自然可以想到，若用更高次的插值多项式来代替f(x，y)，则所得公式的精度会更高。这就是线性多步法的起源思想。 本章前面介绍的方法称为单步法，因为在计算yi+1时，只用到前面yi的值。而对于线性多步法是要利用前面已经算出的若干个值yi-k，yi-1，yi来求yi+1。 现用k次多项式Pk(x)来代替f(x，y(x) 111()()( )( )iiiixxiikkxxy xy xP x dxRx dx(424) 舍去余项 并设yi=y(xi)，而yi+1为y(xi+1)的近似值，于是可得到线性多步法的计算公式1( )iixkkxRRx dx11( )iixiikxyyP k dx

17、 4.1 阿达姆斯(Adams)显式 取q+1个基点xi，xi-1，xi-q，并作牛顿后差插值多项式见式(438)，则0( )( 1)qmmqimtPxfm其中ixxth将式(427)代入式(426)得1100( 1)iiqxmmiiixmqmimimtyydxfmyhf (428) 这里 10( 1),0,1,2,mmtdtmqm(429) 式(428)称为阿达姆斯显式。 对于余项 1112(2)10( )( 1)()( )iixqqxqqtqqRRx dxhydt 112(2)10( 1)( )()(,)qqqtqqi qiRhydtxx 亦即 2(2)15(5)3( )251( )720

18、qqqqRhyRh y 显然当q=3时 (430) m是多项式积分，易算出结果如下： m01234m11/25/123/8251/270例如 130432101(1)(2)613()6 48t ttdtttt 为了易于在电子计算机上实现，常将式(428)中的 用各点的已知函数值表示。特别，当q=2时，有1121123(23165)12(5559379)24iiiiiiiiiiihyyfffhyyffff当q=3时，有 (431) (432) 4.2 阿达姆斯隐式 类似于4.1，取q+1个基点xi+1，xi，xi-q+1，并作牛顿后差插值多项式，则10( )( 1)qmmqimtP xfm(43

19、3) 其中 1ixxth 将式(433)代入式(426)得 1110011001( 1)( 1)( 1),0,1,2,iiqxmmiiixmqmmiimmmtyydtfmtyhdtfmtdt mqm (434)其中 (435) 式(434)称为阿达姆斯隐式。 类似于阿达姆斯显式余项的求法，可得到阿达姆斯隐式的余项为2(2)111( ),qqqqi qiRhyxx (436) 例当q=3时 5(5)319( )720Rh y m的计算结果如下： m01234m1-1/2-1/12-1/24-19/720 若将式(434)中各阶差分mfi+1用各点的已知函数值表示，则可得到便于在电子计算机上实现的

20、数值公式。 例如，当q=2时1111112(58)12(9195)24iiiiiiiiiiihyyfffhyyffff(437) (438) 当q=3时 4.3 阿达姆斯预估校正公式 我们常把阿达姆斯显式及隐式联立使用，即构造所谓阿达姆斯预估校正公式。现以q=2为例构造预估校正公式()111( )111(23165)12(58)12piiiiiciiiiihyyfffhyyfff(439) 并取 ( )11ciiyy 与同阶的龙格库塔方法相比较，阿达姆斯方法计算量小，公式简单，程序易于实现。但它的主要缺点是不能自动开始，开始的前几个值要依赖于其它方法获得。这里介绍两种计算开始值的方法。 (1)

21、 用单步法中的数值方法求出开始值。 (2) 使用y(x)的泰勒展开式 200000()( )()()()()2!yxy xy xy xxxxx 例4 用阿达姆斯方法求初值问题2(0)0,01,0.1yxyyxh(440) 的数值解。 解 首先用泰勒展式求其三个点的值，因为 220( )( )( )12 ( )( )( )2( )( )( )0y xxyxyxy x y xyxyxy x yxx 表 44 设常微分方程组的初值问题为0( )( , )()y xf x yy xs(441) 这里 21212( ),( ),( )( , )( , ),( , ),( , )( ,)TikTkTkyy

22、 xyxyxf x yf x yfx yfx yss ss 初值问题(441)与(41)式形式上完全相似。因此，对于(41)适用的数值计算公式，只要将其中的y0，y，f，都改写成相应的向量形式 s，y，f ，就能写出求解(441)的数值公式。 例如，初值问题(441)的标准四阶龙格库塔公式为112341213243(22)6( ,)(,)22(,)22(,)iiiiiiiiiihkykkkkkf x yhhkf xykhhkf xykkf xh yhk5二阶线性常微分方程边值问题的数值解 设二阶线性常微分方程的边值问题为 y+p(x)y+q(x)y=f(x) y(a)=，y(b)=，axb (

23、442) 其中p(x)，q(x)，f(x)为区间崐a，b上足够光滑的已知函数，且q(x)0，、为已知常数。 在上述条件下，边值问题(442)式存在连续可微的解，且是唯一的。若采用差分方法来解边值问题，其基本步骤是： (1)将区间a，b“离散化”，即给a，b一个分划，此分划常考虑等距; (2)对每一个基点，将各阶导数用差商来近似表示，将微分方程转化为差分方程，进而转化为线性代数方程组; (3)解线性代数方程组，求得各基点上的近似解。 现具体给出求解边值问题(442)的方法步骤。 首先将区间a，b进行等距分划，即令 xi=a+ih，i=0，1，2，n 其中bahn一般称 0nxaxb为边界点，称x

24、1，x2，xn-1为内 其次，在各基点xi上，将y，y用差商来近似表示。这里要求有相同阶数的截断误差，以保证精度协调。我们知道，由(463)式可得到211()()()()2iiiy xy xy xO hh(443) 又对式(458)再求一次导数，注意到( )1( )nndpxdpxdxhds有 230022401( )()(1)()6811()12npxf xsf xhttf x 令t=1，可得到 2410022100122211()()()1211()()( ()2 ()()nnpxf xf xhpxf xf xf xf xhh 取上式右边第一项作为pn(x1)的近似有 由差分与导数的关系可

25、得余项为2(4)( )12hf 所以，有 2112()2 ()()()()iiiiy xy xy xyxO hh(444) 略去(443)和(444)式的截断误差O(h2)，并用yi代替y(xi)可有 11112()22()iiiiiiiyyy xhyyyyxh(445)(446) 将(445)、(446)代入(442)得到近似差分方程为11112022,1,2,1iiiiiiiiinyyyyypq yfhhyyin(447) 其中 pi=p(xi)，qi=q(xi)，fi=f(xi)将式(447)整理后，得 其中 2212212,1,2,1iiiiiiiihapbh qhcpdh f in (449) 这个方程组的系数矩阵是三对角的，可以利用追赶法求解。 根据(445)、(446)式知，对于y(xi)-yi的误差方程只要把fi取成O(h2)即可。现剩下两个问题：其一是线性方程组(448)即差分方程(447) 解的存在唯一性;其二是(448)的解的收敛性，也即当h0时，解是否收敛于微分方程(442)