常微分方程数值解法

第七章 常微分方程数值解法

常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。因此，有必要探讨常微分方程初值问题的数值解法。本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。

第一节 欧拉法

求解常微分方程初值问题

?dy??f(x,y) ?dx??y(x0)?y0 （1）

的数值解，就是寻求准确解y(x)在一系列离散节点x0?x1?x2???xn?? 上的近似值 y0,y1,y2,?,yn,?

?yn?称为问题的数值解，数值解所满足的离散方程统称为差分格式，hi称为步长，实用中常取定步长。

?xi?xi?1显然，只有当初值问题(1)的解存在且唯一时，使用数值解法才有意义，这一前提条件由下 面定理保证。

定理 设函数f?x,y?在区域D:a?x?b,???y???

上连续，且在区域D内满足李普希兹(Lipschitz)条件，即存在正数L，使得对于R内任意两点?x,y1?与?x,y2?,恒有的解y?x?存在并且唯一。 一、欧拉法(欧拉折线法)

若将函数

yxf?x,y1??f?x,y2??Ly1?y2则初值问题(1)

?在点xhn?处的导数y?xn?用两点式代替， 即

y??xn??y(xn?1)?y(xn)，再用yn近似地代替y?xn?,则初值问题(1)变为

?yn?1?yn?hf(xn?yn) (2) ?y?y(x),n?0,1,2,?0?0

1

(2)式就是著名的欧拉(Euler)公式。以上方法称 为欧 拉法或欧拉折线法。欧拉公式有明显的几何意义。从几何上看，求解初值问题(1)就是xy平面上求一条通过点?x0,y0?的曲线y?y?x?，并使曲线上任意一点?x,y?处的切线斜率为

f?x,y?。 欧拉公式的几何意义就是从点

P0x0,y0?出发作一斜率为f?x0,y0?的

直线交直线x?x1于点P1?x1,y1?，P1点的纵坐标y1就是y?x1?的近似值；再从点P1作一斜率为f?x1,y1?的直线x?x2交直线于点P2?x2,y2?, P2点的纵坐标y2就是的近似值y?x2?；如此继续进行，得一条折线P0P1P2?。该折线就是解y?y?x?的近似图形，如图7-1。

图7-1

欧拉法的其它几种解释：

1．

假设y?x?在xn附近展开成泰勒级数

y?xn?1??y?xn??hy??xn???y?xn??hf?xn,y?xn???hh222y???xn???2y???xn???

取h的线性部分，并用yn作为y?xn?的近似值，得 yn?1?yn?hf?xn,yn?

dy2． 对方程

dx?f?x,y?两边从xn到xn?1积分，得

y(xn?1)?y?xn???xn?1xnf(x,y(x))dx (3)

2

用 矩形公式计算上式右侧积分，即 ?xn?1xnf(x,y(x))dx??xn?1xnf?x,y?x??dx

并用yn作为的近似值y?xn?，得yn?1?yn?hf?xn,yn? 故欧拉法也称为矩形法。 二、改进的欧拉法(梯形法)

欧拉法形式简单，低，特别当曲线y=y(x)计算方便，但精度比较的曲率较大时，欧拉法的效果更差。为了达到较高精度的计算公式，对欧拉法进行改进，用梯形公式计算(3)式 右侧积分，即

?xn?1xnf(x,y(x))dx?h2?f?xn,y?xn???f?xn?1,y?xn?1???

并用yn作为y(xn)的近似值，得到改进 的欧拉公式

yn?1?yn?h2?f?xn,yn??f?xn?1,yn?1?? (4)

上述方法称为改进的欧拉法，也称为梯形法。

不难发现，欧拉公式 是 关于yn?1的显式，即只要已知yn,经过一次计算便可得yn?1的值，而改进的 欧拉公式是以yn?1的隐式方程给出，不能直接得到yn?1。隐式方程(4)通常用 迭代法求解，而迭代过程的实质是逐步显式化。

先用欧拉 公式

