第3章 常微分方程的差分方法

第3章 常微分方程的数值解法§1 §2 §3 §4 §5 §6 §7 引 言 欧拉方法 龙格-库塔方法 阿达姆斯方法 算法的稳定性及收敛性 方程组及高阶方程的数值解法 边值问题的数值解法1

§1 引 言在工程和科学技术的实际问题中，常需要解常微 分方程。但常微分方程组中往往只有少数较简单和典 型的常微分方程(例如线性常系数常微分方程等)可 求出其解析解。对于变系数常微分方程的解析求解就 比较困难，而一般的非线性常微分方程就更不用说了。 在大多数情况下，常微分方程只能用近似法求解。这 种近似解法可分为两大类：一类是近似解析法，如级 数解法、逐次逼近法等；另一类则是数值解法，它给 出方程在一些离散点上的近似解。 在具体求解微分方程时，需要具备某种定解条件， 微分方程和定解条件合在一起组成定解问题。定解条2

件有两种：一种是给出积分曲线在初始点的状态，称为 初始条件，相应的定解问题称为初值问题 ；另一种是 给出积分曲线首尾两端的状态，称为边界条件 ，相应 的定解问题则称为边值问题。 例如，弹簧-质量系统的振动问题(图7-1),作一定的 简化后，可用一个二阶常微分方程

d x c x 0 2 dt m来描述。式中，x是质量m离平衡位置(0点)的距离；

2

t 是时间 ；c是弹簧常数。3

当弹簧在振动开始时刻t = t0 时的初始位置x( t0)= x0 和初速度

dx dt

x (t ) x 0 0 t t0

确定时，弹簧的振动规律x(t) 也就唯一确定。这就 是一个常微分方程的初值问题，可写成：2 d x c dt 2 m x 0 x(t 0 ) x0 (t t 0 ) x (t ) x 0 0

xm x c o

图7-1本章先从一阶常微分方程的初值问题：

dy f ( x, y ), x [a, b] (1.1) dx y (a) 出发进行讨论。5

由常微分方程的理论知，只要上式中的函数f(x,y)在区 域 G a x b, y 内连续，且关于y满足李普希

兹(Lipschitz)条件，即存在常数L(它与x, y无关)使

f ( x, y1 ) f ( x, y 2 ) L y1 y 2对 G a x b, y 内任意两个y1和 y2 都成立， 则方程的解必定存在且唯一。下面的分析均假定满足上述 条件。

初值问题(1.1)的数值解法，常采用差分方法，即把一个连续的初值问题离散化为一个差分方程来求解。即将 (1.1)离散化后，求找其解y

= y (x)在一系列离散节点6

a = x0 < x1 < … < xi < … < xn = b上的近似值y0, y1, …, yn。两相邻节点间的距离

hi = xi+1 - xi (i=0,1,2,…,n-1)称为步长。当

hi h (常值)时称为等步长，有

xi x0 ih, (i 1,2,... n)或

x x h, (i 0,1,2,...n 1)i 1 i7

因为初值问题

中的初始条件 y(a) y0 为已知，即 可利用已知的 y 0 来求出下一节点处 y ( x1 )的近似值 y 1

再从 y 1 来求出 y 2 , , 如此继续，直到求出 y n 为止。这 种用按节点的排列顺序一步一步地向前推进的方式求解的差分算法称为“步进式”或“递进式”算法。它

是初值问题数值解法的各种差分格式的共同特点。因此，只要能写出由前几步已知信息

y0 , y1 , , yi来计算的递推公式(即差分格式),即可完全表达该 种算法。8

§22.1 欧拉格式

欧拉方法

对于初值问题(1.1)式，先将其离散化，即把

[a,b] 作 n等分，得各离散节点

xi = a + ih (i=0,1,2, …, n-1)式中 h=(b-a)/n 设y = y (x) 为方程(1.1)的解，则 y = y (xi+1) 在(xi , yi) 点处的泰勒展开式为：

