第九章+常微分方程初值问题数值解法

第九章 常微分方程初值问题

的数值解法

在自然科学的工程技术的许多领域中，常会遇到常微分方程初值问题，但这种问题大多数情况下不存在初等形式的解析解，只能用近似方法来求解。近似解法主要有两类：一类叫做近似解析方法，它能给出解得近似表达式，例如熟知的级数解法和逐次逼近法等；另一类近似解法称为数值解法，它可以给出解在一些离散点上的近似值。数值方法便于电子计算机求解微分方程。本章讨论初值问题常用的数值解法，并介绍有关的基本理论。

假设给定一阶常微分方程的初值问题：

(,),()dy

f x y a x b dx y a a

⎧=≤≤⎪⎨⎪=⎩ （9.1）（9.2） 其中f 为x, y 的已知函数，α为给定的初始值。

由微分方程的理论可知，如果 f (x,y )在区域a x b ≤≤，

y -∞<<+∞内连续，且关于

y 满足李普希兹(Lipschitz)条件，即存

在常数L ，使

|(,)(,)||f x y f x y L y y -≤-

对所有的a x b ≤≤及任何y ，y 均成立，则初值问题（9.1），（9.2）有连续可微的解y (x )存在且唯一。所谓初值问题的数值解，则是问题的解y (x )在一系列点

01212...N N a x x x x x b --=<<<<<=

处的值()n y x 的近似值n y (n =0,1,…,N )。这里相邻两个节点之间的距离1n n n h x x -=-通常称为步长，通常将步长n h 取为常数h 。

§9.1 欧拉方法与改进的欧拉方法

欧拉(Euler)方法是最简单的数值方法，由于它的精确程度较差，已不常用于实际计算。但构造这个方法的基本原理，对于构造一般的数值方法具有普遍意义，因此首先对它进行讨论。

9.1-1 欧拉方法的构造

将方程（9.1）中点n x 处的导数()n y x '用差商近似地表示为

1()()

()n n n y x y x y x h

+-'≈

即在该点有近似等式

1()()

(,())n n n n y x y x f x y x h

+-≈

用近似值n y 代替()n y x ，由上式可以导出其近似值满足的差分方程：

1(,()),

0,1,...,1n n n n y y hf x y x n N +=+=- （9.3）

（9.3）式称为欧拉方法的计算公式或称为欧拉公式。当初始值

0y α=给定时，利用欧拉公式就可以逐次计算出初值问题的数值

解01,,...,N y y y 。

9.1-2 后退的欧拉公式

如果在点1n x +处用向后差商近似(9.1)式中的导数

1()()

()n n n y x y x y x h

+-'≈

即有 111()()(,())n n n n y x y x hf x y x +++≈+

再用近似值i y 代替()i y x (,1)i n n =+，则可导出近似值满足的差分方程

111(,)n n n n y y hf x y +++=+ （9.6）

（9.6）式称为后退的欧拉公式。这是一个两端含有未知数1n y +的方程，这样的方法称为隐式方法。相应地，称右端不隐含1n y +的欧拉公式（9.3）是一个显式方法。隐式方法求1n y +时需解方程，当(,)f x y 关于y 为非线性时，就需要解非线性方程，通常用迭代法求解。

可以证明，隐式的后退欧拉公式（9.6）的局部截断误差为

23,()()2

n h

n h R y x O h ''=-+ （9.7）

事实上，设()n n y y x =时，由（9.6）式有

11111(,)

()(,)

n n n n n n n y y hf x y y x hf x y +++++=+=+ （9.8）

由中值定理得

1111111(,)

(,())(,)[()]

n n n n y n n n f x y f x y x f x y y x η+++++++'=+- （9.9）

将(9.9)式代入（9.8）式，并利用微分方程（9.1），得

11111()()

(,)[()]

n n n y n n n y y x hy x hf x y y x η+++++'=+'+- （9.10）

由微分方程解的泰勒展式（9.4）与（9.10）式相减，得

11

2

11311()[()()]()(,)

2[()]()

n n n n n y n n n y x y h h y x y x y x hf x y y x O h η++++++-'''''=-+--+

以31()()()()n n n y x y x hy x O h +''''=++代入上式，得

11

23

111()()(,)[()]()

2

n n n y n n n y x y h y x hf x y x y O h η+++++-'''=-+-+ 移项整理得

11

2

31()1[()()]1(,)2

n n n y n y x y h y x O h hf x η+++-''=-+'- 注意当h 充分小，使1y hf '<时，有

111

1(,)...1(,)y n y n hf x hf x ηη++'=++'-

即得 2311

()()()2

n n n h y x y y x O h ++''-=-+

9.1-3 改进的欧拉方法

用1n y +与1n y +分别表示用欧拉公式与后退的欧拉公式求得的数值解，由式（9.5）与（9.7）可知，欧拉方法与后退的欧拉方法的局部截断误差的主部只是符号相反，于是显然有

31111

()()()2

n n n y x y y O h +++-+=

由此可以构造如下的计算公式

1

11[(,)(,)]2

n n n n n n h

y y f x y f x y +++=++

(0,1,2,...1)n N =-

（9.11）

（9.11）式称为梯形方法，其局部截断误差为3()O h ，因此它是一个二阶方法，梯形方法也是隐式的，在实际应用中时常与显式的欧拉公式联合使用，构成如下的计算格式：