yn?1?yn?hf?xn,yn??0?

给出yn?1的迭代初值，然后 再用改进的欧拉公式(4)进行迭代，即有

(0)?yn?1?yn?hf(xn,yn)?h?(k?1)(k)y?y?f(xn,yn)?f(xn?1,yn?1)?n?1n (5) 2?k?0,1,2,?????迭代过程进行到连续两次迭代结果之差的绝对值小于给定的精度?即

yn?1?yn?1???k?1?k

为止，这时取

yn?1?yn?1?k?1?

然后再转入下一步计算。

3

???y?下面讨论是否收敛；若收敛，它的极限是否满足(4)式。

kn?1假设f(x,y)满足李普希兹条件 f?x,y1??f?x,y2??L?y1?y2? 则

yn?1?yn?1???k?1??k?h2f?xn?1,yn?1??f?xn?1,yn?1?k??k?1??hL2yn?1?yn?12?k??k?1??hL?????2?yn?1?yn?12?k?1?k?2?hL??????2??yn?1?yn?1?1??0?

?hL??k→∞时，有??2?k由此可见，只要

kn?1hL2?1(这里只要步长h足够小即可)，当

?0，

???y?所以收敛。又因为f(x ,y)对y连续，当k→∞时，对等式

yn?1?k?1??yn?h2?f?xn,yn??fxn?1,yn?1??k???

两端取极限， 得

yn?1?yn?h2?f?xn,yn??f?xn?1,yn?1??kn?1

n?1???y?因 此，只要步长h足够小，就可保证收敛到满足(4)式的y。

三、预估一校正法

改进的欧拉公式在实际计算时要进行多次 迭代，因而计算量较大。在实用上，对于改进的欧拉公式(5)只迭代一次，即先用欧拉公式算出yn?1的预估值y

(0)n?1，再用改进的欧拉公式(4)进行一次迭代得到 校正值yn?1，即

(0)?yn?1?yn?hf(xn,yn)?h?yn?1?yn??f?xn,yn??f?xn?1,yn?1??,n?0,1,2,??2? (6)

预估—校正公式也常写成下列形式：

4

11?y?y?k?k2n1?n?122?k1?hf(xn,yn),n?0,1,2,???k?hf?x?h,y?k?nn1?2? (7)

四、公式的截断误差

定义 若某种微分方程数值解 公式的截断误差是O(h

k?1)，则称这种方法是

k阶方法。为了简化分析，在进行误差分析时，我们假设前一步的结果是准确的，即在 yn=y(xn)的前提下，考虑用yn?1作为y(xn?1)的近似值而产生的截断误差，这种误差称为局部截断误差。由泰勒公式

y?xn?1??y?xn?h??y?xn??by??xn??h22!y???xn???

对于欧拉公式，有

yn?1?yn?hf?xn,yn??y?xn??hy??xn? 于是

y?xn?1??yn?1?Oh2??

2则欧拉公式的截断误差为O(h),所以 欧拉法是一阶方法。

对于预估—校正公式，有

k1?hf?xn,yn??hy??xn?k2?hf?xn?h,yn?k1??hf?xn?h,y?xn??k1??hf?xn,y?xn???hf?hf?xn,y?xn???h2?x?xn,y?xn???k1fy?xn,y?xn?????xn?f?x,y?xn???y??xn?fy?xn,y?xn?????

而

y??x??f?x,y?x??y???x??fx?x,y?x???y??x??fy?x,y?x?? 于是

k2?hy??xn??hy???xn???2

因此

yn?1?yn?12k1?312k2?y?xn??hy??xn??h22y???xn???

3所以y(xn?1)- yn?1= O(h)则预估—校正公式的截断误差为O(h),也即预估—

5

校正法是二阶方法。

可以证明，改进的欧拉公式与预估—校正公式的截断误差相同，均为O(h)。这里略去证明。

例 1：在区间［0,1］上以h=0.1为步长，分别用欧拉法与预估— 校正法求初值