h2 y ( xi 1 ) y ( xi ) hy ( xi ) y ( i ) 2 (2.1) ( xi xi 1 )(i 0,1,2 , n 1)当 y' ' ( i ) 有界且 h 充分小时，可忽略高阶无穷小i

量

h y ' ' ( ) 将上式写成2

2

y ( xi 1 ) y ( xi ) hy' ( xi )或

y ( xi 1 ) y ( xi ) hf ( xi , y ( xi )) (i 0,1,2 , n 1)10

若将 y( xi 1 )和 y ( xi ) 的近似值分别记为 yi 1和 yi , 则得

yi 1 yi hf ( xi , yi )(i 0,1,2 , n 1) (2.2)这就是欧拉(Euler)公式，又称欧拉格式。利用它 可由已知的初值 y 0 出发，逐步算出 y1 , y2 , yn 。这 类形式的方程也称为差分方程。 当假定 yi 为准确，即在 yi y( xi ) 的前提下来估 计误差 y( xi 1) yi 1 ,这种截断误差称为局部截断误 差。由(2.1)和(2.2)可知，欧拉格式在节点x i 1 处的局部截断误差显然为：

R y( x ) yi i 1 2

i 1

h y' ' ( ), ( ), (i 0,1,2, , n 1) x x 2 (2.3) O(h ) i i i i 1 2

如果局部截断误差为O(h ),则称这种数值算法 的精度为P阶。故欧拉格式的精度为一阶。从几何意义 上来看欧拉格式，可如图7-2中所示。由方程(1.1) 知，其积分曲线 y = y (x) 上任意一点(x , y) 的切线斜 率dy/dx 都等于函数 f (x , y) 的值。从初值点P0(即 点(x0 , y0) )出发，作积分曲线y = y (x) 在P0 点上的 切线 P0 P1(其斜率为f (x0 , y0) ),与直线 x = x1 相交于点 P1(即点(x1 , y1) ),得到y1作为y(x1) 的

p 1

近似值，则有12

dy y y (x x ) dx1 0 1 0 x1

y 0 hf ( x0 , y 0 )x x0

y ( x1 ) y ( x0 ) x0 f ( x, y ( x)) dx相比较可知，这时是用切线段 P P 近似代替了曲线 0 1 段 P0 P1 ; P1 点近似代替了P '点； y1 近似代替了1

y(x1) ;hf (x0 , y0) 近似代替了

f ( x, y( x))dxx1 x0

递推继续从P1点出发，作一斜率为f (x1 , y1) 的直线

P P ,与直线x = x2 相交于P2 点(即点(x2 , y2) )1 2

y

P1P0

P21

P

P 2

Pi Pi 1Pi P i 1

Pn 1

Pn Pn

y y (x)

Pn 1

o

x0 x1 x2

xi xi 1 xn 1 xn图7-2

x

与

得到 y2 作为y(x2 ) 的近值；……如此继续，直到Pn 点。这样，得出一条折线P0 P1P2…Pi…Pn 近似代

替积分曲线P'0 P'1P'2…P'i…P'n 。当步数越多时，由于误差的积累，折线P0 P1P2…Pi…Pn可能会越来 越偏离真解 P'0 P'1P'2…P'i…P'n 曲线。

差分是微分的近似计算，所以欧拉格式也可用差商近似代替导数的离散方法得到。在节点xi 处有

y ( xi ) f ( xi , y( xi ))(i 0,1,2,..., n) (2.4)用向前差商

y ( xi 1 ) y ( xi ) h15

近似代替上式中的导数项y'(xi ) , 即

y( x ) y( x ) y ( x ) , (i 0,1,2,..., n 1(2.5) ) hi 1 i i

代入(2.4),可得：

y ( xi 1 ) y ( xi ) hf ( xi , y ( xi )) (i 0,1,2,..., n 1)y(xi+1 )和 y(xi )用其近似值 yi+1 和 yi 代入，则得

y y hf ( x , y )(i 0,1,2,..., n 1)i 1 i i i