第七章 常微分方程模型的数值解法
第七章 常微分方程数值解法简介
微分方程在科学和工程技术中有很广泛的应用。许多实际问题的数学模型都可以用微分方程来描述，归结为常微分方程的定解问题；很多偏微分方程问题，也可以化为常微分方程问题来近似求解，但是求出所需的解绝非易事。实际上，除了极特殊情形外，人们不可能求出微分方程的解析解，只能用各种近似方法得到满足一定精度的近似解。在常微分方程中已经熟悉了级数解法和Picard 逐步逼近法，这些方法可以给出解的近似表达式，称为近似解析方法。另一类方法只给出解在一些离散点上的值，称为数值方法。后一类方法应用范围更广，特别适合用计算机计算，本章主要介绍常用的常微分方程数值解法。
7.1实际问题的微分方程模型
函数是事物的内部联系在数量方面的反映,如何寻找变量之间的函数关系,在实际应用中具有重要意义。在许多实际问题中，往往不能直接找出变量之间的函数关系,但是有时却容易找出变量的改变量之间的关系，从而建立描述问题的微分方程模型。
例7.1.1 将初始温度00150u C =的一碗汤放置于环境温度a u 保持为024C 的桌上，10分钟后测得汤的温度为0100C 。如果汤的温度低于055C 才可以喝，试问再过20分钟后这碗汤能喝了吗？
解：为了解决这一问题，需要了解有关热力学的一些基本规律。热量总是从温度高的物体向温度低的物体传导的；在一定的温度范围内，一个物体的温度变化速度与这个物体的温度和其所在的介质温度的差值成正比。
设物体在t 时刻的温度为()u u t =，从t t t →+?温度从()()u t u t t →+?，注意到热量总是从温度高的物体向温度低的物体传导，因而0a u u >，所以温度差
a u u -恒正，又因物体将随时间而逐渐冷却；则温度的改变量为:
()()(())a u u t t u t k u t t u t
?=+?-=-+?-?
两边除以t ?，并令0t ?→得温度变化速度为：
()a du k u u dt
=--
这里0k >是比例常数。从而得出描述物体冷却过程的微分方程模型为：
0()(0)a du
k u u dt u u ?=--???=?
（7.1.1）
容易求出这个一阶微分方程初值问题的解为：
0()kt
a a u u u u e
-=+- （7.1.2）
根据所给问题的条件知当10t =时，01100u u C ==，得到：
1010()k
a a u u u u e
-=+-
将00150u C =，024a u C =带入得到：
0111ln ln 1.660.0511010
a a
u u k u u -=
≈≈-
从而得这碗汤的温度随时间变化的函数关系为： 0.0510.0510()()24126t
t
a a u u t u u u e
e
--==+-=+ （7.1.3）
于是，将30t =带入计算得，再过20分钟后汤的温度就是0251u C ≈。说明再过20分钟后这碗汤能喝了。
不过并不是所有的微分方程模型都可求出解析解，从而得出实际问题的解释。例如看似简单的微分方程
2
2
dy x y
dx
=+，自德国数学家Wilhelm von Leibniz 提出100
多年后才被法国数学家Joseph Liouville 证明它没有解析解，只能借助于数值的方法求数值解。
例7.1.2 某地区发现一种有免疫性的传染病，为了控制疫情扩散对该地区人群进行隔离处理。为了分析受感染人数的变化规律，需要建立描述传染病传播过程的数学模型。
解：设该地区的总人数为常数N ，任意t 时刻病人、健康人和病人治愈后移出感染系统的移出者的比例分别为(),(),()i t s t r t ，病人的日接触率λ，日治愈率μ，则容易得出从t t t →+?时刻，病人和健康人的改变量为：
[()()]()()()[()()]()()N i t t i t N s t i t t N i t t
N s t t s t N s t i t t λμλ+?-=?-???
+?-=-??
每个方程两边除以t ?，并令0t ?→，化简后得：
00(0),(0)d i s i i d t
d s s i d t i i s s λμλ?
=-??
?
=-?
?
==???
（7.1.4） 其中， ()()()1s t i t r t ++=，对任意的t 。
（7.1.4）就是描述病人和健康人的比例()i t 和()s t 随时间变化的微分方程模
型，这是一个微分方程组的初值问题。但是这一初值问题的解析解是无法求出的，因此不能直接利用()i t 和()s t 的解析式来分析和解决问题。
在数学建模课程中学到的大量数学模型都是微分方程形式给出的，各类微分方程本身和它们的解所具有的特性已在常微分方程及数学物理方程中得以解释。虽然求解微分方程有许多解析方法，但解析方法只能够求解一些特殊类型的方程，在实际应用中人们更关心的是某些特定的自变量在某一个定义范围内的一系列离散点上的近似值。这样一组近似解称为微分方程在该范围内的数值解，寻找微分方程数值解的过程称为微分方程的数值解法。
7.2简单的数值方法与基本概念
7.2.1常微分方程初值问题
设(,)f x y 在区域:,G a x b y ≤≤<∞上连续，求()y y x =满足
0(,),()dy
f x y a x b dx y a y
?=<≤?
?
?=?
（7.2.1） 其中0y 是已知常数,这就是一阶常微分方程的初值问题。为使问题（7.2.1）的解存在、唯一且连续依赖初值0y ，即初值问题（7.2.1）适定，还必须对右端项(,)f x y 加适当限制，通常要求(,)f x y 关于y 满足是已知函数,且满足Lipschitz 条件：存在常数L ，使
1212
(,)(,)f x y f x y L y y -≤- （7.2.2）
对所有[,]x a b ∈和12,(,)y y ∈-∞+∞成立。本章总假定f 满足此条件。
7.2.2 Euler 法及改进的Euler 法 7.2.2.1 Euler 方法的导出与几何意义
最简单的数值解法是Euler 法。将区间[0,]T 作N 等分，小区间的长度/h T N =称为步长，点列i x a ih =+，(0,1,...,)i N =称为节点，0x a =。由已知初值00()y x y =，可算出()y x 在0x x =的导数'00000()(,())(,)y x f x y x f x y ==。下面用三种方法导出Euler 法。本章用()i y x 表示函数()y x 在i x 点的精确值、i y 表示()i y x 的近似值。
（1）幂级数展开法
利用Taylor 展式
2
3
'
''
'''
100000000
()()()()()()
2
6
(,)h
h
y x y x h y x hy x y x y y hf x y R ξ=+=++
+
=++ （7.2.3）
其中01(,)x x ξ∈，并略去二阶小量0R ，得1000(,)y y hf x y =+，1y 就是1()y x 的近似值。利用1y 又可算出2y ，如此下去可算出()y x 在所有节点上的值，一般递推公式为：
1(,)n n n n y y hf x y +=+，0,1,...,1n N =- （7.2.4）
这就是Euler 法。
Euler 法有明显的几何意义。实际上，初值问题（7.2.1）的解是x ，y 平面上过点00(,)x y 的一条积分曲线。按Euler 法，过初始点00(,)x y 作经过此点的积分曲线的切线（斜率为00(,)f x y ），沿切线取点11(,)x y （1y 按（7.2.4）计算）作为11(,())x y x 的近似。然后过11(,)x y 作经过此点的积分曲线的切线（斜率为
11(,)f x y ），沿切线取点22(,)x y （2y 按（7.2.4）计算）作为22(,())
x y x 的近似。
如此下去，即得一以(,)n n x y 为顶点的折线，这就是用Euler 法得到的近似积分曲线（如图7.1.1所示）。从几何上看，h 越小，此折线逼近积分曲线越好，因此也称Euler 法为Euler 折线法。
图7.1.1 Euler 近似积分曲线
（2）数值微分法 利用向前差商近似导数
1()()
()n n n y x y x y x h
+-'≈
1()()()(,)n n n n n n y x y x hy x y h f x y +'≈+=+
从而得出Euler 法的一般递推公式：
1(,)n n n n y y hf x y +=+，0,1,...,1n N =-
（3）数值积分法
将初值问题（7.2.1）写成等价的积分形式：
0()()(,())x x y x y x f t y t dt =+
?
取1x x =时得
10
10()()(,())x x y x y x f t y t dt =+
?
（7.2.5）
用左矩形公式近似右端积分，并用1y 替代1()y x ，即得1000(,)y y hf x y =+，从而也可得出Euler 法的一般递推公式为：
1(,)n n n n y y hf x y +=+，0,1,...,1n N =-。
7.2.2.2 改进的Euler 方法
由Euler 方法的数值积分导出法可知只要给出（7.2.5）式右端定积分的一种近似计算方法，就可得出初值问题（7.2.1）的一种数值求解方法。
如果用右矩形公式近似（7.2.5）式右端积分，则可得1011(,)y y hf x y =+，从而也可得出一般递推公式为：
111(,)n n n n y y hf x y +++=+，0,1,...,1n N =- （7.2.6）
称（7.2.6）为后退Euler 法。
如果用梯形公式近似（7.2.5）式右端定积分，则可得
100011[(,)(,)]2
h y y f x y f x y =+
+
从而得出一般递推公式为：
111[(,)(,)]2
n n n n n n h y y f x y f x y +++=+
+，0,1,...,1n N =- （7.2.7）
称（7.2.7）为改进的Euler 法，显然改进的Euler 法比Euler 法精度更高。后退Euler 法和改进的Euler 法，由于未知数1n y +同时出现在等式的两边，故称为隐式；如果未知数1n y +由已知量直接计算即不只出现在等式右端，则称为显式。对于隐式算法每步计算需要解关于1n y +的方程，而这样的方程往往是非线性的，通常将
初值取为[0]
1
n n y y +=用迭代法求解，一般只需迭代几步即可收敛。一般先用显式公式计算一个初值，再用隐式公式迭代求解。
如果先用显式Euler 公式作预测，算出1(,)n n n n y y hf x y +=+，再将1n y +代入隐式梯形公式的右边作校正，得到111[(,)(,)]2
n n n n n n h y y f x y f x y +++=++。
从而可得
()11(,),(,)2
n n n n n n
n n h y y f x y f
x y h f x y ++=+
++????，0,1,...,1n N =-。
这种方法称为预估-校正法，同时可以看到它是显示格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。
如果在区间11[,]n n x x -+上对初值问题（7.2.1）的方程两边积分，则有
11
11()()(,())n n x n n x y x y x f x y x dx +-+--=
?
并用中矩形求积公式近似对右端的定积分，则得出一般递推公式为：
112(,)n n n n y y hf x y +-=+，0,1,...,1n N =- （7.2.8）
称（7.2.8）为Euler 中点公式。
如果计算1()n y x +的近似值1n y +时只用到前一节点的值n y ，则从初值0y 出发可逐一计算出以后各节点的值；这样的方法称为单步法。而Euler 中点公式法计算
1n y +时需要用到前两个节点的值n y 和1n y -，称这样的方法为双步法或二步法；计
算1n y +时需要用到前面多个节点值的方法称为多步法。多步法附加初值才能逐一计算出以后各节点的值。
7.2.3 截断误差与算法精度的阶
从Euler 法的几何意义上得知，由Euler 法所得的折线明显偏离了积分曲线，可见此方法非常粗糙即误差太大。现在分析一下求解初值问题（7.2.1）的数值方法误差的来源。为使问题简化，不考虑因计算机字长限制引起的摄入误差。
在假设()n n y y x =，即第n 步计算是精确的前提下，称第1n +步计算1()n y x +的截断误差11()n n n R y x y ++=-为局部截断误差。若某算法的局部截断误差为1()p O h +即为1p h +的同阶无穷小，则称该算法有p 阶精度。若局部截断误差
1
2
11()(,)()P p n n n n n R y x y x y h
O h
ψ++++=-=+，
则称1(,)P n n x y h ψ+为误差主项，称(,)n n x y ψ为误差主项系数。 （1）Euler 法的局部截断误差
由Euler 法的一般递推公式和1()n y x +的Taylor 展式得
112
3
2
2
3
3
()()[(,)]
()()[()()][()]2!
3!
()()()()
2!
3!
2
n n n n n n n n n n n n n n n n R y x y y x h y hf x y y x y x y x y x h h h y hy x y x y x h
h h y x O h ++=-=+-+'''''''=+++
+-+'''''''=+
+=
+
所以Euler 法的局部截断误差为2()O h ，即Euler 法为1阶精度算法，其误差主项为
2
()2
n h
y x ''。
（2）后退Euler 法的局部截断误差
同理，由后退Euler 法的一般递推公式、1()n y x +和1'()n y x +的Taylor 展式得
11112
3
2
3
2
3
2
()()[(,)]
()()[()()][()]2!3!()()[()()]
2!
3!
()[()()]2!()(2
n n n n n n n n n n n n n n n n n n n n n n R y x y y x h y hf x y y x y x y x y x h h h y hy x h y x y x y x y x h h h y x y y x h y x h h y x y x h ++++=-=+-+'''''''=++++-++''''''=++
+
+''''''-++++'''''=-
-
2
3
3
)()()
3
2
n n h
h y x O h ''-=-
+
所以后退Euler 法的局部截断误差也是2()O h ，即后退Euler 法也是1阶精度算法，其误差主项为2
()2
n h
y x ''-
。
（3）改进Euler 法的局部截断误差
由改进Euler 法的一般递推公式可得其局部截断误差为
111112
3
2
3
()()[((,)(,))]
2
()['()'()]2
[()()()()]['()'()]
2!
3!
2
()()()[()(2!
3!
2
n n n n n n n n n n n n n n n n n n n n n n n n h R y x y y x h y f x y f x y h y x h y y x y x h
h
h y x hy x y x y x y y x y x h h
h
h hy x y x y x y x y +++++=-=+-+
+=+--
+''''''=++++--
++''''''''=++
+-+ 2
3
4
()()())]
2!
1()()
12
n n n n y x x y x h h y x h O h '''''++
+'''=-
+ 所以改进Euler 法的局部截断误差为3()O h ，即改进Euler 法是2阶精度算法，其误差主项为3
1()12
n y x h
'''-
。
（4）整体截断误差
当然人们更关心的是近似解n y 的误差，即()n n n y x y ε=-，称为整体截断误差。由1()n y x +的Taylor 展式与Euler 法的一般递推公式相减得：
2
3
1112
3
()()()()()[(,)]2!
3!
[(,()(,))]()()2!3!
[(,()(,))]n n n n n n n n n n n n n n n n n n n n n n n
h h y x y y x hy x y x y x y hf x y h
h
h f x y x f x y y x y x h f x y x f x y R εεε+++''''''=-=++
++-+'''''=+-+
+
+=+-+
记max n n
R R =，因(,)f x y 关于y 满足Lipschitz 条件，所以
1(1)n n n n Lh R Lh R εεεε+≤++=++
以此递推，得
2
121
00
0(1)(1)(1)(1)(1)
(1)[(1)1]
n n n n n j
j n
n
Lh R Lh Lh R R
Lh R Lh R Lh Lh Lh
εεεεε---=≤++≤++++≤≤+++=++
+-∑
注意到0n x x nh b =+≤，()/n b a h ≤-，于是
()
()
0(1)L b a L b a n R e
e
L h
εε--≤+
-，1,2,...,n N =。
右端依赖于初始误差0ε和局部截断误差的界R 。对于Euler 法，可取2R Ch =（C
是与n 无关的常数）。若00ε=即00()y y x =无误差，则
()
()
(1)(1)L b a L b a n R C e
e
h L h
L
ε--≤
-=
-
所以()n O h ε=，比局部截断误差低一阶。用同样方法可以证明改进Euler 法的整体截断误差的阶为2()O h ，也比局部截断误差低一阶。 （5）Euler 算法的稳定性
在实际计算中，由于测量误差、舍入误差等因素的影响，初值0y 往往不能精确给出，其误差将依次传递下去。如果传递误差能够被控制，精确来说，传递误差连续依赖于初始误差，则称算法稳定；否则就说不稳定。显然不稳定的算法是不能用的。下面仅考察Euler 法的稳定性。
设从初值0u 和0v 算出的节点值分别为{}n u 和{}n v ，则满足
111(,)n n n n u u hf x u ---=+，111(,)n n n n v v hf x v ---=+，1,2,...,n N =
两式相减，并令n n n e u v =-，则
11111[(,)(,)]n n n n n n e e h f x u f x v -----=+-
从而
111
()
00
(1)(1)n n n n n
L b a e e Lh e Lh e Lh e e
e ----≤+=+≤≤+≤ ，（因nh b a ≤-）。
这说明n e 连续依赖初始误差0e ，即Euler 法稳定。同样可证改进的Euler 法也稳定。 例7.2.1取步长0.02h =,分别用Euler 法、后退Euler 法和中点法求解初值问题：
'
0.9(x)=-,00.112(0)1y y x x y ?≤≤?+?
?=?
。 解 因为步长0.02h =，所以各节点0.02,0,1,2,3,4,5.i x i i == 因为01,y =利用Euler 法的计算公式
10.9,12i i i i
hy y y x +=-
+
可取123450.9820,0.9650,0.9489,0.9337,0.9192y y y y y ===== 利用后退Euler 法的计算公式：
111
0.9,12i i i i hy y y x +++=-
+
可解得1i y +的显示表达：
11
,0.9112i i i y y h x ++=
+
+
于是
123450.9830,0.9669,0.9516,0.9370,0.9232.y y y y y =====
按中点法的计算公式
11 1.8,12i i i i
hy y y x +-=-
+
需要知道两个初值。在此，我们利用后退欧拉法计算的结果10.9830,y =再依次计算得
23450.9660,0.9508,0.9354,0.9218y y y y ====
而该初值问题的解析解是0.45(12),y x -=+用它计算各节点的函数值可得
123450.9825055161,0.9659603719,0.9502806582,0.9353925462,0.9212307783.
y y y y y =====
把上述三种方法计算的结果同准确对照,可以看出它们确实都是准确值的近似值,只是误差不一样。Euler 法的误差偏小,后退欧拉Euler 法偏大,中点法最小。
7.3 线性多步法
用Euler 法计算节点0n x x nh =+的近似值n y 时只用到前一节点的值1n y -，是线性的单步法。为了提高解的精度，需要构造线性多步法，其一般形式为
k
k
j
n j j n j j j y h f α
β++===∑∑ （7.3.1）
其中(,)n j n j n j f f x y +++=，j α和j β是常数，且0k α≠，0α和0β不同时为0。按（7.3.1）计算n k y +时要用到前面k 各节点的值11,,,n n n k y y y ++- ，因此（7.3.1）是多步法或称为k -步法。又因为方程（7.3.1）关于,n j n j y f ++是线性的，所以称为线性多步法。显然，若0k β=，则线性多步法（7.3.1）是显式
的；若0k β≠，则线性多步法（7.3.1）是隐式的。用线性多步法进行计算时，除需要给定0y 外，还要附加初值121,,,k y y y - ，这可用其它方法计算。由于多步法每计算一步用到的信息更多，因此希望构造出精度更高的算法。
7.3.1 数值积分法
将微分方程'()(,)y x f x y =在1[,]n n x x +上积分得
11()()(,())n n
x n n x y x y x f x y x dx ++=+
?
（7.3.2）
适当取1k +各节点，作被积函数(,())f x y x 的k 次Lagrange 插值多项式
,()n k L x ，并用,()n k L x 近似代替(,())f x y x ，就可得到形如（7.3.1）的线性多步
法。插值节点的不同取法就可得出不同的线性多步法。
（1）Adams 外插法
取1,,,n n n k x x x -- 为节点，构造(,())f x y x 的Lagrange 插值多项式
,()n k L x ，则
,,(,())()()n k n k f x y x L x r x =+ （7.3.3）
其中,()n k r x 是插值余项。将（7.3.3）代入（7.3.2），得
111,,()()()()n n n
n
x x n n n k n k x x y x y x L x dx r x dx +++=+
+
?
?
（7.3.4）
舍去余项
1,,()n n
x n k n k x R r x dx +=
?
（7.3.5）
并用j y 代替()j y x ，即得
11,()n n
x n n n k x y y L x dx ++=+
?
（7.3.6）
由此可见，Adams 方法的局部截断误差为1,,()n n
x n k n k x R r x dx +=
?
。
下面给出Adams 方法（7.3.6）的具体形式。假设前1(0)k k +>个节点处的函
数值已知，即()(,1,2,,)j y x j n n n n k =---L 的近似值j y 已算出，从而函数值
(),(,1,2,,)j j
j f f
x
y j n n n n k ==---L 也已知。这样就可以利用1k +组数据点：
(()),,n n n x f x y ，(())111,,n n n x f x y ---，…，(()),,n r n r n r x f x y ---
构造被积函数()(),f x y x 的k 次插值多项式
()()(),0
,k
n k n j
n j j j L x f x
y l x --==
∑
其中()(0,1,,)j l x j r =L 是拉格朗日插值基函数。从而可得
()()110,i i
k
x n n n j n j j x j y y f x y l x dx ++--==+
∑
?
（7.3.6）
记(),n j n j n j f f x y ---= ，()11
,i i
x kj j x l x dx h
α+=
?
则有
10
k
n n rj n j j y y h f α+-==+∑. （7.3.7）
这就是求解初值问题的Adams 显式公式，1n y +是关于(,1,2,,)j y j n n n n k =---L 的线性表达式，所以它是线性1k +步法。
在上述Adams 显式公式的推导中，选用了,n x 1,n x -…, n k x -作为1k +插值节点，但构造的k 次插值多项式(),n k L x 是代替区间[]1n n x x +上的未知函数
()(),f
x y x ，因此属于“外插”，称为Adams 外插法，也称为Adams-Bashorth 法。
显然当0k =时Adams 外插法就是Euler 法。
（2）Adams 内插法
如果将Adams 外插法推导过程中的节点改为1,n x +,n x 1,n x -…, 1n k x -+。这时，公式（7.3.7）相应地变为
110
r
n n rj n j j y y h f β+-+==+∑. （7.3.8）
由于式（7.3.8）右端第二项含有111(,)n n n f f x y +++=(可能是非线性表达式)，所以式（7.3.8）属于隐式公式，称为Adams 隐式公式，且插值区间包含了积分区间
[]1n n x x +，因此属于“内插”，称为Adams 内插法，也称为Adams-Moulton 法。显
然当0k =时Adams 内插法就是改进的Euler 法。
当1,2,3,4k =时，表7.3.1和表7.3.2分别给出了Adams 显式公式和隐式公式的系数值。
表7.3.1 Adams 外插法系数值
j
0 1 2 3 4 5 0j α
1 12j α 3 -1 212j α 23 -16 5 324j α 55 -59 37 -9 4720j α 1901 -2774 2616 -1274 251 51440j α
4277
-7923
9982
-7298
2877
-475
表7.3.2 Adams 内插法系数值
j
0 1 2 3 4 5 0j β
1 12j β 1 1 212j β 5 8 -1 324j β 9 19 -5 1 4720j β 251 646 -264 106 -19 51440j β
475
1427
-798
482
-173
27
利用插值多项式的余项可以求出Adams 方法的局部截断误差，对指定的k ，表7.3.3列出了它们的局部截断误差主项。
表7.3.3 Adams 方法的局部截断误差主项
k
1
2 3
显式公式 ()3
'''
512i h y x ()4(4)
38
i h y
x ()5(5)
251720i h y
x 隐式公式
()3
(3)
112
i h y
x -
()4
(4)
124
i h y
x -
()5
(5)
19
720
i h y
x -
应该指出，用数值积分法只能构造一类特殊的多步法，其系数满足
1k α=，01α=-，0j α=当0,j k ≠。下面介绍更一般的待定系数法。
7.3.2 待定系数法
为了分析一般线性多步法的局部截断误差，令
[(),][()'()]k
j
j j L y x h y x jh h y x jh α
β==
+-+∑ （7.3.9）
设()y x 是初值问题的解，将()y x jh +和'()y x jh +在点x 用Taylor 公式展开，代入（7.3.9）按h 的同次幂合并同类项，得
2()
012[(),]()'()''()()p p p L y x h c y x c hy x c h y x c h y x =+++++
其中
001112011
1
12122()11(2)(2)
!(1)!k k k p p p p p k k c c k c k k
p p ααααααβββαααβββ--=+++??
=+++-+++??
??
?=+++-+++-??
，1,2,p =
若()y x 有2p +次连续微商，则可选取足够大的k 和j α，j β使010p c c c ==== ，而10p c +≠，即选j α，j β满足
01120111121202()011(2)(2)0
!
(1)!k k k p p p p k k k k k p p ααααααβββαααβββ--+++=??
+++-+++=??
???+++-+++=-??
（7.3.10） 此时有
1
(1)
2
1[(),]()()p p p p L y x h c h
y
x O h
++++=+
令'()(,())y x f x y x =，则
00
1
(1)
2
1[(),][()'()]
[()(,()]()()
k
n j
n j n j k j
n j n n j p p p p n L y x h y x jh h y x jh y x jh h f x jh y x jh c h
y
x O h
α
βα
β==++++=
+-+=
+-++=+∑∑
于是由满足线性方程组（7.3.10）的j α，j β得到的线性多步法（7.3.1）的局部截断误差为1(1)21,1()()()p p p p n k p n R c h y x O h O h +++++=+=，可以证明此线性多步法的整体
截断误差的阶是()p O h ，所以此线性多步法为p 阶k 步法。显然p 的大小和k 有关。
因为线性多步法（7.3.1）可以相差一个非零乘数，所以不妨设1k α=。当0k β=时n k y +可用12,,,n k n k n y y y +-+- 直接表示，称为显方法。反之，当0k β≠时，求n k y +需要解一个方程（一般用迭代法），称为隐方法。用待定系数法构造线性多步法的一个基本要求是选取j α，j β使局部截断误差的阶尽可能高。
（1）Milne 方法
作为待定系数法的一个应用，下面讨论一般的二步法。此时2k =，21α=，其余5个系数01012,,,,ααβββ由01230c c c c ====确定，即满足方程组：
0011101221123112102()0
1(4)(2)0
211
(8)(4)062c c c c αααβββαββαββ=++=??=+-++=???=+-+=??=+-+=??，解之得 1000102
0(1)
1(15)122(1)3
1(5)12
ααβαβαβα=-+??
?=-+?
??=-???=+??，其中0α为任意常数。 所以一般二步法为
201002010(1)[(5)8(1)(15)]12
n n n n n n h y y y f f f ααααα++++--+=
++--+ （7.3.11）
且由1
1
121211(2)(2
)!
(1)!
p
p
p p p k k c k k
p p αααβββ--=
+++-
+++- 得
41120111(16)(8)(1)246
24
c αββα=+-+=-+ 51120111(32)(16)(1713)
120
24
360
c αββα=
+-+=-
+
所以当01α≠-时40c ≠，此时的二步法（7.3.11）是3阶2步法；当01α=-时，
40c =，但50c ≠，此时方法（7.3.11）化为
221(4)3
n n n n n h y y f f f +++=+
++ （7.3.12）
这是4阶2步法，是具有最高阶的2步法，称为Milne 方法，它的余项是
()
()()55
61445i h y
x O h +。这一方法也可用Simpson 公式导出。此外，若取0
α=，
则二步法（7.3.11）为2步Adams 内插法；若取05α=-，则二步法（7.3.11）为
显式方法。
（2）Hmaming 方法
用待定系数法容易求出，4阶3步的公式
121113(9)(2).8
8
n n n n n n y y y h f f f +-+-=
-+
+-
（7.3.13）
这就是著名的4阶Hmaming 公式，它的余项为 ()
()()55
6140i h y
x O h -
+。
更多常用的线性多步法和多步法计算中的问题等可以阅读相关参考文献，在此不再赘述。
关于线性多步法的稳定性、收敛性和误差估计，以及绝对稳定性和相对稳定性等基本理论问题也请参阅相关的参考文献。
7.4 非线性单步法——Runge-Kutta 法
Euler 法是最简单的单步法，单步法不需要附加初值，所需的存储量小，改变步长灵活，但是线性单步法的阶最多是2。本节介绍非线性（关于f ）高阶单步法，主要介绍Runge-Kutta 法。
7.4.1 Taylor 展开法
设初值问题（7.2.1）的解充分光滑。将1()i y x +在i x 点用Taylor 公式展开：
'
2''
()
(1)
(1)
1111()()()()()()2!
!
(1)!
p p p p n n n n n n y x y x hy x h y x h y
x h
y
p p ξ+++=++
++
+
+
（7.4.1） 其中n ξ是介于n x 与1n x +之间的常数，()(),1,2,...,1k n y x k p =+，是未知函数()y x 在n x 点的k 阶导数，但它的值可以利用微分方程本身计算：
'''
'''
''''''''
'''2'''()
()(,())
()(,())(,())()()(,())2(,())()(,())(())(,())()()n n n n x n n y n n n n xx n n xy n n n yy n n n y n n n p n y x f x y x y x f x y x f x y x y x y x f x y x f x y x y x f x y x y x f x y x y x y
x ?=?=+??=+???++?
??=? （7.4.2） 舍去（7.4.1）中的含n ξ的项，则得到如下求解初值问题的非线性单步法：
'
2
''()
1112!
!
p p n n n
n
n
y y hy h y h y p +=++
++
（7.4.3）
其中()(),1,2,...,k n y x k p =按（7.4.2）计算。这一方法的局部截断误差为()1p O h +，
因此它是非线性p 阶单步法。由于需要计算(,)f x y 的高阶偏导数，计算工作量太大，所以一般不用（7.4.3）作数值计算，但可用它计算附加初值。
7.4.2 Runge-Kutta 法
单步法的一般形式为11(,,,)n n n n n y y h x y y h ?++=+，
1(,,,)n n n x y y h ?+称为增量函数。显然，Euler 法的1(,,,)(,)n n n n n x y y h f x y ?+=，Taylor 展开法的增量函数
'
''
1
()
'
1
(1)
11111(,,,)2!
!
2!
!
p p p p n n n n n n
n n n
x y y h y hy h
y f hf h
y p p ?---+=+
++
=+
++
。
Taylor 展开法是用f 在同一点(,)n n x y 的高阶导数来计算1(,,,)n n n x y y h ?+，计算量较大。Runge-Kutta 方法，简称R-K 方法，是用f 在区间1[,]n n x x +上的若干函数值来计算1(,,,)n n n x y y h ?+，使单步法局部截断误差的阶和Taylor 展开法的相等。将初值问题在区间1[,]n n x x +上写成积分形式：
11()()(,())n n
x n n x y x y x f x y x dx ++=+
?
（7.4.4）
在区间[,]n n x x h +上取m 个点123
n m
n t x t t t x
h =≤≤≤≤≤+ ，
若知道
(,()),1,2,i i
i K f t y t i m == ，则用它们的一次组合去近似（7.4.4）中的(,())f x y x ，
从而得到数值方法：
1m
n n i i i
y y h c K +=+∑ （7.4.5）
问题是如何计算i K （因()i y t 未知）和选取组合系数i c ，使方法（7.4.5）局部截断误差的阶足够高。一个直观的想法是：设1(,())n n K f x y x =已经知道，由Euler 法，2121111211()()()(,())()()y t y t t t f t y t y t t t K =+-=+-，于是
222212112211(,())(,()())(,())n K f t y t f t y t t t K f t y t t K ==+-=+-。
同样利用Euler 法又可以算出
333312113222211322(,())(,()()())(,()())
n K f t y t f t y t t t K t t K f t y t t K t t K ==+-+-=+-+-如此继续下去便解决了,1,2,,i K i m = 的计算问题。
为了便于推导，令1,1,2,,i i n i t t a h x a h i m =+=+= ，则Runge-Kutta 法的一般形式如下：
11122122211333113221122,11(...)(,)(,)
(,())......(,(...))
n n m m n n n n n n m n m n m m m m m y y h c K c K c K K f x y K f x a h y hb K K f x a h y h b K b K K f x a h y h b K b K b K +--=++++??
=??=++??
=+++???=+++++?? （7.4.6） 其中，01,1,2,,i a i m ≤≤= ，
12...1m c c c +++= （7.4.7）
212
3132312,1.........m m m m m b a b b a b b b a -=??
+=???
?+++=?
（7.4.8） 系数{}{},i ij a b 和{}i c 按下列原则确定：将,1,2,,i K i m = 关于h 展开，以之代入增量函数11122(,,,)(...)n n n m m x y y h c K c K c K ?+=+++中进行合并化简，令h 的k 次幂(1,2,...,)k h k p =的系数与Taylor 展开法增量函数
'
''
1
()
'
1
(1)
11111(,,,)2!
!
2!
!
p p p p n n n n n n
n n n
x y y h y hy h
y f hf h
y p p ?---+=+
++
=+
++
的
同次幂的系数相等，结合（7.4.7）和（7.4.8）求出系数{}{},i ij a b 和{}i c ，如此得到的算法（7.4.6）称为m 级p 阶Runge-Kutta 法。
下面推导一些常用的计算方案。将()n y x h +展开到h 的三次幂：
'
2''
3(3)
4
12
3
2
11()()()()()()()
2!
3!
()(,())((,())(,())(,()))
2
[(,())2(,())(,())(,())(,())6
(n n n n n n n n n x n n y n n n n xx n n xy n n n n yy n n n n y n y x y x h y x hy x h y x h y
x O h h
y x hf x y x f x y x f x y x f x y x h
f x y x f x y x f x y x f x y x f x y x f x +=+=++
+
++=+++++++4
23
,())((,())(,())(,()))]()
11()[()()]()(,,)
2
6
n x n n y n n n n n y n T n n y x f x y x f x y x f x y x O h y x h f hF h Ff G O h y x h x y h ?++=++
+
++=+其中，x y F f ff =+，22xx xy yy G f ff f f =++。 而且有展式
1(,)n n K f x y f
==
222112212223
212112
2
3
22(,)(,)
1()(2)()
2
1()
2
n n n n x y y xy yy K f x a h y b K f x a h y a K f ha f K f h a f K f K f O h f ha F h a G O h =++=++=++++++=++
+
2
2
3
33311322323231(,)()()2
n n y K f x a h y b K b K f ha F h a b f F a G O h =+++=+++
+
于是由Runge-Kutta 法得
11122332
2
2
3
123223323232233()(,,)
1[()()(2())()]
2
n n n n n n y y y h c K c K c K y h x y h y h c c c f h a c a c F h a b c f F a c a c G O h ?+=+++=+=++++++
+++比较(,,)T n n x y h ?和(,,)n n x y h ?的h 同次幂系数，可得以下具体计算方案。
（1）一级Runge-Kutta 法
令m =1，此时1231,0c c c ===，容易得出一级Runge-Kutta 法：
1n n y y hf
+=+
可见一级Runge-Kutta 法就是Euler 法。 （2）二级Runge-Kutta 法 令m =2，此时1231,0c c c +==，于是
2
2
3
22221(,,)()
2
n n x y h f ha c F h a c G O h ?=++
+
与(,,)T n n x y h ?比较0,h h 的系数，则得
121c c +=，2212
a c =
，
这有无穷多组解，从而有无穷多个二级二阶Runge-Kutta 算法。两个常见的方法分别是：
当取120,1c c ==时，得212
a =
，则
111(,)22n n n n n y y hf x h y hf +=++
+
称为中点法，这是一种修正的Euler 法。
当取1212
c c ==
时，得21a =，则
1111
((,))((,)(,))2
2n n n n n n n n n n n n y y h f f x h y hf y h f x y f x y hf ++=+
+++=+
++ 这是由1n n n y y hf +=+与改进Euler 法1
111
((,)(,))2
n n n n n n y y h f x y f x y +++=++构成的
预估-正法，也是一种改进的Euler 法。
（3）三级Runge-Kutta 法
令m =3，比较(,,)n n x y h ?与(,,)T n n x y h ?，令02,,h h h 的系数相等，则得
1231c c c ++=，223312
a c a c +=
，22223313
a c a c +=
，232316
a b c =
四个方程不能完全确定六个系数，含有两个自由参数，因此有无穷多个三级三阶Runge-Kutta 算法。两个常见的方法有： 若取12313,0,4
4
c c c =
==
，2332122,,3
3
3
a a
b =
=
=
，则得算法为：
1
1312
132(3)4
(,)11(,)33
22(,)33n n n n n n n n h y y K K K f x y K f x h y hK K f x h y hK +?
=++??
=???=++??
?
=++??
（7.4.9） 称此为Heun 三阶方法。
若取123121,,636
c c c =
=
=
，23321,1,2
2
a a
b =
==，则得算法为：
1
1231213
12(4)6
(,)
11(,)
22
(,2)
n n n n n n n n h y y K K K K f x y K f x h y hK K f x h y hK hK +?=+++??
=??
?=++??=+-+? （7.4.10） 称次为Kutta 三阶方法。当f 与y 无关时，这就是Simpson 公式。
常微分方程边值问题的数值解法
第8章 常微分方程边值问题的数值解法 引 言 第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。 只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微分方程的边值问题(boundary-value problem). 为简明起见，我们以二阶边值问题为 则边值问题(8.1.1)有唯一解。 推论 若线性边值问题 ()()()()()(),, (),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤?? ==? (8.1.2) 满足 (1) (),()p x q x 和()f x 在[,]a b 上连续； (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。 求边值问题的近似解，有三类基本方法： (1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解； (2) 有限元法(finite element method)；
(3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。 差分法 8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法 设二阶线性常微分方程的边值问题为 (8.2.1)(8.2.2) ()()()(),,(),(), y x q x y x f x a x b y a y b αβ''-=<
一阶常微分方程解法总结
第 一 章 一阶微分方程的解法的小结 ⑴、可分离变量的方程： ①、形如 )()(y g x f dx dy = 当0)(≠y g 时，得到 dx x f y g dy )() (=，两边积分即可得到结果； 当0)(0=ηg 时，则0)(η=x y 也是方程的解。 例1.1、 xy dx dy = 解：当0≠y 时，有 xdx y dy =，两边积分得到)(2ln 2为常数C C x y += 所以)(112 12 C x e C C e C y ±==为非零常数且 0=y 显然是原方程的解； 综上所述，原方程的解为)(12 12 为常数C e C y x = ②、形如0)()()()(=+dy y Q x P dx y N x M 当0)()(≠y N x P 时，可有 dy y N y Q dx x P x M ) () ()()(=，两边积分可得结果； 当0)(0=y N 时，0y y =为原方程的解，当0(0=） x P 时，0x x =为原方程的解。 例1.2、0)1()1(2 2 =-+-dy x y dx y x 解：当0)1)(1(2 2 ≠--y x 时，有 dx x x dy y y 1 122-=-两边积分得到 )0(ln 1ln 1ln 22≠=-+-C C y x ，所以有)0()1)(1(22≠=--C C y x ； 当0)1)(1(2 2 =--y x 时，也是原方程的解； 综上所述，原方程的解为)()1)(1(2 2 为常数C C y x =--。 ⑵可化为变量可分离方程的方程： ①、形如 )(x y g dx dy = 解法：令x y u =，则udx xdu dy +=，代入得到)(u g u dx du x =+为变量可分离方程，得到
常微分方程解题方法总结.doc
常微分方程解题方法总结 来源：文都教育 复习过半， 课本上的知识点相信大部分考生已经学习过一遍 . 接下来， 如何将零散的知 识点有机地结合起来， 而不容易遗忘是大多数考生面临的问题 . 为了加强记忆， 使知识自成 体系，建议将知识点进行分类系统总结 . 著名数学家华罗庚的读书方法值得借鉴， 他强调读 书要“由薄到厚、由厚到薄”，对同学们的复习尤为重要 . 以常微分方程为例， 本部分内容涉及可分离变量、 一阶齐次、 一阶非齐次、 全微分方程、 高阶线性微分方程等内容， 在看完这部分内容会发现要掌握的解题方法太多， 遇到具体的题 目不知该如何下手， 这种情况往往是因为没有很好地总结和归纳解题方法 . 下面以表格的形 式将常微分方程中的解题方法加以总结，一目了然，便于记忆和查询 . 常微分方程 通解公式或解法 ( 名称、形式 ) 当 g( y) 0 时，得到 dy f (x)dx ， g( y) 可分离变量的方程 dy f ( x) g( y) 两边积分即可得到结果； dx 当 g( 0 ) 0 时，则 y( x) 0 也是方程的 解 . 解法：令 u y xdu udx ，代入 ，则 dy 齐次微分方程 dy g( y ) x dx x u g (u) 化为可分离变量方程 得到 x du dx 一 阶 线 性 微 分 方 程 P ( x)dx P ( x) dx dy Q(x) y ( e Q( x)dx C )e P( x) y dx
伯努利方程 解法：令 u y1 n，有 du (1 n) y n dy ， dy P( x) y Q( x) y n（n≠0,1）代入得到du (1 n) P(x)u (1 n)Q(x) dx dx 求解特征方程：2 pq 三种情况： 二阶常系数齐次线性微分方程 y p x y q x y0 二阶常系数非齐次线性微分方程 y p x y q x y f ( x) (1)两个不等实根：1, 2 通解： y c1 e 1x c2 e 2x (2) 两个相等实根：1 2 通解： y c1 c2 x e x (3) 一对共轭复根：i , 通解： y e x c1 cos x c2 sin x 通解为 y p x y q x y 0 的通解与 y p x y q x y f ( x) 的特解之和. 常见的 f (x) 有两种情况： x （ 1）f ( x)e P m ( x) 若不是特征方程的根，令特解 y Q m ( x)e x；若是特征方程的单根，令特 解 y xQ m ( x)e x；若是特征方程的重根， 令特解 y*x2Q m (x)e x； （2）f (x) e x[ P m ( x) cos x p n ( x)sin x]
常微分方程数值解法的误差分析教材
淮北师范大学 2013届学士学位论文 常微分方程数值解法的误差分析 学院、专业数学科学学院数学与应用数学 研究方向计算数学 学生姓名李娜 学号 20091101070 指导教师姓名陈昊 指导教师职称讲师 年月日
常微分方程数值解法的误差分析 李娜 （淮北师范大学数学科学学院，淮北，235000） 摘要 自然界与工程技术中的很多现象，往往归结为常微分方程定解问题。许多偏微分方程问题也可以化为常微分方程问题来近似求解。因此，研究常微分方程的数值解法是有实际应用意义的。数值解法是一种离散化的数学方法，可以求出函数的精确解在自变量一系列离散点处的近似值。随着计算机计算能力的增强以及数值计算方法的发展，常微分方程的数值求解方法越来越多，比较成熟的有Euler 法、后退Euler法、梯形方法、Runge—Kutta方法、投影法和多步法，等等．本文将对这些解的误差进行分析,以求能够得到求解常微分数值解的精度更好的方法。 关键词：常微分方程, 数值解法, 单步法, 线性多步法, 局部截断误差
Error Analysis of Numerical Method for Solving the Ordinary Differential Equation Li Na (School of Mathematical Science, Huaibei Normal University, Huaibei, 235000) Abstract In nature and engineering have many phenomena , definite solution of the problem often boils down to ordinary differential equations. So study the numerical solution of ordinary differential equations is practical significance. The numerical method is a discrete mathematical methods, and exact solution of the function can be obtained in the approximation of a series of discrete points of the argument.With the enhanced computing power and the development of numerical methods,ordinary differential equations have more and more numerical solution,there are some mature methods. Such as Euler method, backward Euler method, trapezoidal method, Runge-Kutta method, projection method and multi-step method and so on.Therefore, numerical solution of differential equation is of great practical significance. Through this paper, error of these solutions will be analyzed in order to get a the accuracy better way to solve the numerical solution of ordinary differential. Keywords:Ordinary differential equations, numerical solution methods, s ingle ste p methods, l inear multi-step methods, local truncation error
常微分方程数值解法
i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。
常微分方程初值问题数值解法
常微分方程初值问题数值解法 朱欲辉 （浙江海洋学院数理信息学院, 浙江舟山316004） [摘要]：在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法.然而在生产实 际和科学研究中所遇到的问题往往很复杂, 在很多情况下都不可能给出解的解析表达式. 本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如Euler法、改进的Euler法、Runge－Kutta法以及线性多步法中的Adams显隐式公式和预测校正 公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例子, 分别用不同的方法计算出近似解, 从得出的结果对比各种方法的优缺点. [关键词]：常微分方程;初值问题; 数值方法; 收敛性; 稳定性; 误差估计 Numerical Method for Initial-Value Problems Zhu Yuhui (School of Mathematics, Physics, and Information Science, Zhejiang Ocean University, Zhoushan, Zhejiang 316004) [Abstract]:In the course about ordinary differential equations, the methods for analytic solutions of some typical equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems can’t be e xpressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods include Euler method, improved Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of theoretical analysis. [Keywords]:Ordinary differential equation; Initial-value problem; Numerical method; Convergence; Stability;Error estimate
各类微分方程的解法大全
各类微分方程的解法 1.可分离变量的微分方程解法 一般形式:g(y)dy=f(x)dx 直接解得∫g(y)dy=∫f(x)dx 设g(y)及f(x)的原函数依次为G(y)及F(x),则G(y)=F(x)+C为微分方程的隐式通解 2.齐次方程解法 一般形式:dy/dx=φ(y/x) 令u=y/x则y=xu,dy/dx=u+xdu/dx,所以u+xdu/dx=φ(u),即du/［φ(u)-u］=dx/x 两端积分,得∫du/［φ(u)-u］=∫dx/x 最后用y/x代替u,便得所给齐次方程的通解 3.一阶线性微分方程解法 一般形式:dy/dx+P(x)y=Q(x) 先令Q(x)=0则dy/dx+P(x)y=0解得y=Ce－ ∫P(x)dx,再令y=u e－∫P(x)dx代入原方程解得u=∫Q(x) e∫P(x)dx dx+C,所以y=e－∫P(x)dx［∫Q(x)e∫P(x)dx dx+C］ 即y=Ce－∫P(x)dx +e－ ∫P(x)dx∫Q(x)e∫P(x)dx dx为一阶线性微分方程的通解 4.可降阶的高阶微分方程解法 ①y(n)=f(x)型的微分方程 y(n)=f(x) y(n-1)= ∫f(x)dx+C1 y(n-2)= ∫［∫f(x)dx+C1］dx+C2 依次类推,接连积分n次,便得方程y(n)=f(x)的含有n个任意常数的通解②y”=f(x,y’) 型的微分方程 令y’=p则y”=p’,所以p’=f(x,p),再求解得p=φ(x,C1) 即dy/dx=φ(x,C1),所以y=∫φ(x,C1)dx+C2 ③y”=f(y,y’) 型的微分方程
令y ’=p 则y ”=pdp/dy,所以pdp/dy=f(y,p),再求解得p=φ(y,C 1) 即dy/dx=φ(y,C 1),即dy/φ(y,C 1)=dx,所以∫dy/φ(y,C 1)=x+C 2 5.二阶常系数齐次线性微分方程解法 一般形式:y ”+py ’+qy=0,特征方程r 2+pr+q=0 6.二阶常系数非齐次线性微分方程解法 一般形式: y ”+py ’+qy=f(x) 先求y ”+py ’+qy=0的通解y 0(x),再求y ”+py ’+qy=f(x)的一个特解y*(x) 则y(x)=y 0(x)+y*(x)即为微分方程y ”+py ’+qy=f(x)的通解 求y ”+py ’+qy=f(x)特解的方法: ① f(x)=P m (x)e λx 型 令y*=x k Q m (x)e λx ［k 按λ不是特征方程的根,是特征方程的单根或特征方程的重根依次取0,1或2］再代入原方程,确定Q m (x)的m+1个系数 ② f(x)=e λx ［P ｌ(x)cos ωx+P n (x)sin ωx ］型 令y*=x k e λx ［Q m (x)cos ωx+R m (x)sin ωx ］［m=max ﹛ｌ,n ﹜,k 按λ+i ω不是特征方程的根或是特征方程的单根依次取0或1］再代入原方程,分别确定Q m (x)和R m (x)的m+1个系数
常微分方程数值解法
第八章 常微分方程数值解法 考核知识点： 欧拉法，改进欧拉法，龙格-库塔法，单步法的收敛性与稳定性。 考核要求： 1． 解欧拉法，改进欧拉法的基本思想；熟练掌握用欧拉法，改进欧拉法、求微 分方程近似解的方法。 2． 了解龙格-库塔法的基本思想；掌握用龙格-库塔法求微分方程近似解的方 法。 3． 了解单步法的收敛性、稳定性与绝对稳定性。 例1 用欧拉法，预估——校正法求一阶微分方程初值问题 ? ??=-='1)0(y y x y ，在0=x （0，1）0.2近似解 解 （1）用1.0=h 欧拉法计算公式 n n n n n n x y y x y y 1.09.0)(1.01+=-+=+，1.0=n 计算得 9.01=y 82.01.01.09.09.02=?+?=y （2）用预估——校正法计算公式 1,0)(05.01.09.0)0(111)0(1=???-+-+=+=++++n y x y x y y x y y n n n n n n n n n 计算得 91.01=y ，83805.02=y 例2 已知一阶初值问题 ???=-='1 )0(5y y y 求使欧拉法绝对稳定的步长n 值。 解 由欧拉法公式 n n n n y h y h y y )51(51-=-=+ n n y h y ~)51(~1-=+
相减得01)51()51(e h e h e n n n -==-=-Λ 当 151≤-h 时，4.00≤ 第八章 常微分方程的数值解法 一．内容要点 考虑一阶常微分方程初值问题:?????==0 0)() ,(y x y y x f dx dy 微分方程的数值解：设微分方程的解y (x )的存在区间是[a,b ]，在[a,b ]内取一系列节 点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ；(一般采用等距节点，h=(b-a)/n 称为步长)。在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。 用数值方法，求得f(x k )的近似值y k ，再用插值或拟合方法就求得y(x)的近似函数。 (一)常微分方程处置问题解得存在唯一性定理 对于常微分方程初值问题:?????==0 0)() ,(y x y y x f dx dy 如果： (1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。 (2) ),(y x f 对于y 满足利普希茨条件，即 2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程?????==0 0)() ,(y x y y x f dx dy 的解存在且唯一，这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。 定义：任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+，其中),,(h y x k k Φ称为算法的增量函数。 收敛性定理：若一步方法满足： (1)是p 解的. (2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件. (3) 初始值y 0是精确的。则),()()(p h O x y kh y =-kh =x -x 0,也就是有 0x y y lim k x x kh 0h 0 =--=→)( (一)、主要算法 1．局部截断误差 局部截断误差：当y(x k )是精确解时，由y(x k )按照数值方法计算出来的1~ +k y 的误差y (x k+1)- 1~ +k y 称为局部截断误差。 注意：y k+1和1~ +k y 的区别。因而局部截断误差与误差e k +1=y (x k +1) -y k +1不同。 如果局部截断误差是O (h p+1),我们就说该数值方法具有p 阶精度。 常微分方程数值解法 【作用】微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题，大体上可以按以下几步： 1. 根据实际要求确定要研究的量(自变量、未知函数、必要的参数等)并确定坐标系。 2. 找出这些量所满足的基本规律(物理的、几何的、化学的或生物学的等等)。 3. 运用这些规律列出方程和定解条件。基本模型 1. 发射卫星为什么用三级火箭 2. 人口模型 3. 战争模型 4. 放射性废料的处理通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来” 的于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。 1. 改进Euler 法： 2. 龙格—库塔( Runge—Kutta )方法： 【源程序】 1. 改进Euler 法： function [x,y]=eulerpro(fun,x0,x1,y0,n);%fun 为函数，(xO, x1)为x 区间，yO 为初始值，n 为子 区间个数 if nargin<5,n=5O;end h=(x1-xO)/n; x(1)=xO;y(1)=yO; for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(fun,x(i),y(i)); y2=y(i)+h*feval(fun,x(i+1),y1); y(i+1)=(y1+y2)/2; end 调用command 窗口 f=i nlin e('-2*y+2*x A2+2*x') [x,y]=eulerpro(f,O,,1,1O) 2 x +2x , (0 < x < , y(0) = 1 求解函数y'=-2y+2 2. 龙格—库塔( Runge—Kutta )方法： [t,y]=solver('F',tspan ，y0) 这里solver为ode45, ode23, ode113，输入参数F是用M文件定义的微分方程y'= f (x, y)右端的函数。tspan=[t0,tfinal]是求解区间，y0是初值。 注：ode45和ode23变步长的，采用Runge-Kutta算法。 ode45表示采用四阶-五阶Runge-Kutta算法，它用4阶方法提供候选解，5阶方法控制误差，是一种自适应步长(变步长)的常微分方程数值解法，其整体截断误差为(△ 口人5解 决的是Nonstiff(非刚性)常微分方程。 实验4 常微分方程的数值解 【实验目的】 1．掌握用MATLAB软件求微分方程初值问题数值解的方法； 2．通过实例用微分方程模型解决简化的实际问题； 3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。 【实验内容】 题3 小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。 模型及其求解 火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=9.8m/s2。 在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。因此有如下二式： a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8 dh/dt=v 又知初始时刻t=0，v=0，h=0。记x(1)=h，x(2)=v，根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下： function [ dx ] = rocket( t,x ) a=[(32000-0.4*x(2)^2)/(1400-18*t)]-9.8; dx=[x(2);a]; end ts=0:1:60; x0=[0,0]; [t,x]=ode45(@rocket,ts,x0); h=x(:,1); v=x(:,2); a=[(32000-0.4*(v.^2))./(1400-18*t)]-9.8; [t,h,v,a]; 数据如下： t h v a 0 0 0 13.06 1.00 6.57 13.19 13.30 2.00 26.44 26.58 1 3.45 3.00 59.76 40.06 13.50 4.00 106.57 53.54 13.43 5.00 16 6.79 66.89 13.26 6.00 240.27 80.02 12.99 7.00 326.72 92.83 12.61 8.00 425.79 105.22 12.15 9.00 536.99 117.11 11.62 10.00 659.80 128.43 11.02 11.00 793.63 139.14 10.38 12.00 937.85 149.18 9.71 13.00 1091.79 158.55 9.02 14.00 1254.71 167.23 8.33 15.00 1425.93 175.22 7.65 16.00 1604.83 182.55 6.99 17.00 1790.78 189.22 6.36 18.00 1983.13 195.27 5.76 19.00 2181.24 200.75 5.21 20.00 2384.47 205.70 4.69 21.00 2592.36 210.18 4.22 22.00 2804.52 214.19 3.79 23.00 3020.56 217.79 3.41 24.00 3240.08 221.01 3.07 25.00 3462.65 223.92 2.77 26.00 3687.88 226.56 2.50 27.00 3915.58 228.97 2.27 -293- 第十五章 常微分方程的解法 建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如 22x y dx dy +=。于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。 §1 常微分方程的离散化 下面主要讨论一阶常微分方程的初值问题，其一般形式是 ?????=≤≤=0 )() ,(y a y b x a y x f dx dy （1） 在下面的讨论中，我们总假定函数),(y x f 连续，且关于y 满足李普希兹(Lipschitz)条 件，即存在常数L ，使得 |||),(),(|y y L y x f y x f ?≤? 这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。 所谓数值解法，就是求问题（1）的解)(x y 在若干点 b x x x x a N =<<<<=L 210 处的近似值),,2,1(N n y n L =的方法，),,2,1(N n y n L =称为问题（1）的数值解， n n n x x h ?=+1称为由n x 到1+n x 的步长。今后如无特别说明，我们总取步长为常量h 。 建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： （i ）用差商近似导数 若用向前差商h x y x y n n ) ()(1?+代替)('n x y 代入（1）中的微分方程，则得 )1,,1,0())(,() ()(1?=≈?+N n x y x f h x y x y n n n n L 化简得 ))(,()()(1n n n n x y x hf x y x y +≈+ 如果用)(n x y 的近似值n y 代入上式右端，所得结果作为)(1+n x y 的近似值，记为1+n y ， 则有 )1,,1,0() ,(1?=+=+N n y x hf y y n n n n L （2） 这样，问题（1）的近似解可通过求解下述问题 ?? ?=?=+=+) () 1,,1,0(),(01a y y N n y x hf y y n n n n L （3） 得到，按式（3）由初值0y 可逐次算出N y y y ,,,21L 。式（3）是个离散化的问题，称为差分方程初值问题。 i.常微分方程初值问题数值解法 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法--差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<<<<=L （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-L 方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t L 上的差分解1,,N u u L 。 二阶微分方程： 时为非齐次 时为齐次，0)(0)()()()(22≠≡=++x f x f x f y x Q dx dy x P dx y d 二阶常系数齐次线性微分方程及其解法： 2 122,)(2,,(*)0)(1,0(*)r r y y y r r q pr r q p qy y p y 式的两个根、求出的系数； 式中的系数及常数项恰好是，，其中、写出特征方程：求解步骤： 为常数； ，其中?'''=++?=+'+''式的通解：出的不同情况，按下表写、根据(*),321r r 二阶常系数非齐次线性微分方程 型 为常数； 型，为常数 ，]sin )(cos )([)()()(,)(x x P x x P e x f x P e x f q p x f qy y p y n l x m x ωωλλλ+===+'+'' 二阶常系数非齐次线性微分方程的一般形式是 ''+'+=y py qy f x () （1） 其中p q ,是常数。 方程（1）的通解为对应的齐次方程 0=+'+''qy y p y （2） 的通解Y 和方程（1）的一个特解*y 之和。即 *y Y y +=.我们已解决了求二阶常系数齐 次线性方程通解的问题，所以，我们只需讨论求二阶常系数非齐次线性微分方程的特解* y 的方法。 下面我们只介绍当方程（1）中的)(x f 为如下两种常见形式时求其特解*y 的方法。 一、 f x e P x x m ()()=?λ型 由于方程（1）右端函数f x ()是指数函数e x λ?与m 次多项式P x m ()的乘积，而指数 函数与多项式的乘积的导数仍是这类函数，因此，我们推测： 方程（1）的特解应为 y e Q x x *?=λ()( Q x ()是某个次数待定的多项式 ) y e Q x e Q x x x *??'=+'λλλ()() y e Q x Q x Q x x *?"=?+'+''λλλ[()()()]22 代入方程（1），得 e Q x p Q x p q Q x e P x x x m λλλλλ???''++'+++≡?[()()()()()]()22 消去e x λ?，得 ''++'+++≡Q x p Q x p q Q x P x m ()()()()()()22λλλ （3） 讨论 01、如果λ不是特征方程 r pr q 20++=的根。 即 02≠++q p λλ 由于P x m ()是一个m 次的多项式，欲使（3）的两端恒等，那未Q x ()必为一个m 次多项式，设为 Q x b x b x b x b m m m m m ()=++++--0111Λ 将之代入（3），比较恒等式两端x 的同次幂的系数，就得到以b b b b m m 01 1,,,,Λ-为未知数的m +1个线性方程的联立方程组，解此方程组可得到这m +1个待定的系数，并得到特解 y e Q x x m *?=λ() 02、如果λ是特征方程 r pr q 20++=的单根。 即 λλ20++=p q ，但 20λ+≠p 欲使（3）式的两端恒等，那么'Q x ()必是一个m 次多项式。 因此，可令 Q x x Q x m ()()=? 并且用同样的方法来确定)(x Q 的系数b b b b m m 0 11,,,,Λ-。 03、如果λ是特征方程 r pr q 20++=的二重根。 即 λλ20++=p q ，且 20λ+=p 。 欲使（3）式的两端恒等，那么''Q x ()必是一个m 次多项式 因此， 可令 Q x x Q x m ()()=?2 并且用同样的方法来确定)(x Q 的系数b b b b m m 011,,,,Λ-。 目录 1 引言 (1) 2 二阶常系数常微分方程的几种解法 (1) 2.1 特征方程法 (1) 2.1.1 特征根是两个实根的情形 (2) 2.1.2 特征根有重根的情形 (2) 2.2 常数变异法 (4) 2.3 拉普拉斯变化法 (5) 3 常微分方程的简单应用 (6) 3.1 特征方程法 (7) 3.2 常数变异法 (9) 3.3 拉普拉斯变化法 (10) 4 总结及意义 (11) 参考文献 (12) 二阶常微分方程的解法及其应用 摘要：本文通过对特征方程法、常数变易法、拉普拉斯变换法这三种二阶常系数常微分方程解法进行介绍，特别是其中的特征方程法分为特征根是两个实根的情形和特征根有重根的情形这两种情况，分别使用特征值法、常数变异法以及拉普拉斯变换法来求动力学方程，现今对于二阶常微分方程解法的研究已经取得了不少成就，尤其在二阶常系数线性微分方程的求解问题方面卓有成效。应用常微分方程理论已经取得了很大的成就，但是，它的现有理论也还远远不能满足需要，还有待于进一步的发展，使这门学科的理论更加完善。 关键词：二阶常微分方程;特征分析法；常数变异法；拉普拉斯变换 METHODS FOR TWO ORDER ORDINARY DIFFERENTIAL EQUATION AND ITS APPLICATION Abstract:This paper introduces the solution of the characteristic equation method, the method of variation of parameters, the Laplasse transform method the three kind of two order ordinary differential equations with constant coefficients, especially the characteristic equation method which is characteristic of the root is the two of two real roots and characteristics of root root, branch and don't use eigenvalue method, method of variation of constants and Laplasse transform method to obtain the dynamic equation, the current studies on solution of ordinary differential equations of order two has made many achievements, especially in the aspect of solving the problem of two order linear differential equation with constant coefficients very fruitful. Application of the theory of ordinary differential equations has made great achievements, however, the existing theory it is still far from meeting the need, needs further development, to make the discipline theory more perfect. Keywords:second ord er ordinary differential equation; Characteristic analysis; constant variation method; Laplasse transform 1 引言 数学发展的历史告诉我们，300年来数学分析是数学的首要分支，而微分方程常微分方程数值解法
常微分方程数值解法
常微分方程的数值解
15第十五章 常微分方程的解法
(整理)常微分方程数值解法
(完整版)专题一(二阶常微分方程解法)
二阶常微分方程的解法及其应用