2x?dy?y??问题?dxy的数值解

??y(0)?13解：该方程为贝努利方程，其精确解为y?1?2x欧拉公式的具体形式为

?2xnyn?1?yn?h?y??nyn?????

预估—校正公式的具体形式为

?11?yn?1?yn?k1?k2?22?2xnk1?hf(yn?)?yn???2?xn?h????k?hfy?k??21?nyn?k1??? ?取步长h?0.1,x0?0,y0?1,计算结果见下表：表7-1

表7-1

n xn y 欧拉法 预估-矫正法 11.0959091.1840971.2662011.3433601.4164021.4859561.5525141.6164751.678166准确解 11.0954451.1832161.2649111.3416411.4142141.4832401.5491931.6124521.67332001234567891000.10.20.30.40.50.60.70.80.911.11.1918181.2774381.3582131.4351331.5089661.5803381.6497831.717779 1 1.784771 1.737867 1.732051

6

近似解与准确解比较，欧拉法的结果大致只有两位有效数字，而预估—校正法的结果则 有3位有效数字。

第二节 龙格—库塔法

前面讨论的欧拉法与改进的欧拉法都是一步 法，即计算yn?1时，只用到前一步值。龙格—库塔(Runge-Kut ta)法(简称为R-K方法)是一类高精度的一步法，这类方法与泰勒级数法有着密 切的关系。 一、泰勒级数法

设有初值问题

?dy??f(x,y)?dx??y(x0)?y0

由 泰勒展开式

y?xn?1??y?xn?h??y?xn??hy??xn??hkh22!y???xn????k!y?k??xn??O?hk?1? （1）

若令

yn?1?y?xn??hy??xn??h22!???hkk!y?k??xn? (2)

则

y?xn?1??yn?1?Oh?k?1?

即公式(2)为k阶方法。

从理论上讲，只要解y(x)有任意阶导数，泰勒展开方法就可以构造任意阶求yn?1的公 式，但由于计算这些导数是非常复杂的。如

y??x??f?x,y?x???fy???x??fx?y?fy?fx?ffyy???fxx?2ffxy?fxfy?ffy?f22fyy??

所以这种方法实际上不能用来解初值问题。 二、龙格 —库塔方法(R-K方法)

R-K方法不是通过求导数的方法构造近似公式，而是通过计算 不同点上的函数值，并对这些函数值作线性组合，构造近似公式，再把近似公式与解的泰勒展开

7

式进行比较，使前面的若干项相同，从而使近似公式达到一定的阶数。我们先分析欧拉法与预估—校正法。对于欧拉法

?yn?1?yn?k1??k1?hf(xn,yn)

每步计算f的值一次，其截断误差为O(h)。对于预估—校正法。 

11?y?y?k?k2n1?n?122?k1?hf?xn,yn???k?hf?x?h,y?k?nn1?2?

3每步计算f的值两次，其截断误差为O(h)。

2下面对预估—校正法进行改进，将该公式写成更一般的形式 

yn?1?yn?R1k1?R2k2k1?hf?xn,yn?k2?hf?xn?ah,yn?bh? (3)

其中R1,R2,a,b为待定常数。选择这 些常数的原则是在yn?y(xn)的前提下，使y(xn?1)?yn?1)的阶尽量 高。为此，作泰勒展开

k2?hf(xn?ah,yn?bk)?h[f?ahf2x?bk1fy??]x?hf?h(af?bffy)??

其中

f,fx,fy,?都是在(xn,yn)处的函数值。将k1,k2代入yn?1得

