第十五章常微分方程的解法 建立微分方程只是解决问题的第一步,通常需要求出方程的解来说明实际现象,并 加以检验。如果能得到解析形式的解固然是便于分析和应用的,但是我们知道,只有线 性常系数微分方程,并且自由项是某些特殊类型的函数时,才可以肯定得到这样的解, 而绝大多数变系数方程、非线性方程都是所谓“解不出来”的,即使看起来非常简单的 方程如 a-y+x2,于是对于用微分方程解决实际问题来说,数值解法就是一个十 分重要的手段。 §1常微分方程的离散化 下面主要讨论一阶常微分方程的初值问题,其一般形式是 dy dx f(x,y)a≤x≤b y(a)=yo 在下面的讨论中,我们总假定函数f(x,y)连续,且关于y满足李普希兹 Lipschitz)条 件,即存在常数L,使得 If(x, y)-f(x,y)ksLLy-y 这样,由常微分方程理论知,初值问题(1)的解必定存在唯一。 所谓数值解法,就是求问题(1)的解y(x)在若干点 处的近似值yn(n=1,2,…,N)的方法,yn(n=1,2,…,N)称为问题(1)的数值解 bn=xn+-xn称为由xn到xn+1的步长。今后如无特别说明,我们总取步长为常量h 建立数值解法,首先要将微分方程离散化,一般采用以下几种方法 (i)用差商近似导数 若用向前差商 代替y(xn)代入(1)中的微分方程,则得 y(xmD)-y(u) ≈f(xn,y(xn)(n=0,1,…) 化简得 y(x+)=y(x, )+hf(x,,y(x,)) 如果用y(xn)的近似值yn代入上式右端,所得结果作为y(xn+)的近似值,记为yn+1, 则有 y1=yn+的f(xn,yn)(n=01… (2) 这样,问题(1)的近似解可通过求解下述问题 ∫yn+=yn+竹/(xn,yn)(n=01…) (3) lyo=y(a) 得到,按式(3)由初值y可逐次算出y1,y2…。式(3)是个离散化的问题,称为差 分方程初值问题
-179- 第十五章 常微分方程的解法 建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并 加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线 性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以肯定得到这样的解， 而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的 方程如 2 2 y x dx dy = + ，于是对于用微分方程解决实际问题来说，数值解法就是一个十 分重要的手段。 §1 常微分方程的离散化 下面主要讨论一阶常微分方程的初值问题，其一般形式是 ⎪ ⎩ ⎪ ⎨ ⎧ = = ≤ ≤ 0 ( ) ( , ) y a y f x y a x b dx dy (1) 在下面的讨论中，我们总假定函数 f (x, y) 连续，且关于 y 满足李普希兹(Lipschitz)条 件，即存在常数 L ，使得 | f (x, y) − f (x, y) |≤ L | y − y | 这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。 所谓数值解法，就是求问题（1）的解 y(x) 在若干点 a = x0 < x1 < x2 <L< xN = b 处的近似值 y (n 1,2, ,N) n = L 的方法， y (n 1,2, ,N) n = L 称为问题（1）的数值解， n n n h = x − x +1 称为由 n x 到 n+1 x 的步长。今后如无特别说明，我们总取步长为常量h 。 建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： （i）用差商近似导数 若用向前差商 h y x y x n n ( ) ( ) +1 − 代替 '( ) n y x 代入（1）中的微分方程，则得 ( , ( )) ( 0,1, ) ( ) ( ) 1 ≈ = L + − f x y x n h y x y x n n n n 化简得 ( ) ( ) ( , ( )) n 1 n n n y x ≈ y x + hf x y x + 如果用 ( ) n y x 的近似值 n y 代入上式右端，所得结果作为 ( ) n+1 y x 的近似值，记为 n+1 y ， 则有 ( , ) ( 0,1, ) yn+1 = yn + hf xn yn n = L （2） 这样，问题（1）的近似解可通过求解下述问题 ⎩ ⎨ ⎧ = + = + = ( ) ( , ) ( 0,1, ) 0 1 y y a y y hf x y n n n n n L （3） 得到，按式（3）由初值 0 y 可逐次算出 y1, y2 ,L。式（3）是个离散化的问题，称为差 分方程初值问题
需要说明的是,用不同的差商近似导数,将得到不同的计算公式。 (i)用数值积分方法 将问题(1)的解表成积分形式,用数值积分方法离散化。例如,对微分方程两端 积分,得 y(xn)-y(x)2=f(xyx)dx(m=01…) 右边的积分用矩形公式或梯形公式计算。 (i) Taylor多项式近似 将函数y(x)在xn处展开,取一次 Taylor多项式近似,则得 y(,1)=y(x,)+hy (x,)=y(x)+hf(n, y(xu)) 再将y(xn)的近似值yn代入上式右端,所得结果作为y(xn+1)的近似值yn+1,得到离 散化的计算公式 ym1=yn+hf(xn, yu) 以上三种方法都是将微分方程离散化的常用方法,每一类方法又可导出不同形式的 计算公式。其中的 aylor展开法,不仅可以得到求数值解的公式,而且容易估计截断 误差 §2欧拉( Euler)方法 2.1Euer方法 Euler方法就是用差分方程初值问题(3)的解来近似微分方程初值问题(1)的解, 即由公式(3)依次算出y(x)的近似值yn(n=12…)。这组公式求问题(1)的数值 解称为向前 Euler公式。 如果在微分方程离散化时,用向后差商代替导数,即y(xm)≈(xn)-y(x) h 则得计算公式 )(n=01…) (5) yo=y(a) 用这组公式求问题(1)的数值解称为向后Euer公式。 向后 Euler法与 Euler法形式上相似,但实际计算时却复杂得多。向前 Euler公式 是显式的,可直接求解。向后 Euler公式的右端含有ynt,因此是隐式公式,一般要用 迭代法求解,迭代公式通常为 ym1=y, +hf(x,, yu) y+=yn+hf(xn+1,y()(k=0,2,…) 2.2 Euler方法的误差估计 对于向前Euer公式(3)我们看到,当n=1,2,…时公式右端的yn都是近似的, 所以用它计算的yn+会有累积误差,分析累积误差比较复杂,这里先讨论比较简单的 所谓局部截断误差。 假定用(3)式时右端的yn没有误差,即yn=y(xn),那么由此算出 ym1=y(x,)hf(xm,(x)
-180- 需要说明的是，用不同的差商近似导数，将得到不同的计算公式。 （ii）用数值积分方法 将问题（1）的解表成积分形式，用数值积分方法离散化。例如，对微分方程两端 积分，得 ∫ + + − = = 1 ( ) ( ) ( , ( )) ( 0,1, ) 1 n n x x y xn y xn f x y x dx n L (4) 右边的积分用矩形公式或梯形公式计算。 （iii）Taylor 多项式近似 将函数 y(x) 在 n x 处展开，取一次 Taylor 多项式近似，则得 ( ) ( ) '( ) ( ) ( , ( )) n 1 n n n n n y x ≈ y x + hy x = y x + hf x y x + 再将 ( ) n y x 的近似值 n y 代入上式右端，所得结果作为 ( ) n+1 y x 的近似值 n+1 y ，得到离 散化的计算公式 ( , ) n 1 n n n y = y + hf x y + 以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的 计算公式。其中的 Taylor 展开法，不仅可以得到求数值解的公式，而且容易估计截断 误差。 §2 欧拉（Euler）方法 2.1 Euler 方法 Euler 方法就是用差分方程初值问题（3）的解来近似微分方程初值问题（1）的解， 即由公式（3）依次算出 ( ) n y x 的近似值 y (n = 1,2,L) n 。这组公式求问题（1）的数值 解称为向前 Euler 公式。 如果在微分方程离散化时，用向后差商代替导数，即 h y x y x y x n n n ( ) ( ) '( ) 1 1 − ≈ + + ， 则得计算公式 ⎩ ⎨ ⎧ = + = + + + = ( ) ( , ) ( 0,1, ) 0 1 1 1 y y a y y hf x y n n n n n L （5） 用这组公式求问题（1）的数值解称为向后 Euler 公式。 向后 Euler 法与 Euler 法形式上相似，但实际计算时却复杂得多。向前 Euler 公式 是显式的，可直接求解。向后 Euler 公式的右端含有 n+1 y ，因此是隐式公式，一般要用 迭代法求解，迭代公式通常为 ⎪⎩ ⎪ ⎨ ⎧ = + = = + + + + + + ( , ) ( 0,1,2, ) ( , ) ( ) 1 1 ( 1) 1 (0) 1 y y hf x y k L y y hf x y k n n n k n n n n n （6） 2.2 Euler 方法的误差估计 对于向前 Euler 公式（3）我们看到，当n = 1,2,L时公式右端的 n y 都是近似的， 所以用它计算的 n+1 y 会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的 所谓局部截断误差。 假定用（3）式时右端的 n y 没有误差，即 ( ) n n y = y x ，那么由此算出 ( ) ( , ( )) n 1 n n n y = y x + hf x y x + （7）
局部截断误差指的是,按(7)式计算由xn到xn+这一步的计算值yn+与精确值y(xn) 之差ν(xm+1)-yn+1。为了估计它,由 Taylor展开得到的精确值y(xn+)是 y(mD=y(x,)+hy'(x,)+y(x,)+O(h') (8) (7)、(8)两式相减(注意到y=f(x,y))得 h2 y(x-+ 2"(x,)+O(h)=O(h) (9) 即局部截断误差是h2阶的,而数值算法的精度定义为 若一种算法的局部截断误差为O(h),则称该算法具有p阶精度。 显然p越大,方法的精度越高。式(9)说明,向前 Euler方法是一阶方法,因此 它的精度不高。 §3改进的Euer方法 3.1梯形公式 利用数值积分方法将微分方程离散化时,若用梯形公式计算式(4)中之右端积分, ∫∫(x(x)k=U(xy(x)+f(x…y(x: 并用yn,yn代替y(xn),y(xn),则得计算公式 yu=yn+lf(n,y,)+f(n+l,ynD) 这就是求解初值问题(1)的梯形公式 直观上容易看出,用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法 梯形公式也是隐式格式,一般需用迭代法求解,迭代公式为 ym4=y,+hf(x,y) ym =yn+lf(u,y)+f(n, yuD) (10) (k=0,1,2,…) 由于函数∫(x,y)关于y满足 Lipschitz条件,容易看出 其中L为 Lipschitz常数。因此,当0<-<1时,迭代收敛。但这样做计算量较大。 如果实际计算时精度要求不太高,用公式(10)求解时,每步可以只迭代一次,由此导 出一种新的方法一改进 Euler法。 32改进 Euler法 按式(5)计算问题(1)的数值解时,如果每步只迭代一次,相当于将 Euler公式 与梯形公式结合使用:先用 Euler公式求yn的一个初步近似值yn1,称为预测值,然 后用梯形公式校正求得近似值yn1,即 181
-181- 局部截断误差指的是，按（7）式计算由 n x 到 n+1 x 这一步的计算值 n+1 y 与精确值 ( ) n+1 y x 之差 1 1 ( ) n+ − n+ y x y 。为了估计它，由 Taylor 展开得到的精确值 ( ) n+1 y x 是 ' '( ) ( ) 2 ( ) ( ) '( ) 3 2 1 y x O h h y x y x hy x n+ = n + n + n + （8） （7）、（8）两式相减（注意到 y'= f (x, y) ）得 ' '( ) ( ) ( ) 2 ( ) 3 2 2 1 1 y x O h O h h y x y n+ − n+ = n + ≈ （9） 即局部截断误差是 2 h 阶的，而数值算法的精度定义为： 若一种算法的局部截断误差为 ( ) p+1 O h ，则称该算法具有 p 阶精度。 显然 p 越大，方法的精度越高。式（9）说明，向前 Euler 方法是一阶方法，因此 它的精度不高。 §3 改进的 Euler 方法 3.1 梯形公式 利用数值积分方法将微分方程离散化时，若用梯形公式计算式(4)中之右端积分， 即 [ ( , ( )) ( , ( ))] 2 ( , ( )) 1 1 1 ∫ ≈ + + + + n n n n x x f x y x f x y x h f x y x dx n n 并用 1 , n n+ y y 代替 ( ), ( ) n n+1 y x y x ，则得计算公式 [ ( , ) ( , )] 2 n+1 = n + n n + n+1 n+1 f x y f x y h y y 这就是求解初值问题（1）的梯形公式。 直观上容易看出，用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。 梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = = + + = + + + + + + ( 0,1,2, ) [ ( , ) ( , )] 2 ( , ) ( ) 1 1 ( 1) 1 (0) 1 k L f x y f x y h y y y y hf x y k n n n n n k n n n n n (10) 由于函数 f (x, y) 关于 y 满足 Lipschitz 条件，容易看出 | | 2 | | ( 1) 1 ( ) 1 ( ) 1 ( 1) 1 − + + + + + − ≤ − k n k n k n k n y y hL y y 其中 L 为 Lipschitz 常数。因此，当 1 2 0 < < hL 时，迭代收敛。但这样做计算量较大。 如果实际计算时精度要求不太高，用公式（10）求解时，每步可以只迭代一次，由此导 出一种新的方法—改进 Euler 法。 3.2 改进 Euler 法 按式（5）计算问题（1）的数值解时，如果每步只迭代一次，相当于将 Euler 公式 与梯形公式结合使用：先用 Euler 公式求 n+1 y 的一个初步近似值 n+1 y ，称为预测值，然 后用梯形公式校正求得近似值 n+1 y ，即
ym=yn+的f(xn,yn) 预测 h (11) yn+[f(xn,yn)+∫(xn1,yn)]校正 式(11)称为由 Euler公式和梯形公式得到的预测一校正系统,也叫改进 Euler法 为便于编制程序上机,式(11)常改写成 yu+hf(xn, yn) hf(x, +h, y,) 改进Euer法是二阶方法。 §4龙格一库塔( Runge.-Kuta)方法 回到 Euler方法的基本思想一用差商代替导数一上来。实际上,按照微分中值定理 应有 y(xmD-y(x) (xn+h),0<b<1 注意到方程y=f(x,y)就有 y(x,+1)=y(xu)+hf(x,+ Ch, y(x,+Bh)) 不妨记K=f(xn+团mh,y(xn+团mh),称为区间[xn,xn+1]上的平均斜率。可见给出一种 斜率K,(13)式就对应地导出一种算法。 向前 Euler公式简单地取f(xn,yn)为K,精度自然很低。改进的 Euler公式可理 解为K取f(xn,yn),f(xn+1,丙n+1)的平均值,其中n1=yn+hf(xn,yn),这种处 理提高了精度。 如上分析启示我们,在区间[xn,xn1内多取几个点,将它们的斜率加权平均作为 K,就有可能构造出精度更高的计算公式。这就是龙格一库塔方法的基本思想。 41首先不妨在区间[xn,xn+1]内仍取2个点,仿照(13)式用以下形式试一下 yn+I=y+h(nk,+n,,) k,=f(n,y,) (14) k2=f(x, +ah, y,+ Bhk,),0<a,B< 其中λ,2α,B为待定系数,看看如何确定它们使(14)式的精度尽量高。为此我们 分析局部截断误差y(xn+1)-yn+1,因为yn=y(xn),所以(14)可以化为
-182- ⎪ ⎩ ⎪ ⎨ ⎧ = + + = + + + + + 校正 预测 [ ( , ) ( , ) ] 2 ( , ) 1 1 1 1 n n n n n n n n n n f x y f x y h y y y y hf x y （11） 式（11）称为由 Euler 公式和梯形公式得到的预测—校正系统，也叫改进 Euler 法。 为便于编制程序上机，式（11）常改写成 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = + = + + = + + ( ) 2 1 ( , ) ( , ) n 1 p q q n n p p n n n y y y y y hf x h y y y hf x y （12） 改进 Euler 法是二阶方法。 §4 龙格—库塔（Runge—Kutta）方法 回到 Euler 方法的基本思想—用差商代替导数—上来。实际上，按照微分中值定理 应有 '( ), 0 1 ( ) ( ) 1 = + < < + − y x θh θ h y x y x n n n 注意到方程 y'= f (x, y) 就有 ( ) ( ) ( , ( )) y xn+1 = y xn + hf xn +θh y xn +θh （13） 不妨记 K f (x h, y(x h)) = n +θ n +θ ，称为区间[ , ] n n+1 x x 上的平均斜率。可见给出一种 斜率 K ，（13）式就对应地导出一种算法。 向前 Euler 公式简单地取 ( , ) n n f x y 为 K ，精度自然很低。改进的 Euler 公式可理 解为 K 取 ( , ) n n f x y ， ( , ) n+1 n+1 f x y 的平均值，其中 ( , ) n 1 n n n y = y + hf x y + ，这种处 理提高了精度。 如上分析启示我们，在区间[ , ] n n+1 x x 内多取几个点，将它们的斜率加权平均作为 K ，就有可能构造出精度更高的计算公式。这就是龙格—库塔方法的基本思想。 4.1 首先不妨在区间[ , ] n n+1 x x 内仍取 2 个点，仿照（13）式用以下形式试一下 ⎪ ⎩ ⎪ ⎨ ⎧ = + + < < = + = + + ( , ), 0 , 1 ( , ) ( ) 2 1 1 1 1 1 2 2 α β α β λ λ k f x h y hk k f x y y y h k k n n n n n n （14） 其中 λ1,λ2 ,α,β 为待定系数，看看如何确定它们使（14）式的精度尽量高。为此我们 分析局部截断误差 1 1 ( ) n+ − n+ y x y ，因为 ( ) n n y = y x ，所以（14）可以化为
)+h(nk,+2k2) k,=f(n,y(xn))=y(x,) k,=f(,+ah,y(x,)+Bhk,) (15) =f(m,y(x,))+ahf(xn,y(x,)) + Bhk, (xmn, y(xn))+O(h) 其中k2在点(xn,y(xn)作了 Taylor展开。(15)式又可表为 n=y(xn)+(+A2)hy(x)+2mh2(1+2)+O(h2) 注意到 h V(m41=y(x,)+hy'(,)+y(x)+O(h') 中y=f,y=f1+∥y,可见为使误差y(xn+)-yn+=O(h2),只须令 11+2=1,2a 1 B (16) 2 a 待定系数满足(16)的(15)式称为2阶龙格一库塔公式。由于(16)式有4个未知数 而只有3个方程,所以解不唯一。不难发现,若令A1=2=,a=B=1,即为改 进的 Euler公式。可以证明,在[xn,xn+内只取2点的龙格一库塔公式精度最高为2 424阶龙格一库塔公式 要进一步提高精度,必须取更多的点,如取4点构造如下形式的公式 yu+1=y,+h(k,+12k2+1 k3+A4k) xn,yn k2=f(r+ah,y,+B,hk,) (17) Ch, y,+,hk,+B, hk2) h,y,+Bhk,+B,+Phk,) 其中待定系数λ1,a1,月共13个,经过与推导2阶龙格一库塔公式类似、但更复杂的计 算,得到使局部误差y(xn+)-yn1=O(h)的11个方程。取既满足这些方程、又较简 单的一组A,a1,B,可得
-183- ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ + + = + = + + = = + = + + ( , ( )) ( ) ( , ( )) ( , ( )) ( , ( ) ) ( , ( )) '( ) ( ) ( ) 2 1 2 1 1 1 1 1 2 2 hk f x y x O h f x y x hf x y x k f x h y x hk k f x y x y x y y x h k k y n n n n x n n n n n n n n n β α α β λ λ （15） 其中 2 k 在点( , ( )) n n x y x 作了 Taylor 展开。（15）式又可表为 ( ) ( ) '( ) ( ) ( ) 2 3 yn+1 = y xn + 1 + 2 hy xn + 2 h f x + ff y +O h α β λ λ λ α 注意到 ' '( ) ( ) 2 ( ) ( ) '( ) 3 2 1 y x O h h y x y x hy x n+ = n + n + n + 中 y'= f ， x y y' '= f + ff ，可见为使误差 ( ) ( ) 3 y xn+1 − yn+1 = O h ，只须令 1 λ1 + λ2 = ， 2 1 λ2α = ， =1 α β （16） 待定系数满足（16）的（15）式称为 2 阶龙格—库塔公式。由于（16）式有 4 个未知数 而只有 3 个方程，所以解不唯一。不难发现，若令 2 1 λ1 = λ2 = ，α = β = 1，即为改 进的 Euler 公式。可以证明，在[ , ] n n+1 x x 内只取 2 点的龙格—库塔公式精度最高为 2 阶。 4.2 4 阶龙格—库塔公式 要进一步提高精度，必须取更多的点，如取 4 点构造如下形式的公式： ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ = + + + + = + + + = + + = + = + + + + ( , ) ( , ) ( , ) ( , ) ( ) 4 3 4 1 5 2 6 3 3 2 2 1 3 2 2 1 1 1 1 1 1 1 2 2 3 3 4 4 k f x h y hk hk hk k f x h y hk hk k f x h y hk k f x y y y h k k k k n n n n n n n n n n α β β β α β β α β λ λ λ λ （17） 其中待定系数λi αi β i , , 共 13 个，经过与推导 2 阶龙格—库塔公式类似、但更复杂的计 算，得到使局部误差 ( ) ( ) 5 y xn+1 − yn+1 = O h 的 11 个方程。取既满足这些方程、又较简 单的一组λi αi β i , , ，可得