浅谈常微分方程的数值解法及其应用[文献综述]

文献综述

信息与计算科学

浅谈常微分方程的数值解法及其应用

一、前言部分

微分方程差不多是和微积分同时先后产生的，苏格兰数学家耐普尔创立对数的时候，就讨论过微分方程的近似解.牛顿在建立微积分的同时，对简单的微分方程用级数来求解.

后来瑞士数学家雅各布?贝努利、欧拉、法国数学家克雷洛、达朗贝尔、拉格朗日等人又不断地研究和丰富了微分方程的理论.

微分方程的理论逐步完善的时候，利用它就可以精确地表述事物变化所遵循的基本规律，只要列出相应的微分方程，有了解方程的方法.微分方程也就成了最有生命力的数学分支.总之，力学、天文学、几何学等领域的许多问题都导致微分方程.在当代，甚至许多社会科学的问题亦导致微分方程，如人口发展模型、交通流模型等.因而微分方程的研究是与人类社会密切相关的. [1]

“常微分方程”是理学院数学系所有专业学生的重要专业基础课之一，也是工科、经济等专业必学内容之一.其重要性在于它是各种精确自然科学、社会科学中表述基本定律和各种问题的根本工具之一，换句话说，只要根据实际背景，列出了相应的微分方程，并且能（数值地或定性地）求出这种方程的解，人们就可以预见到，在已知条件下这种或那种“运动”过程将怎样进行，或者为了实现人们所希望的某种“运动”应该怎样设计必要的装置和条件等等.例如，我们要设计人造卫星轨道，首先，根据力学原理，建立卫星运动的微分方程，列出初始条件，然后求出解，即卫星运行轨道.随着物理科学所研究的现象在广度和深度两方面的扩展，微分方程的应用范围更广泛. [2]从数学自身的角度看，微分方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展.从这个角度说，微分方程变成了数学的中心. [3]总之，微分方程从它诞生起即日益成为人类认识并进而改造自然、社会的有力工具，成为数学科学联系实际的主要途径之一.文章就常微分的数值解法以及应用展开简单的论述。

二、主体部分

2.1微分方程概念介绍

2.1.1 微分方程概况

由一元函数得到的方程.即：称含有自变量，未知函数及其导数的关系式

22(,,,,...,)0n n dy d y d y F x y dx dx dx

=. （1） 为常微分方程.其中出现的最高阶导数的阶数，叫做常微分方程的阶.例如 dy dx

=x ，dy y dx = ，是一阶常微分方程. 22sin 0d g dt p