yn?1?yn?(R1?R2)hf?h(aR2fx?bR2ffy)???y(xn)?h(R1?R2)y?(xn)?h(aR2fx?bR2ffx??

22与泰勒展开式(2)进行比较，要使得

??R1?R2???aR2???bR?2??y?xn?1??yn?1?oh??只要 四个参数满足

3?11212 (4)

若R1?R2?12,a?b?1,即得预估—校正公式。满足(4)式的R1,R2,a,b可以

8

有各种不同的取法，但不管如何取法，都要计算两次f的值(即 计算f在两个不同点的函数值)，截断误差都是Oh??。满足条件(4)的一族公式(3)统称为 二

3阶龙格—库塔公式。人们容易想到，如果不增加计算函数值的次数，能否适当地选择这 四个参数，使近似公式的局部截断误差的阶再提高，比如达到

Oh??。为此，把k42多展 开 一项，有

?bffyk2?hf?h2?afx??h32?a2fxx?2abffxy?bf22fyy?Oh???4 （4）所

以

yn?1?y?xn??h?R1?R2?y??xn??hh32?aR2fx?bR2ffy??2?a2R2fx?2abR2ffxy?bR2f22fyy?Oh???4

而y?xn?1?在xn的泰勒展开式

y?xn?1??y?xn??hy??xn???y?xn??hf?h2h22y???xn??h3h36y????xn??Oh2??422?fx?ffy??36?fxx?2ffxy?f2fyy?fxfy?ffy?Oh???4

由于

y?xn?1?展开式的h项中

fxfy?ffy是无法通过选择参数R1,R2,a,b来消

去的，所以不论四个参 数如何选择，都不能使局部截断误差y?xn?1??yn?1达到O?h?。要想提 高近似公式的阶，只能继续增加计算f的值的次数。如

4果每步计算三次f的值，可将公式写 成下列形式：

?yn?1?yn?R1k1?R2k2?R3k3??k1?hf?xn,yn? (5) ???k?hfx?ah,y?bkn2n211?2?k3?hf?xn?a3h,yn?b31k1,b32k2??与二阶龙格—库塔公式的讨论方法类似，要使

y?xn?1??yn?1?Oh??，只

4 9

?R1?R2?R3?1?a?b21?2?a3?b31?b32?1aR?aR?需8个参数满足? 33?222?122?a2R2?a3R3?3?1?a2b32R3??6?方程组包含6个方程，8个未知量，其解不唯一。满足条件(6)的一族公式(5)统称为三阶龙格—库 塔公式。一个比较简单的三阶龙格—龙塔公式是

141?y?y?R?R?R3n12?n?1666?k?hf?x,y??1nn?11???k2?hf?xn?h,yn?k1?22?????k3?hf?xn?h,yn?k1?2k2?

截断误差为O?h?的四阶龙格—库塔公式 是常用的公式，每步都要计算四次

5f的值。它的一般形式是

?yn?1?yn?R1k1?R2k2R3k3R4k4?k?hf(xn,yn)?1? ?k2?h(xn?a2h,yn?b21k1) (6)

?k?hf?x?ah,y?bk?bk?n3n311322?3??k4?hf?xn?a4h,yn?b41k1?b42k2?b43k3?

(6)式中13个待定常数需满足下列11个方程的方程组

10

表7-6

j bj 0 1 1 -1/2 2 -1/12 3 -1/24 当k=0时，(7)式为  yn?1?yn?hfn?1 当k=1时，(7)式为梯形公式 当 k=3时，(7)式为 即

yn?1?yn?h24(9fn?1?19fn?5fn?1?fn?2?1)

yn?1?yn?h(fn?1?12?fn?1?112?fn?1?2124?fn?1)3

(8)

类似于亚当斯显式公式，可以导出隐式公式(8.3.8)的局部截断误差为

y(xn?1)?yn?1??19720hy5(5)(xn)???O(h)5

所以隐式公式(8)是四阶方法。 7.3.3亚当斯预估一校正公式

 把亚当斯显式公式与隐式公式联立使用，构造亚当斯预估一校正公式。以四阶亚当斯方法为例有下列公式：

yn?1?yn?(0)h24

(55fn?59fn?1?37fn?2?9fn?3)h24

(9)

yn?1(k?1)?yn?(9f(xn?1,yn?1?19fn?5fn?1?fn?2)(k)迭代过程进行到连续 两次迭代结果之差的绝对值小于给定的精度ε为止。

(k?1)这时取yn?1?yn?1，然后转入下一步计算。只迭代1次的公式

yn?1?yn?(0)h24h(55fn?59fn?1?37fn?2?9fn?3)

(.10)

yn?1(k?1)?yn?24(9f(xn?1,yn?1?19fn?5fn?1?fn?2)(k)称为亚当斯预估一校正公式。第1式称为预估公式，第2式称为校正公式。其局部截断误差为

16

y(xn?1)?yn?1??19720hy5(5)(xn)?O(h)5

下面我们具体讨论误差估计问题。设 由(10)式求得准确解 y(xn+1)的近似解为yn+1，由于

y(xn?1)?yn?1?(0)25172019720hy5(5)(xn)

y(xn?1)?yn?1??hy5(5)(xn)上面二式相减得

y(xn?1)?yn?1?(0)270720 所以

hy5(5)(xn)

y(xn?1)?yn?1yn?1?yn?1(0)??19270

于是有

y(xn?1)?yn?1??19270(yn?1?yn?1)(0)

若

?19270(yn?1?yn?1?ε

(0)

则可持续求解yn+2,?，否则,应缩小步长h重新计算。

与同阶的龙格—库塔方法相比较，亚当斯方法计算量小，公式简单，程序易于实现。但是由于亚当斯方法在计算yn+1时，不仅用到前面一步的信息，而且还要用到更前面几步的信息，因此它不能自动开始。开始的前几个值必须借助于一步法获得。

例：分别用四阶亚当斯显式公式与预估一校正公式求解初值问题

dx dy?y?2xy

y(0)?1 0?x?1,h?0.1

解：用第二节例1，按标准四阶龙格—库塔方法求得的结果y1,y2,y3作为开始值，然后用显式公式(6)与预估一校正方法(10) 进行计算，计算结果如表7-7所示。

17

表7-7

xn 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

yn R-K方法 0 1.09544553 1.18321675 1.26491223 显式方法 1.34155176 1.41404642 1.48301819 1.54891888 1.61211634 1.67291704 1.73156976 预估计-校正法 准确解 1.341641136 1.41421383 1.48323983 1.54919338 1.61245154 1.67332000 1.73205072 1 1.09544512 1.18321596 1.26491106 1.34164079 1.41421356 1.48323970 1.54919334 1.61235155 1.67332005 1.73205081 近似解与准确解进行比较，显式方法的结果有4位有效数字，而预估一校正法的结果则有7位有效数字。

第四节 线性多步法

线性多步法的一般公式为

yn?1?0yn?1yn?1=+

+?+?kyn?k+h(??1fn?1??0fn?...??kfn?k)

(1)

当β-1 =0时，公式(1)是显式的；当β-1 ≠0时 ，公式(1)是隐式的。“线性”是指yn?1是yn,yn?1, ? , yn?k, fn?1,fn,?,fn?k的线性组

合，“多步”是指计算yn+1时 ，不仅用到前一步的结果yn,而且用到更前面几步的结果。所有形如(8.4.1)式的公式都 可以利用泰勒展开的方法构造出来。作为例子，我们推导形如

yn?1?0yn?1yn?1?2yn?2=

+

+

+h(

?0fn??1fn?1??2fn?2??3fn?3)

(2)

的数值计算公式。

18

假设 yi=y(xi), i=n, n-1, n-2 fi=f(xi,yi)i=n,n-1,n-2,n-3 将上述各项在xn处泰勒展开 Yn-i=y(xn-ih)=y(xn)-ih

fn?iy(xn)'+

(ih)2!2y(xn)-

''(ih)3!3y(xn)+?, i=1,2

'''=

y(xn?i)y(xn)''=

y(xn?ih)'

2 =-

ihy?(xn)+

(ih)2!y???(xn)-

(ih)3!3y(4)(xn)+?, i=1,2,3

代入(2)式 得

yn?1=(?0??1??2)y(xn)+(??0?2?2??0??1??2??3)hy?(xn)+

(?1?4?2?2?1?4?2?6?3)(

h22!y??(xn)+ h3??1?8?2?3?1?12?2?27?3)

3!y???(xn)+

4(

?1?16?2?4?1?32?2?108?3)

h4!y(4)(xn)+... (3)

将y(xn?1)在xn处泰勒展开

y(xn?1)?y(xn)?hy?(xn)?h22!y??(xn)?h33!y???(xn)?h33!y???(xn)?... (4)

比较(3)式与(4)式，让h的同次幂的系数相等，得  ?0??1??2?1

??1?2?2??0??1??2??3?1

?1?4?2?2?1?4?2?6?3?1 (5) ??1?8?2?3?1?12?2?27?3?1 ?1?16?2?4?1?32?2?108?3?1 ??

如果(5)的前7个方程求出未知数?0,?1,?2,?0,?1,?2,?3,则其局部截断误差为

O(h).当然，也可以从前k

7个方程解出未知数(k≤7)，其局部截断误差为

19

O(h).比如从前5 个方程解出

k7个未知数，因此，有2个自由未知数，若令

?1??2?0,则从(5)的前5 个方程解得 

?0?1,?0?5524,?1??5924,?2?3724,?3??924

对应的公式是

yn?1?yn?h24(55fn?59fn?1?37fn?2?9fn?3)

它正好是四阶阿达姆斯显式公式。

类似地， 若设计算公式为

yn?1= ?0yn+?1yn?1+?2yn?2+h(??1fn?1??0fn??1fn?1??2fn?2)

为使这类公式为四阶方法，系数应满足 ?0??1??2?1

??1?2?2???1??0??1??2?1

?1?4?2?2??1?2?1?4?2?1 (7) ??1?8?2?3??1?8?1?12?2?1 ?1?16?2?4??1?4?1?32?2?1 特别地，取?1??2?0，则从(7)式解得 

?0?1,??1?924,?0??h241924,?1??524,?2?124

对应的公式是

yn?1?yn?(9fn?1?19fn?5fn?1?fn?2)

它正好是四 阶亚当斯隐式公式。

第五节 方程组与高阶方程的数值解法

前面几节介绍了一阶微分方程初值问题的各种数值解法，这些解法同样适用 于一阶微分方程组与高阶方程。为了避免书写的复杂，下面仅就两个未知函数的方程组和二阶方程为例说明各种方法的计算公式。 7.5.1一阶方程组

设有一阶微分方程组初值问题

20

dydxdz?f(x,y,z),y(x0)?y0

(1)

?g(x,y,z),z(x0)?z0dx

1、欧拉法的计算公式 yn?1?yn?hf(xn,yn,zn) zn?1?zn?hg(xn,yn,zn)

y0?y(x0),z0?z0(x0), n?0,1,2,? (2) 2、改进的欧拉法的计算公式

对n=0,1,2,?，计算

y(0)n?1?yn?hf(xn,yn,zn) (0)

zn?1?zn?hg(xn,yn,zn)

y(k?1)(k)(k)n?1?yn?h2[f(xn,yn,zn)?f(xn?1,yn?1,zn?1)] (3)

z(k?1)?hn?12[f(x(k)(k)

?znn,yn,zn)?f(xn?1,yn?1,zn?1)]

k?0,1,2,?

当连续两次迭代结果之差的绝对值小于给定的精度ε,即 y(k?1)n?1?y(k)n?1?ε z(k?1)n?1?z(k)n?1?ε

时，取y(k?1)n?1?yn?1,z(k?1)n?1?zn?1，然后转入下一步运算。

3、四阶标准龙格—库塔公式

yn?1?yn?16(k1?2k2?2k3?k4) z1n?1?zn?6(m1?2m2?2m3?m4)

k1?hf(xn,yn,zn) m1?hg(xn,yn,zn)

k12?hf(xhn?2,yn?k12,zn?m2) (4)

mk12?hg(xhn?2,yn?2,zn?m12)

21

k3?hf(xhn?2,yk2n?2,zn?m22)

mk2m23?hg(xn?h2,yn?2,zn?2)

k4?hf(xn?h,yn?k3,zn?m3) m4?hg(xn?h,yn?k3,zn?m3)

n?0,1,2,?

4、四阶亚当斯显式公式 h

yn?1?yn?24(55fn?59fn?1?37fn?2?9fn?3)

zhn?1?zn?24(55gn?59gn?1?37gn?2?9gn?3) (5)

n?4,5,6,?  其中 

 fn?i?f(xn?i,yn?,zn?i)

gn?i?g(xn?i,yn?,zn?i) i?0,1,2,3 需要注意，前几步的值由其它单步法求得。 7.5.2二阶方程

设有二阶微分方程初值问题

y???g(x,y,y?) y(x0)?y0,y?(x0)?y?0 (6)

令z=y′，则(6)化为一阶微分方程初值问题

dydx?z

dzdx?g(x,y,z) (7)

y(x0)?y0,z(x0)?y?0

(7)式中

dydx?z?f(x,y,z),于是应用四阶龙格—库塔公式(4),有k1?hzn?y?n

22

令

k2?h(zn?m12m22)?hy?n?hm12hm23

2k3?h(zn?)?hy?n?

k4?h(zn?m3)?hy?n?hm则

yn?1?yn?16(k1?2k2?2k3?k4)?yn?hy?n?16(m1?m2?m3)

又因为

zn?y?n?1,zn?1?y?n?1

zn?1?zn?16(m1?2m2?2m3?m4)

所以

y?n?1?zy?n?16(m1?2m2?2m3?m4)

于是得初值问题(6)的计算公式为

yn?1?yn?hy?n?1616(m1?m2?m3)

y?n?1?y?n?(m1?2m2?2m3?m4)m1?hg(xn,yn,y?n)

m2?hg(xn?h2h2,yn?121212hy?n,y?n?141212m1)

12m2)m3?hg(xn?,yn?hy?n?hm1,y?n?

m4?hg(xn?h,yn?hy?n?hm2,y?n?m3)

第六节 边值问题的数值解法

在具体求解微分方程时，需要附加某种定解条件。微分方程与定解条件一起组成定解问题。对于高阶微分方程，定解条件通常有两种给法，第一种是给出积分 曲线在初始时刻的性态，称这类条件为初始条件，相对应的定解问题称为初值问题；第二种是给出积分曲线在首末两端的性态，称这类条件为边值条件，相

23

对应的定解问题为边值问题 。 本节以二阶微分方程

y″=f(x,y,y′),a≤x≤b

为例讨论边值问题，其边值条件 可分为三类：

第一边值条件

y(a)=a, y(b)=β

第二边值条件

 y′(a)=a, y′(b)= β

第三边值条件

y′(a)-a0y(a)=a , y′(b)- β0y(b)=β

其中 a0≥0,β0≥0,a0 +β0>0.

本节介绍解线性方程边值问题的差分方法以及 适用于非线性方程边值问题的试射法

。

7.6.1解线性方程边值问题的差分方法

二阶线性方程的一般形式为

 y″+p(x)y′+q(x)y=f(x), a≤x≤b (1)

首先将区间［a,b］进行等距划分，分点

 xi?a?ih, i=0,1,2,?,n 其中

h?b?an

一般地称x0,?a与xn＝b为边界点，称x1,x2,?,xn?1为内部节点。

其次，在各个xi节点上，将y′,y″用差商近似表示。这里要求有相同阶数的截断误差，以保证精度协调。对于内部节点，二阶导数用二阶中心差商表示，得

yi?1?2yi?yi?1h2?y??(xi)?O(h), i=1,2,?,n-1

2一 阶导数用一阶中心差商表示，得

24

yi?1?yi?12h?y?(xi)?O(h), i=1,2,?,n-1

2 假设yi?y(xi),，则

y??(xi)?yi?1?2yi?yi?1h2,

y?(xi)?yi?1?yi?12h,于是得方程(1)的差分方程

yi?1?2yi?yi?1h2+

piyi?1?yi?1qiy?fi+

2h

i=1,2,?,n-1 (2)

其中pi?p(xi),qi?q(xi),fi?f(xi)．

将(2)整理，可以写成下列形式

aiyi?1?biyi?ciyi?1?di  i=1,2,?,n-1 (3)

其中 

ai?1?1212hpi2bi??2?hqici?1?2hpi di?hfi

(3)是含有n+1个未知数yi (i=0,1,2,?,n)的线性方程组，方程的个数为n-1。要使方程组(3)有唯一解，还需要有两上边值条件补充两个方程。对于第一边值条件，直接就得到两个方程

y0?a yn??

于是得到第一边值问题的差分方程组

?y1??d1?a1a2?????y2d2????b1c1???????an?2?a2b2? (4) y??d?n?2??n?2??y????n?1??dn?1?cn?1??c2 25

这个方程组是三对角方程组，可以利用追赶法求解。

 对于第二及第三 边值条件，由于条件中包含了导数，所以边值条件也必须用差商来近似表示。为使截断误差达到O(h)，需要用到牛顿等距插值公式，得

y??1h(?y0?1212?y0)?22?3y0?4y1?y22hyn?2?4yn?1?3yn2h

??(?yn?yn?yn)?2

将这两个近似公式代入边值条件中，就得到两个方程，再与(3)联立，就可得到对应的差分方程组， 用追赶法解出yi(i=0,1,?,n)。

对于边值问题的收敛性，即考虑当h→0时，差分 方程组的解是否收敛于微分方程的准确解。以第一边值问题为例，给出如下结论，不加证明。 定理 若ai?0,ci?0,bi?ai?ci, i=1,2,?,n-1，则差分方程组(.4)存在唯一解。且当h→0时，(8.6.4)的解收敛于第一边值问题的准确解。 例 用差分法解边值问题 y″-y=x

y(0)=0, y(1)=1, 0≤x≤ 1, h=0,1 解： h?

?2?y1??0.1?10??????2y20.2?10????????????????2?y8??0.8?10??y???0.991??? ?9????

110,节点xi?i10,边值问题的差分方程可写成下列形式

?

用追赶法解方程组，结果如表7-8．

表7-8

i 1 差分解yi 0.070489 准确解y(xi) 0.070467 26

2 3 4 5 6 7 8 9 7.6.2试射法

0.142684 0.218305 0.299109 0.386904 0.483568 0.591068 0.711479 0.847005 0.142641 0.218244 0.299033 0.386819 0.483480 0.590985 0.711411 0.846963 设二阶微分方程第一边值问题  y″=f(x,y,y′)

(5)

y(a)=a,y(b)=β

试射法的基本思想是把边值问题化为初值问题来解。具体作法是通过反复调整初始时刻的斜率y′( a)的值m,使得初值问题

y″=f(x,y,y′) y(a)=a,y′(a)=m

的解满足另一个边值条件y(b)=β,也就是从初值问题(6)的经过点(a,a),而且有不同斜率的积分曲线中，去寻找一条通过(b，β)点的曲线。

首先凭经验或按照实际存在的运动规律选取m的两个预测值m1,m2，再分别按照两个斜率求解相应的初值问题(6)，可以得到y(b)的两个结果β1,,β2。如果β1,,β2都不满足给定的精度，就用线性插值的方法校正m1,与m2 ，得到新的斜率值

m3?m1?m2?m1(6)

?2??1(???1)

然后再按斜率值m3计算初值问题(6)，又得新的结果y(b)=β3。继续这一过程，直到计算结果y(b)与β相 当接近为止。

值得注意的是，用线性插值的依据是不足的。如果有更好的插值公式可利用

27

的话，则可能使测试的次数有效地减少。

28

本文来源：https://www.bwwdw.com/article/x4z.html