常微分方程的差分方法

(i 0,1,2,...n,1) (2.9)
21
不管是显式欧拉格式〔2.2〕，还是隐式欧拉格式 〔2.6〕，它们都是单步格式或称为一步格式。因 为它们在计算yi+1时只用到前一步所得结果yi一个 信息；而格式〔2.8〕那么除了yi外，还需用到更 前一步所得信息yi-1，即需调用前两步的信息，因 此〔2.8〕称为两步欧拉格式，或称为中点欧拉格 式。
y(3)(i)
y(xi1)2hf(xi,y(xi))
(i 0,1,2,...n,1)
20
y(x)y ,(x)和 y(x)分别用其近似值代入，则得
i1
i
i1
yi1 yi1 2hf(xi,yi) (2.8) (i0,1,2,...n,1)
显然，其局部截断误差为
h3 R
y ( (3) )
i3
i
O(h3)
第章常微分方程的差 分方法
1
§1 引 言 在工程和科学技术的实际问题中，常需要解常微分方程。但常微分方程组中往往只有少数较简单和典型
的常微分方程〔例如线性常系数常微分方程等〕可求出其解析解。对于变系数常微分方程的解析求解就比 较困难，而一般的非线性常微分方程就更不用说了。在大多数情况下，常微分方程只能用近似法求解。这 种近似解法可分为两大类：一类是近似解析法，如级数解法、逐次逼近法等；另一类那么是数值解法，它 给出方程在一些离散点上的近似解。
在具体求解微分方程时，需要具备某种定解条件，微分方程和定解条件合在一起组成定解问题。定解 条
2
件有两种：一种是给出积分曲线在初始点的状态，称为 初始条件，相应的定解问题称为初值问题 ；另一种是 给出积分曲线首尾两端的状态，称为边界条件 ，相应 的定解问题那么称为边值问题。
例如，弹簧-质量系统的振动问题〔图7-1〕，作一定的 简化后，可用一个二阶常微分方程
预测 y y hf(x , y )
校正
yi1i1
y i
i
ii
h
f (x , y ) f (x ,
2
ii
i1
(i 0,1,2,...n, 1)
y i1
)〔2.11〕
26
可以证明，〔2.11〕格式的精度为二阶。它是一种一步 显式格式。〔2.11〕也可写成
y yhf(x,y)f(x,yh(fx,y))
i1
i2
ii
i1 i
ii
(i0,1,2,..n. ,1)
（2.12）
或将它写成以下平均化的形式：
差分是微分的近似计算，所以欧拉格式也可用
差商近似代替导数的离散方法得到。在节点xi 处有
y (x i) f(x i,y (x i)i) 0 (,1 ,2 ,.n ) .〔.2.4,〕
用向前差商 y(xi1)y(xi) h
15
近似代替上式中的导数项y'(xi ) , 即
y(xi)y(xi 1)h y(xi),(i0 ,1 ,2 ,.n . .1 ,)〔2.5〕
i 1
i
7
因为初值问题中的初始条件 y(a)y 为，即可利 0
用的 来求出y下一节点处 0
的近y似(x值) 再从 y来
1
1
求出y 1
如此y继, 续,，直到求出 为止。这y n种用按节 2
点的排列顺序一步一步地向前推进的方式求解的差分
算法称为“步进式〞或“递进式〞算法。它是初值问
题数值解法的各种差分格式的共同特点。因此，只要
y(xi ) hf (xi1, y(xi1))
(i 0,1,2,...,n 1)
y(xi+1 )和 y(xi )用其近似值 yi+1 和 yi 代入，那么得
yi1 yi hf(xi1,yi1) (i0,1,2,...n,1) 〔2.6〕
18
其局部截断误差为：
Ri
h2 2
y(i
)
x i1 i xi (2.7)
9
y(x )y(x)hy(x)h2 y()
i1
i
i
2
Байду номын сангаас
i
(2.1)
(xx )i(0,1,2 ,n1)
i
i1
当 y''( ) 有界且 h 充分小时，可忽略高阶无穷小
i
量 h2 y' ' ( ) 将上式写成
2
i
y(xi1)y(xi)hy'(xi)
或 y(xi1)y(xi)hf(xi,y(xi))
(i0,1,2,n1)
i1
i2
ii
i1 i1 〔2.10〕
(i0,1,2,...n,1)
使误差相互抵消，提高误差阶数，从而提高算法的精 度。事实上，格式(2.10)的局部截断误差为O(h3)，即 其精度为二阶。(2.10)称为梯形格式，是一种隐式格 式。从几何角度上来看，梯形格式是取[xi,xi+1]的两 端点的平均斜率。从图7-3可见，用梯形格式得到的P 点比用显式欧拉格式得到的A点和用隐式欧拉格式得到 的B点都要合理。
y,y。, 这y类形式
12
n
当假定
y
i
为准确，即在
y i
y(xi)
的前提下来估
x 计误差 y( )y
i1
i1
，这种截断误差称为局部截断误差。
x 由〔2.1〕和〔2.2〕可知，欧拉格式在节点 处的局 i1
部截断误差显然为：
11
R y y(x )
i
i1
i1
h2 y''()
2
i
(,xi i
x ) i1
19
再如，用中心差商表示的三点数值微分公式
y(x)
1[y(x
)y(x
h2 )]
y ( (3) )
i
2h
i1
i1
6
i
(i0,1,2,..n.,1)
代入
y ( x i) f( x i,y ( x i)( ) i 0 ,1 ,2 ,.n ) ..,
可得
y(xi1)y(xi1)2hf(xi,
y(xi)
)h3 3
10
假设将y(x ) i 1
和y( x ) i
的近似值分别记为y i 1
和y i
，
那么得
y i 1 y i h ( x i,f y i ) i ( 0 , 1 , 2 ,n 1 ) (2.2)
这就是欧拉〔Euler〕公式，又称欧拉格式。利用它可
由的初值 出发y，逐步算出 0
的方程也称为差分方程。
xx0
y(x) 1
y(x ) x1
0
x0
f
(x,
y(x))dx
相比较可知，这时是用切线段 P P 近似代替了曲线
段
P
P
； P1 点近似代替了 P
01
'点； y1 近似代替了
0
y(x1)
；1 hf
(x0
,
y0) 近似代替了1
x1
x0
f
(x,
y(x))dx
递推继续从P1点出发，作一斜率为f (x1 , y1) 的直线
d 2x c x
dt 2 m
x(t ) x
0
0
0 (t
t) 0
x(t ) x
0
0
4
x m
x co
图7-1
本章先从一阶常微分方程的初值问题：
dy f (x, y),x[a,b]
dx
〔1.1〕
y(a)
出发进行讨论。
5
由常微分方程的理论知，只要上式中的函数f(x,y)在区
域 G a x b , y 内连续，且关于y满足李普希 兹(Lipschitz)条件，即存在常数L(它与x, y无关)使
(x0 , y0) 〕出发，作积分曲线y = y (x) 在P0 点上的切
线
〔其斜率为f (x0 , y0) 〕，与直线 x = x1 相交
于点 PP10〔P即1 点(x1 , y1) 〕，得到y1作为y(x1) 的近似
值，那么有
12
dy
y y (x x )
y hf(x,y )
1
0
d x 1
0
0
00
比较〔2.3〕,〔2.7〕和〔2.9〕可知，两 步欧拉格式比显式或隐式欧拉格式具有更高的精度， 因为它的局部截断误差是O(h3)。
由〔2.3〕和(2.7)可见，显式欧拉格式与隐式 欧拉格式的局部截断误差的符号正好相反，因此可 以设想取〔2.2〕和(2.6)的平均，即两式相加除以 2，得
22
y yhfx,yfx ,y
25
2.2 改进的欧拉格式
显式欧拉格式计算工作量比较小，但精度低。梯形格式 虽提高了精度，但为隐式格式，需用迭代法求解，计算工 作量大。综合这两种格式可得到改进的欧拉格式。
为预先测用值欧，拉它格的式精〔度2不.2高〕，求接出着一用个梯初形步格的式近〔似2值.10y〕i 1对，它称
校正一次，即迭代一次，求得yi+1 ，称为校正值。这种预 测-校正方法称为改进的欧拉格式：
O(h2) (i0,1,2,...n,1)
(2.6)称为向后欧拉格式，又称为隐式欧拉格式。因 为在此式的右边也包含未知的yi+1，所以是yi+1的一 个函数方程，故称它为隐式格式。〔2.2〕的右边那 么没有未知的yi+1，因此是一种显式格式。隐式格式 的计算当然比显式格式要困难得多，一般情况下，只 能用迭代法求解，计算量比较大。
能写出由前几步信息
y,y, ,y
01
i
来计算的递推公式〔即差分格式〕，即可完全表达该 种算法。
8
§2 欧拉方法 2.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) 点处的泰勒展开式为：
6
a = x0 < x1 < … < xi < … < xn = b
上的近似值y0, y1, …, yn。
两相邻节点间的距离
hi = xi+1 - xi (i=0,1,2,…,n-1)
称为步长。当 h h 〔常值〕时称为等步长，有 i
xi x0i,h (i1 ,2,.n .).
或
x x h ,(i 0 ,1 ,2 ,.n .1 .)
欧拉格式〔2.2〕的程序图见图7-4
23
y A P
y y(x)
B
O
x
i
图7-3
xx i 1
24
输入 a、b、n、a
h (b a) / n xi a yi a
i0,1 ,2,..n .,
xi1xi h yi1yi hf(xi,yi)
输出 , x i 1 y i 1
i n
i n stop
图7-4
d2x dt2
c m
x
0
来描述。式中，x是质量m离平衡位置〔0点〕的距离； t 是时间 ；c是弹簧常数。
3
当弹簧在振动开始时刻t = t0 时的初始位置x( t0)= x0
和初速度
dx
x(t ) x
dt
0
0
t t0
确定时，弹簧的振动规律x(t) 也就唯一确定。这就
是一个常微分方程的初值问题，可写成：
代入〔2.4〕,可得：
y(xi1)y(xi)h(fxi,y(xi)) (i0,1,2,..n.,1)
y(xi+1 )和 y(xi )用其近似值 yi+1 和 yi 代入，那么得
y y h (x f ,y )i (0 ,1 ,2 ,.n . 1 .),
i 1
i
ii
16
此即〔2.2〕〔欧拉格式〕。
显然，欧拉格式具有递推性，在计算yi+1时只要用到前一 步所得结果yi一个信息就够了，因此是一种单步格式或称一步 格式。其中〔2.5〕是一个数值微分公式。故用其他数值微分 公式也可导出略异于〔2.2〕的其他形式的算式来。例如，用 向后差商表示的两点数值微分公式
P P ，与直线x = x2 相交于P2 点〔即点(x2 , y2) 〕 12
13
y
P0
P1
P1
P2 P2
Pi Pi1 Pn 1 Pn
y y(x)
Pi
Pi
1
Pn
1
Pn
o
x0 x1 x2 x i xi1 x n1 xn
x
图7-2 与
14
得到 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 曲线。
f(x ,y 1 )f(x ,y 2) L y 1 y 2
对 G a x b , y 内任意两个y1和 y2 都成立，
那么方程的解必定存在且唯一。下面的分析均假定满足上 述条件。
初值问题〔1.1〕的数值解法，常采用差分方法，即把 一个连续的初值问题离散化为一个差分方程来求解。即将 〔1.1〕离散化后，求找其解y = y (x)在一系列离散节点
y(x )1y(x )y(x)hy()
i1
h
i1
i
2i
(i0,1, ,n1)
代入
y ( x) f( x,y ( x)( i ) 0 ,, 1 ,2 ,.n .1 ) .,
i 1
i 1
i 1
17
可得
y(xi1) y(xi ) hf (xi1, y(xi1))
h2 2
y(i ),xi1
i
xi
(,i0,1,2,,n1)
O(h2)
(2.3)
如果局部截断误差为O(hp1)，那么称这种数值算
法的精度为P阶。故欧拉格式的精度为一阶。从几何意
义上来看欧拉格式，可如图7-2中所示。由方程〔1.1〕
知，其积分曲线 y = y (x) 上任意一点(x , y) 的切线斜
率dy/dx 都等于函数 f (x , y) 的值。从初值点P0〔即点