θθ+=是二阶常微分方程.设)(x y ?=定义于 区间J 上，有直到n 阶的导数，将它代入（1），使（1）变成关于x 的恒等式，即

()()(,(),,...,)0,n n d x d x F x x x J dx dx

???=∈. 就称y ＝()x ?为（1）的一个定义于J 上的解，并称J 为该解的定义区间. [4]

如果一个微分方程中出现多元函数的偏导数，或者说如果未知函数和几个变量有关，而且方程中出现未知函数对几个变量的导数，那么这种微分方程就是偏微分方程.

2.2微分方程产生的历史背景

微分方程差不多是和微积分同时先后产生的，苏格兰数学家耐普尔创立对数的时候，就讨论过微分方程的近似解.牛顿在建立微积分的同时，对简单的微分方程用级数来求解。后来瑞士数学家雅各布?贝努利、欧拉、法国数学家克雷洛、达朗贝尔、拉格朗日等人又不断地研究和丰富了微分方程的理论。

微分方程的形成与发展是和力学、天文学、物理学，以及其他科学技术的发展密切相关的.数学的其他分支的新发展，如复变函数、李群、组合拓扑学等，都对微分方程的发展产生了深刻的影响，当前计算机的发展更是为常微分方程的应用及理论研究提供了非常有力的工具. [5]

牛顿研究天体力学和机械力学的时候，利用了微分方程这个工具，从理论上得到了行星运动规律.后来，法国天文学家勒维烈和英国天文学家亚当斯使用微分方程各自计算出那时尚未发现的海王星的位置.这些都使数学家更加深信微分方程在认识自然、改造自然方面的巨大力量.

微分方程的理论逐步完善的时候，利用它就可以精确地表述事物变化所遵循的基本规律，只要列出相应的微分方程，有了解方程的方法.微分方程也就成了最有生命力的数学分支.总之，力学、天文学、几何学等领域的许多问题都导致微分方程.在当代，甚至许多社会

科学的问题亦导致微分方程，如人口发展模型、交通流模型等.因而微分方程的研究是与人类社会密切相关的. [6]

2.3 微分方程发展现状及其基本功能

在数学学科内部的许多分支中，微分方程是常用的重要工具之一，微分方程进一步发展的需要，有推动着其它数学分支的发展；相反，微分方程每一步进展都离不开其他数学分支的支援.数学的其他分支的新发展，如复变函数、李群、组合拓扑学等，都对微分方程的发展产生了深刻的影响.当前计算机的发展更是为微分方程的应用及理论研究提供了非常有力的工具.时至今日，可以说微分方程在所有自然科学领域和众多社会科学领域都有着广泛的应用，如自动控制、各种电子学装置的设计、弹道的计算、飞机和导弹飞行的稳定性的研究、化学反应过程稳定性的研究等.只要能够列出相应的微分方程，有了解方程的方法，利用它就可以精确地表述事物变化所遵循的基本规律.从微积分理论形成以来，人们一直用微分方程来描述、解释或预见各种自然现象，不断的取得了显著的成效.

[7] 2.4常微分方程的数值求解方法

2.4.1 Euler 法

Euler 法是最简单的数值方法，[,]a b 为求解良态初值问题'(,)y f t y =，0()y a y =的区间。实际上，下面的过程不是要找到满足该初值问题的可微函数，而是要生成点集{(,)}k k t y ，并且将这些点作为近似解，即()k k y t y ≈。如何构造“近似满足微方程”的“点集”呢？首先为这些点选择横坐标，为方便起见，将区间[,]a b 划分为M 个等距子区间，并选择网络点

k t a kh =+， k=0，1，……，M 其中b a h M

-=

（1） 值h 称为步长。然后近似解

'(,)y f t y = 在0[,]M t t 上， 00()y t y = （2） 设()y t ，'()y t 和''()y t 连续，；；利用泰勒定理将()y t 在0t t =处展开，对每个值t ，存在一个0t 和t 之间的值1c ，使得 ''2

'

10000()()()()()()2y c t t y t y t y t t t -=+-+ （3） 将'00()(,())y t f t y t =和10h t t =-代人等式（3），得到1()y t 的表示：

''2

11000()()()(,())2

y c h y t y t hf t y t =++ （4） 如果步长 h 足够小，则可以忽略 2 次项（包含2h 的项），得到

1000(,)y y hf t y =+ （5） 这就是欧拉近似。

重复该过程，就能得到近似解曲线()y y t =的一个点序列。欧拉方法的一般步骤是 1k k t t h +=+ ， 1(,)k k k k y y hf t y +=+ 其中 k = 0，1，……，M-1[8]（6）

2.4.2 泰勒级数法

泰勒级数法有着广泛的应用，并且是比较求解初值问题的各种不同数值方法的标准，它可设计为任意指定的精度。下面首先将泰勒定理用新的公式表示，使之适合于求解微分方程。 定理9.5（泰勒定理）设1()N y t C

+∈ 0[,]t b ，且()y t 在不动点0[,]k t t t b =∈处有N 次

泰勒级数展开：

1()()(,())()N k k N k k y t h y t hT t y t O h ++=++ （1） 其中，

()11()(,())!j N

j k N k k j y t T t y t h j -==∑ （2） ()1()(,())j j y t f t y t -=表示函数f 关t 的（1j -）次全导数。求导公式可以递归地计算： ''''(3)''''22(4)''2'''''''''3

232()()()2()2()

33()33()(33)(2)

t y t y tt yt y yy tt yt yy y t y ttt ytt yyt ty y yy yyy ttt ytt yyt yyy y tt yt yy y t f

y t f f y f f f

y t f f y f y f y f f f f f f f f f y f f y f y f y f y f y y f y f f f f f f f f f f f f f ==+=+=+++=++++=++++++=++++++23()()()t y yt yy y t y f f f f f f f f f f +++++ （3）

并且一般有

()(1)()(,())N N y

t P f t y t -= （4） 其中P 为导数算子

()P f t y

??=+??

区间0[,]M t t 上的初值问题'

()(,)y t f t y =的近似数值解可由各子区间1[,]k k t t +上的公式

（1）来推导。N 次泰勒方法的一般步骤为 3

23211...2!3!!

N N k k d h d h d h y y d h N +=+++++ （5） 其中在各步0,1,......,1k M =-有()(),1,2,......,j j k d y t j N ==。

N 次泰勒方法的最终全局误差是1

()N O h +阶的，因此可选择所需大小的N ，使得误差足够小。如果N 是固定，则理论上可以推导出步长h ，使之满足任意的最终全局误差。然而在实际运算中，通常用h 和/2h 计算两个近似结果集，然后比较其结果[9]。

2.4.3 龙格—库塔方法

泰勒方法的优点是最终全局误差的阶为()N

O h ，并且可以通过选择较大的 N 来得到较小的误差。然而泰勒方法的缺点是，需要先确定 N ，并且要计算高阶导数，它们可能十分复杂。每个龙格一库塔（Runge-Kutta ）方法都由一个合适的泰勒方法推导而来，使得其最终全局误差为()N O h 。一种折中方法是每步进行若干次函数求值，从而省去高阶导数计算。这种方法可构造任意 N 阶精度的近似公式。最常用的是N= 4 的龙格一库塔方法，它适用于一般的应用，因为它非常精确、稳定，且易于编程。许多专家声称，没有必要使用更高阶的方法，因为提高的精度与增加的计算量相抵消。如果需要更高的精度，则应该使用更小的步长或某种自适应方法。

4 阶龙格一库塔方法（RK4）可模拟N=4 的泰勒方法的精度。它基于如下对1k y +，的计算：

111223344k k y y w k w k w k w k +=++++ （1）

其中1k ，2k ，3k 和4k 形如

1211132213243415263(,)

(,)

(,)

(,)k k k k k k k k k hf t y k hf t a h y b k k hf t a h y b k b k k hf t a h y b k b k b k ==++=+++=++++ （2）

通过与 N = 4 阶的泰勒级数方法的系数匹配，使得局部误差为5()O h ，龙格和库塔得出了如下方程组：

11

232

4563

123421324322221324311

2

1

3

b a b b a b b b a w w w w w a w a w a w a w a w a =+=++=+++=++=++= （3）

3332132433134152631234315262223134152641361

4

1

()61()8

1()12

124

w a w a w a w a b w a b a b w a a b w a a b a b w a b w a b a b w a b b ++=++=

++=++== 该方程组有11个方程和13个未知量，必须补充两个条件才可以求解。最有用的选择是

112

a =

，20b = （4） 其余变量的解为 23134561234111,1,,,0,0,12221111,,,6336

a a

b b b b b w w w w = ========== （5） 将式（4）和（5）中的值代入式（2）和式（1），得到标准的4N =阶龙格—库塔方法，其描述如下。自初始点00(,)t y 开始，利用

12341(22)6

k k h f f f f y y ++++=+

（6） 生成近似值序列，其中[10] 1213243(,)

(,)22(,)22

(,)

k k k k k k k k f f t y h h f f t y f h h f f t y f f f t h y hf ==++=++=++ （7） 2.4.4 预报—校正方法

欧拉方法、休恩方法、泰勒方法以及龙格一库塔方法都称为单步长方法，因为它们只利用前一个点的信息来计算下一个点，即计算11(,)t y 时只使用了初始点00(,)t y 。一般地，只有k y 用来1k y +。当计算出若干个点之后，就可以利用几个已计算出的点来计算下一个点。以亚当斯一巴什福斯4步法的推导为例，计算1k y +需要3k y -，2k y -，1k y -和k y 。该方法不

是自启动的，要生成点{(,):4}k k t y k ≥，必须先给出其4个初始点00(,)t y ，

11(,)t y ，22(,)t y ，33(,)t y （可用前面各节中的方法完成）。

多步法的一个优点是，可以确定它的局部截断误差（local truncation error ，简称 L.T.E.），并可以包含一个校正项，用于在每一步计算中改善解的精确度。该方法还可以确定步长是否小到能得到1k y +的精确值，同时又大到能够免除不必要的和费时的计算。使用预报子和校正子的组合在每一步只需要进行两次函数(,)f t y 求值

[11]。

2.4.4.1亚当斯一巴什福斯一莫尔顿方法

亚当斯一巴什福斯一莫尔顿方法（Adams —Bashforth —Moulton ）是由基本微积分定理推导出的多步法：

11()()(,())k k t k k t y t y t f t y t dt ++=+

? (1)

预报子使用基于点332211(,),(,),(,)k k k k k k t f t f t f ------和(,)k k t f 的(,())f t y t 的拉格朗日多项式逼近值，并在区间1[,]k k t t +上对式（1）积分，这个过程产生亚当斯一巴什福斯预报子：

1321(9375955)24

k k k k k k h p y f f f f +---=+-+-+ （2） 校正子的推导类似。这时可以实用刚刚计算出的值1k p +。基于点22(,)k k t f --，11(,)k k t f --(,)k k t f 和新的点11111(,)(,(,))k k k k k t f t f t p +++++=构造(,())f t y t 的一个新的拉格朗日多项式逼近，然后在区间1[,]k k t t +上对该多项式积分，即可得到亚当斯一莫尔顿校正子

[12]：

1211(5199)24

k k k k k k h y y f f f f +--+=+-++ （3） 2.4.4.2 米尔恩—辛普森方法

米尔恩—辛普森方法是预报子基于区间31[,]k k t t -+上的对(,())f t y t 的积分：

1

313()()(,())k k t k k t y t y t f t y t dt +-+-=+

? （4）

预报子使用(,())f t y t 基于332211(,),(,),(,)k k k k k k t f t f t f ------和(,)k k t f 的拉格朗日多项式逼近，在区间31[,]k k t t -+上对它积分，得到米尔恩预报子：

13214(22)3

k k k k k h p y f f f +---=+-+ (5) 校正子的推导类似。此时值1k p +已知，基于点11(,)k k t f --，(,)k k t f 和新点11111(,)(,(,))k k k k k t f t f t p +++++=构造(,())f t y t 的新的拉格朗日多项式，然后在区间

11[,]k k t t -+上对该多项式积分，结果为大家所熟悉的辛普森公式[13]：

1111(4)3k k k k k h y y f f f +--+=+

++ （6） 2.5 常微分方程的应用

2.5.1 微分方程在传染病模型中的应用

随着卫生设施的改善、医疗水平的提高以及人类文明的不断发展，诸如霍乱、天花、等曾经肆虐全球的传染性疾病已经得到有效的控制。但是一些新的、不断变异着的传染病毒却悄悄向人类袭来。20世纪80年代四分险恶的艾滋病毒开始肆虐全球，至今仍在蔓延；2003年春来历不明的非典病毒突袭人间，给人们的生命财产带来极大的危害。长期以来，建立传染病的数学模型来描述传染病的传播过程，分析受感染人数的变化规律，探索制止传染病蔓延的手段等，一直是各国有关专家和官员关注的问题。这里我们介绍的是按照一般传播机理建立的数学模型如下。 )()()()())()()(()(21313211t p t p t p t dp t p t p t p b dt

t dp αγ-+-++=， )()()()()(2212t p d d t p t p dt

t dp βα++-=， （1） )()()()(3323t p d t p dt

t dp γβ+-=， )3,2,1(0,,,:00033022011=≥====i p p p p p p p t i 其中

对上述常微分方程组进行分析求解，就可以了解不同时刻传染病的动力学特征(比如：传染病病情的发展趋势，即各类人的人口数量)(t p i 的分布情况)[14].

2.5.2常微分方程在捕鱼业的应用

设某水域现有鱼量x ，由于受资源限制所能容纳的最大鱼量m x ，高自然增长率r ，捕捞增长率k ，按人口的逻辑模型建立微分方程。

kx x x rx dt dx m

--=)1( 要保持鱼量平衡

0=dt dx ，设平衡点为0x ，解得m x r k r x -=0 设),(x f dt

dx =考虑)(x f 在0x 的泰勒展式 )(0))(()(000x x x x x f x f -+-'= r k x f -=')(0

当)(0x f '＞0时 )(x f 与0x x -同号 0x 为不稳定平衡点

当)(0x f '＜0时 )(x f 与0x x -异号 0x 为稳定平衡点

)(0x f '＜0即 r ＞k 设)1()(1m

x x rx x f -= kx x f =)(2 由于k ＜r 曲线)(1x f 与)(2x f 有交点，因)(1x f 在原点切线为rx y = 解得，易知当2

0m x x =

时，取得最大捕捞量 r k 21=， m x r rx x f 4

21)(02== 最大捕捞量为m x r 4 [15] 三、总结部分

微分方程是在解决实际问题的过程中产生的，微分方程的研究又促进实际问题的解决， 同时也促进其他学科的发展.微分方程在物理、 工程、 力学、天文学、生物学、 医学、 经济学等诸多领域都有重要作用.目前，微分方程的实际背景广、应用性强的特点已受到广泛关注.由于同一类型的微分方程往往可以用来描述许多性质上颇不相同的自然现象，对一些重要的微分方程开展研究，可以有多方面的应用前景，并可望在新兴学科或边缘学科的开发中及时地发挥作用. [16]

微分方程在几百年的发展史上， 在每个时期最新的技术科学中建立了自已的落脚点， 得到了实际的应用，从而刺激它飞跃地发展，在完善理论的同时继续深入其它新兴技术学科

领域近几十年来，世界科学技术进入核能、火箭、人造卫星、数字时代，微分方程定性理论及方法不论在应用上、理论上均不断地扩展着自身的领域，显示出前所未有的强大生命力.它的理论和方法，过去和现在都对力学、天文、物理、化学、生物、各种技术科学(如自动控制、无线电电子学等)及若干社会科学( 如人口理论、经济预测等)提供了有力的工具，后者反过来也不断地向它提出新的间题，刺激着它不断地向前发展. [17]

四、参考文献

[1] 张良勇,董晓芳. 常微分方程的起源与发展[J]. 高等函授学报(自然科学版),2006(3): 34-38.

[2] 黄赛.常微分方程发展的主要历史沿革[J].教育与职业.2006(2):36-41.

[3] 黄焕福.常微分方程课程建设初探[J].教育与职业,2007(14):146-148.

[4] 林建平.常微分方程早期发展概观[J]. 南京工程学院学报,2001(2):1-4.

[5] 周仲旺.几类特殊的常微分方程[J].潍坊学院学报.2003(6): 7.

[6] 邵晓锋,徐卫卫,李龙星.试析常微分方程模型的归结方法[J]. 黄冈职业技术学院报,

2007(1):37-39.

[7] 张同斌,管军军,焦万堂. 基本微分方程的应用研究[J].河南科学,2007(6):888-890.

[8] 胡建伟,汤怀民.微分方程数值解方法[M].科学出版社,2003(4).

[9] John H. Mathews,Kurtis D.Fink.Numerical Methods Using MATLAB Fourth Edition[M].电子工业大学出版社.2005(12).

[10] 王高雄,周之铭等.常微分方程[M].高等教育出版社,2008(1).

[11] John H. Mathews,Kurtis D.Fink.Numerical Methods Using MATLAB Fourth Edition[M].电子工业大学出版社.2005(12).

[12] J.W.Thomas.偏微分方程的数值解法[M].世界图书出版公司,1995(8).

[13] 郭晓梅.常微分方程的数值解法[J].网络财富，2009(8).

[14] 姜启源,谢金星,叶俊.数学模型[M].高等教育出版社.2003(8).

[15] James Stewart. Calculus[M].Fifth Edition. Beijing:高等教育出版社.2004(7): 1141- 1161.

[16] 欧阳瑞,孙要伟.常微分方程在数学建模中的应用[J]. 宿州教育学院学报,2008(2):

[17] James Stewart. Calculus[M].Fifth Edition. Beijing:高等教育出版社.2004(7):

1141- 1161.