常微分方程常用数值解法

第一章 绪论

1.1 引言

常微分方程是现代数学的一个重要分支，是人们解决各种实际问题的有效工具。微分方程的理论和方法从17世纪末开始发展起来，很快成了研究自然现象的强有力工具，在17到18世纪，在力学、天文、科学技术、物理中，就已借助微分方程取得了巨大的成就。1864年Leverrer根据这个方程预见了海王星的存在，并确定出海王星在天空中的位置。现在，常微分方程在许多方面获得了日新月异的应用。这些应用也为常微分方程的进一步发展提供了新的问题，促使人们对微分方程进行更深入的研究，以便适应科学技术飞速发展的需要。

研究常微分方程常用数值解是数学工作者的一项基本的且重要的工作。在国内外众多数学家的不懈努力，使此学科基本上形成了一套完美的体系。微分方程的首要问题是如何求一个给定方程的通解或特解。到目前为止，人们已经对许多微分方程得出了求解的一般方法。由于在生产实际和科学研究中所遇到的微分方程问题比较复杂，使这些问题的解即使能求出解析表达式，也往往因计算量太大而难于求出，而对于一些典型的微分方程则可以运用基本方法求出其解析解，并可以根据初值问题的条件把其中的任意常数确定下来。

由于求通解存在许多困难，人们就开始研究带某种定解条件的特解。首先是Cauchy对微分方程初始解的存在惟一性进行了研究。目前解的存在惟一性、延拓性、大范围的存在性以及解对初始解和参数的延续性和可微性等理论问题都已发展成熟。与此同时，人们开始采取各种近似方法来求微分方程的特解，例如求微分方程数值解的Euler折线法、Runge-Kutta法等，可以求得若干个点上微分方程的近似解。最后，由于当代高科技的发展为数学的广泛应用和深入研究提供了更好的手段。用计算机结合Matlab软件求方程的精确解、近似解，对解的性态进行图示和定性、稳定性研究都十分方便有效。

本章先介绍常微分的一般概念、导出微分方程的一些典型例子及求解微分方程的思路分析。从而得到常微分方程的常用数值解法。

1.2 常微分方程的概念

1.常微分方程的定义

含有未知量的等式称为方程，它表达了未知量所必须满足的某些条件。一般说来，凡含有自变量、未知函数以及未知函数的导数或微分的方程称为微分方程。如果微分方程中的未知函数只依赖于一个自变量，则称为常微分方程；如果未知函数依赖于两个或多个的自变量，并且在方程中出现偏导数，则称为偏微分方程。

在在一个微分方程中，未知函数最高阶导数的阶数，称为方程的阶数。如果一个微分方程关于未知函数及其各阶导数都是线性的，则称它为线性微分方程，否则称之为非线性微分方程。

本论文主要介绍常微分方程，也简称微分方程。

以y错误！未找到引用源。为未知函数，x错误！未找到引用源。为自变量的一阶常微分方程的一般形式可表示为微分方程，

F(x,y,y,)0, (1.2.1) 将（1.2.1）中y,错误！未找到引用源。解出，则得到方程

y,f(x,y), (1.2.2)

或

M(x,y)dxN(x,y)dy0, (1.2.3)

也称（1.2.1）为一阶隐式微分方程，（1.2.2）为一阶显式微分方程，（1.2.3）为一阶微分方程的微分形式。

错误！未找到引用源。阶隐式方程的一般形式为

F(x,y,y,,,y(n))0, (1.2.4)

错误！未找到引用源。阶显式方程的一般形式

y(n)f(x,y,y,,,y(n1)), (1.2.5)

方程（1.2.4）中，如果函数错误！未找到引用源。对未知函数错误！未找到引用源。和它的各阶导数y,y,,,y(n1)错误！未找到引用源。都是一次的，则称其为线性常微分方程，否则，称其为非线性微分方程。

以y错误！未找到引用源。为未知函数，错误！未找到引用源。为自变量的n错误！未找到引用源。阶线性微分方程具有如下形式：

y(n)p1(x)y(n1)pn1(x)y,pn(x)yf(x), (1.2.6)

2.常微分方程的解

设函数y(x)是定义在区间错误！未找到引用源。上的错误！未找到引用源。阶可微导数。如果把错误！未找到引用源。代入方程（1.2.4）后能使其成为恒等式，即

Fx,x,'x,..,nx0,xI

则称错误！未找到引用源。是微分方程（1.2.4）在区间错误！未找到引用源。上的一个解。

例如，yekx是微分方程dyky0在,的一个解.ytanx是微分方程dx

y'1y2在区间,的一个解. 22

如果关系式Fx,y0决定的隐函数y(x)是方程（1.2.4）的解，则我们称Fx,y0是（1.2.4）的隐式解。例如一阶微分方程

xdxydy0

有隐式解

x2y2c0.

我们把含有n个相互独立的任意常数c1,c2,...cn的解

y(x,c1,c2,..,cn)

称为n阶微分方程（1.2.4）的通解。

在通解之中当一组任意常数确定时，所得到确定的解称为特解。例如，yc1cosxc2sinx是二阶线性方程y''y0的通解，而ysinx,ycosxsinx都是其特解，其中c1,c2是任意常数。一般地，方程的特解可由其通解中任意常

数取确定的常数导出，且方程的通解不一定表示方程的所有解。

3.常微分方程初值问题

为了确定微分方程的一个特解，我们可以给出这个微分方程的所满足的定解条件，常见的定解条件是初始条件，即方程（1.2.4）在某一点x0所满足的条件： yx0y0,y'x0y'0,...,ynx0yn0. （1.2.7） 微分方程（1.2.4）连同初始条件（1.2.7）一起称为初始值问题。

一阶常微分方程初值问题是求解函数y(x)，且满足 dyf(x,y),axb, dx （1.2.8）

y(a)y.

其中f(x,y)错误！未找到引用源。为已知函数，y为给定的初值。 定理1.2.1 假设f(x,y)错误！未找到引用源。在矩形区域

(x,y)x[a,b],y(,)内连续，且关于变元yLipschitz连续，即存在正常数L，使得对任意x[a,b]及，成立不等式 f(x,y1)f(x,y2)Ly1y2, (1.2.9) 其中常数错误！未找到引用源。称为Lipschitz常数,则处置问题存在惟一解yxC[a,b]，yx连续依赖于初值y.

由常微分方程的基本理论，我们有：

换句话说，若是yx如下初值问题 dyf(x,y),axb, dx (1.2.10)

y(a)y.

的解错误！未找到引用源。，则

limyxyx,x[a,b].错误！未找到引用源。 0

因此问题（1.2.8）是适定的。

定理1.2.2 当假定函数f(x,y)错误！未找到引用源。满足定理1.2.1中的

条件时，也即微分方程（1.2.8）的解y(x)错误！未找到引用源。是适定的。

综上，高阶常微分方程初值问题一般为 dnydydn1yf(x,y,,,n1),axb,dxndxdx i (1.2.11) dy(a)y,i0,1,,n1.i1dxi

其中y1,,yn错误！未找到引用源。为给定值。

引进新的变量函数

dk1y(x)yk(x),axb,k1,2,,n. (1.2.12) k1dx

则初值问题（1.2.11）化成了一阶常微分方程组初值问题

dy1dxy2,

,

dyn1axb (1.2.13) yn,dx

dyndxf(x,y1,,yn),

yi(a)yi,i1,2,,n.

通过求解（1.2.13）得到（1.2.11）的解y(x).

1.3 常微分方程常用数值解法的思路

一般说来，对于某些典型微分方程可以通过初等积分法求解，而对于复杂微分方程，需要得到解在若干个点上的近似值或者便于计算的近似表达式即可。本论文讲研究常微分方程常用数值解法及其matlab程序设计。

对于一阶常微分方程的初值问题为

dyf(x,y), (1.3.1) dxy(x0)y0.

求微分方程初值问题（1.3.1）的数值解，就是求解函yy(x)在一系列离

散点错误！未找到引用源。xi(i1,2,,n)错误！未找到引用源。上精确值y(x1),y(x2),,y(xn)错误！未找到引用源。的近似值y1,y2,,yn错误！未找到引用源。。

由于计算的复杂及繁冗，在计算常微分方程初值问题的数值解是一般使用Matlab软件进行编程计算，按如下步骤：

（1）引入点列xi错误！未找到引用源。，其中xixi1hi,i1,2,,h，错误！未找到引用源。称为步长。为了便于使用计算机进行编程计算，一般取步长为定值，即hih错误！未找到引用源。，

xix0ih,i0,1,， (2.2)

（2）使用常微分方程常用数值解方法，即输入由yi1错误！未找到引用源。计算出yi(i1,2,,n)错误！未找到引用源。的递推公式。

（3）利用（2）中的公式逐步求出近似解y1,y2,,yn.错误！未找到引用源。

第二章 单步法

所谓单步法是指这类方法在计算yn1时，只用到前一步的值xn1,xn,yn，然后逐步往下计算。单步法的精度不高，这是它的缺点。但我们只要给定初值就可以计算，所以还是比较方便的。以下我们介绍两种典型的单步法：Euler法和Runge-Kutta法。

2.1 Euler法

为了计算出初始值问题 dyf(x,y), dx （2.1.1）

y(x0)y0.

在[x0,x0b]这一区间上若干个点上解的近似解，我们将此区间n等分，令hb,xkx0kh,k1,2,...,n,yk为（2.1.1）的解，即xxk当时精确值y(xk)的n

近似解。我们的任务是要用尽可能简单的方法求出与y(xk)尽可能接近的值yk。Euler法又称Euler折线法，Euler折线法是数值解法最简单的一种。

Euler折线法的基本思想是利用微分中值定理对（2.1.1）的解yyx函数进行近似。由于yx是可微函数，故

yx0hyx0y'h,

其中，是介于x0和x0h之间的一个值。当h比较小时，y'和y'x0相差不大，故用x0代替上式中的，就得到了的近似值

y1y0fx0,y0h.

当y1知道以后，再用类似方法来求yx在点xk的近似值yk.

ykyk1fxk1,yk1h,k1,2,..,n. （2.1.2） 这样就得到了（2.1.1）的解在点x0,x1,...,xn的近似解，从几何上看，Euler折线法就是在局部范围内用切线上的值去替代解曲线上的值。

Euler公式在实际计算时较少采用，但由于它的结构简单，易于分析，在理论上具有非常重要的意义。

2.1.1 向前Euler法

Euler折线法的计算量小，但误差大，后来人们对Euler折线法进行了改进： 由方程（2.1.1），在节点xn处成立

y'xnfxn,yxn. （2.1.3） 将yxn1在xn处利用Taylor展开式展开，得

h2

y''n, yxn1y(xn)hy'xn2

其中将nxn,xn1.将式（2.1.3）代入上式得

h2

y''n. （2.1.4） yxn1y(xn)hyxn,yxn2

在式（2.1.4）中略去高阶项

h2

y''n （2.1.5） Exn,h 2

并用yn和yn1分别代替式（2.1.4）中的y(xn)和yxn1，可得Euler公式（也称

向前Euler公式：

yn1ynhfxn1,yn,n0,1,..,N1. （2.1.6）

2.1.2 向后Euler法

向后Euler法和向前Euler法相差不大，只是将yxn1在xn处利用Taylor展开式展开，得

h2

y''n, yxn1y(xn)hy'xn2

其中将nxn1,xn.将式（2.1.3）代入上式得

h2

y''n. （2.1.7） yxn1y(xn)hyxn,yxn 2

在式（2.1.7）中略去高阶项，并用yn和yn1分别代替式（2.1.7）中的y(xn)和yxn1，可得向后Euler公式：

yn1ynhfxn1,yn1,n0,1,..,N1. （2.1.8） 其中y0y(a)ya，局部截断误差为

h2

Exn,hy''n1O(h2). 2

因此，向后Euler公式同向前Euler公式，为一阶方法。

向前Euler公式（2.1.6）和向后Euler公式（2.1.8）相比，向前Euler公式关于yn1是显式格式。而向后Euler公式中，yn1隐含在方程中，称这样的格式为隐式格式。

2.1.3 中点差分公式法

若将展开式（2.1.4）和（2.1.7）多展开一项，得

h2h3

1y''ny"'(n), yxn1y(xn)hyxn,yxn26

h2h3

y''ny"'(n2). yxn1y(xn)hyxn,yxn26

1其中n,n2xn1,xn1.将上述两式相减并整理得

h3

yxn1y(xn1)2hyxn,yxny"'(n), 3

其中nxn1,xn1.略去高阶项，并用yn1,yn和yn1代替上式中y(xn1),yxn 和y(xn1)，即得中心差分公式

yn1yn12hfxn,yn,n0,1,..,N1. （2.1.9） 其局部截断误差为

h3

y"'nO(h3). Exn,h3

即中心差分公式为二阶方法。

2.1.4 梯形公式

如果将（2.1.1）中微分方程两端从xn到xn1积分，则得

b y(xn1)y(xn)xn1

xnfx,ydx.

如果用梯形求积公式计算上式中的积分，则有

hh3

f"(n,y(n)). yxn1y(xn)[fxn,yxnfxn1,yxn1212

略去高阶项可得梯形公式：

yn1ynh[fxn,ynfxn1,yn1],n0,1,..,N1. (2.1.10) 2

其中y0y(a)ya，局部截断误差为

h3

Exn,hy"'nO(h3). 12

显然，梯形公式为二阶隐式公式。

2.1.5 改进Euler法

一般来说，隐式格式比显示格式具有较好的数值稳定性。但是，由于在隐式格式里，yn1往往满足的是一个非线性方程，因而需要使用迭代法得到的一个yn1近似值，然后代入隐式格式作校正，并以这个校正值作为yxn1的近似值。印版地，可以选择适当的显示格式计算预测值，我们把这样的格式称为预估一校正公式。预估一校正公式往往既便于计算又具有较好的数值稳定性。

如，对于梯形公式，可用向前Euler公式计算出的预测值，再用梯形公式修正得到：

y0n1ynhf(xn,yn),



 h0

yn1ynf(xn,yn)f(xn1,yn1).

2



上述公式为改进Euler公式。

其总体截断误差为O(h3).在步长h取得足够小的时侯，我们用改进Euler方法解得数值与理论值非常接近。

2.2 Runge-Kutta法

Runge-Kutta法是一种间接使用Taylor展开式来获取高阶方法的数值方法，也是一种特殊的单步法，具有相当实用价值。其基本做法是：首先，用函数fx,y在xn,yn附近的m个点上的函数值的线性组合来代替yx的导数。然后，通过Taylor展开式，确定其中的组合系数。这样既可以避免计算高阶导数，又可以获得较高的方法阶。

设m为一个正整数（代表使用fx,y的函数值的个数）， cii1,2,...,m和

ai,biji2,3,...,m,1ji为一些特定参数，则数值方法



yyh(cKcK..cK),

n1122mmn1

K1f(xn,yn),

K2f(xna2h,ynhb21K1), (2.2.1)

......

m1



Kmf(xnamh,ynhbmjKj).

j1称为m级Runge-Kutta公式。

2.2.1 二阶Runge-Kutta法

显然，公式（2.2.1）是显式的。如果对于该代数方程组的一组解，对应的方法阶为p，则得到一个m级p阶显式Runge-Kutta公式。

当m2时，公式（2.2.1）变为

yn1ynh(c1K1c2K2),



K1f(xn,yn),

Kf(xah,yhbK).

n2n2112

其局部截断误差为

Exn,hO(h3). 经计算得参数方程组：



c1c21,

1

c2a2, （2.2.2）

2

1cb.2212

1

1.若取c1c2,a21,b211,则对应两二级二阶方法（又称变形Euler公式）：

2

h

yy(K1K2),nn1

2



K1f(xn,yn),

Kf(xh,yhK).

nn1

2

132

2.若取c1,c2,a2b21,则对应两级二阶方法（又称两级Huen公式）：

443

h

yy(K13K2),nn1

4



K1f(xn,yn),

22

K2f(xnh,ynhK1).

33

2.2.2 三阶Runge-Kutta法

1.三级Kutta公式（三级三阶格式）

h

yy(K14K2K3),nn1

6



K1f(xn,yn),



hhKf(x,yK),nn1

222

Kf(xh,yhK2hK).

nn123

2.三级Huen公式（三级三阶格式）

hyy(K13K3),nn1

4



K1f(xn,yn),

 hh

Kf(x,yK),nn12

33



Kf(x2h,y2hK).3nn233

这两个常用方法在实际使用中能达到较低的精确要求。

2.2.3 古典Runge-Kutta法

古典Runge-Kutta公式是精度高、最实用的一种常微分方程常用数值解法。



yn1

ynh(K12K22K3K4),

6

K1f(xn,yn), 

Kxhh

2f(n2,yn2K1),

K3

f(xnh2

,ynhK2),K4fxnh,ynhK3.

2.2.4 一级隐式中点公式

设m为一个正整数（代表使用fx,y的函数值的个数），ai,biji2,3,...,m,1ji为一些特定参数，则数值方法

yn1ynh(c1K1c2K2..cmKm),m

K1

f(xna1h,ynhb1jKj

),

j1



mK2f(x

na2h,ynhb2jKj), (2.2.2)

j1......

m1Kmf(xnamh,ynhbmjKj).

j1称为m级隐式Runge-Kutta公式。 当m1,p2时，有一级隐式中点公式为

yn1ynhK1,



Khh

1f(xn2,yn2K1).2.2.5 二级隐式中点公式

当m1,p2时，有二级隐式中点公式为

cii1,2,...,m和

h

yy(K1K2),nn1

2



3132

K1f(xnh,ynhK1hK2),

6412

33231Kf(xh,yhKhK2).2nn1

6124

相对显式公式而言，对于同样的m，隐式格式不仅可获取较高的精度，而且具有较好的数值稳定性。

第三章 多步法

所谓多步法指：这类方法在计算yn1时，除了用到前一步的值xn1,xn,yn，之外，还要用到

xnp,ynp(p1,2,,k;k0)， 这前面k步的值，这个算法的代表就是阿达姆斯(Adams)方法。

本章介绍线性多步法。线性多步法多使用于求解步骤较多的情况。一般地，

k步线性多步法公式可写成

yn1jyxnjhjfxnj,ynj. （3.1）

j0

j1

k1k1

当10时，式（3.1）为显式格式。否则为隐式格式。 1.显式格式

首先，将常微分方程（2.1.1）两端在区间[xnp,xn1]上求定积分，得 yxn1yxnp

xn1xnp

ft,ytdt.

如果在上面等式中的定积分用关于结点xn,xn1,...,xnq的数值积分近似，可得如下

显式格式

yn1ynphjfxnj,ynj.

j0q

其中

1xn1

jljtdt,j0,1,...,q,

hxnp

而ljtj0,1,...,q为函数fx,yx关于结点xn,xn1,...,xnq的Lagrange插值函数的基函数。

当q3时，四步Adams显式格式

yn1yn

h

[55fxn,yn59fxn1,yn137fxn2,yn29fxn3,yn3]. 24

25155

hynnO(h6). 720

其局部截断误差为

En1

2.隐式格式

如果在等式（3.1）中的定积分用关于结点xn1,xn,...,xn1q的数值积分近似，可得如下隐式格式

yn1ynphjfxn1j,yn1j.

j0q

其中

j

1xn1

ljtdt,j0,1,...,q, xnph

而ljtj0,1,...,q为函数fx,yx关于结点xn1,xn,...,xn1q的Lagrange插值函数的基函数。

当q3时，三步Adams隐式格式 yn1yn

h

[9fxn1,yn119fxn,yn5fxn1,yn1fxn2,yn2]. 24

1955

hynnO(h6). 720

其局部截断误差为

En1

Adams显式与隐式法相比较而言：同阶数情况下，隐式格式的局部截断误差

系数的绝对值比显式格式的小，隐式格式的稳定性相对显式格式的要大。

总结

本论文总结了单步法（向前Euler法，向后Euler法，中心差分法，梯形法，改进Euler法，Runge-Kutta法）及线性多步法（Adams法）几种比较常用的常 微分方程数值解法，并通过数值算例对同一问题的解的误差计算比较，总结各种算法的优缺点。一般而言，对Euler法而言，求解常微分方程初值问题数值时计算简单、程序简洁并且计算量小，但是，计算精度却很低。中心差分法及梯形法较好地提高计算精确度，同时又保持计算简单，程序简洁的优点。改进的Euler法比Euler法得到的计算结果要好很多。我们用改进Euler方法解得数值与理论值非常接近。对Runge-Kutta法来说，越高阶的方法计算精度越高，但相应的计算量也越大。古典Runge-Kutta公式是精度高、最实用的一种常微分方程常用数值解法。而线性多步法多使用于求解步骤较多的情况。

我们要把常微分方程常用数值解法掌握并灵活运用并加以推广应用，来解决我们的实际生活中所遇到问题，为我国的发展服务，充分体现了常微分方程常用数值解法的重大意义。作为一名数学工作者，我们有责任也有义务继续探索常微分方程组的数值解法，并且在不断探索其解法的前提下，开发出更好更优的软件来解决实际算法问题。但是，其现有理论也还远远不能满足需要，有待于进一步

研究，使这门学科的理论更加完善。

数值实例

1.分别用改进Euler公式和古典Runge-Kutta公式解常微分方程初值问题

2x

,0x1,yy y

y(0)1

取步长h0.1，计算结果保存于下表1

中，精确解为y(x) 解：fun函数：

function f=fun(x,y); f=y-2*x/y;

精确解：

x=0:0.1:1; y=sqrt(1+2*x) y =

Columns 1 through 9

1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125

Columns 10 through 11 1.6733 1.7321

1.应用改进Euler法，运用matlab软件，程序为： format long h=0.1; x=0:h:1.0; y=zeros(size(x)); y(1)=1;

y1=sqrt(1+2*x); N=length(x); for n=1:N-1

K1=feval(@fun,x(n),y(n));

K2=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K1); y(n+1)=y(n)+h*K1; end y

2.应用古典Runge-Kutta法，运用matlab软件，程序为： format long h=0.1; x=0:h:1.0; y=zeros(size(x)); y(1)=1;

y1=sqrt(1+2*x); N=length(x); for n=1:N-1

K1=feval(@fun,x(n),y(n));

K2=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K1); K3=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K2); K4=feval(@fun,x(n)+h,y(n)+h*K3);

y(n+1)=y(n)+h/6*(K1+2*K2+2*K3+K4); end y

表1

2.分别用向前Euler

公式、向后Euler公式、中心差分公式和梯形公式解常微分方程初值问题

yy,0x1, y(0)1

取步长h0.1。

解：M文件：

x=0:0.1:1;

y=-exp(x);

1.应用向前Euler法，运用matlab软件，程序为：

h=0.1;

x=0:h:1;

y=-exp(x);

y1=zeros(size(x));

n=length(x);

y1(1)=1;

for i=1:n-1;

y1(i+1)=(1+(h/2))*y1(i);

end

error=abs(y-y1)

2.应用向后Euler法，运用matlab软件，程序为：

h=0.1;

x=0:h:1;

y=-exp(x);

y2=zeros(size(x));

n=length(x);

y2(1)=1;

for i=1:n-1;

y2(i+1)=(1/(1-(h/2)))*y2(i);

end

error=abs(y-y2)

3.应用中心差分法，运用matlab软件，程序为：

h=0.1;

x=0:h:1;

y=-exp(x);

y3=zeros(size(x));

n=length(x);

y3(1)=1;

y3(2)=-exp(h);

for i=1:n-1;

y3(i+2)=y3(i)+h*y3(i+1);

end

error=abs(y-y3)

4.应用梯形法，运用matlab软件，程序为：

h=0.1;

x=0:h:1;

y=-exp(x);

y4=zeros(size(x));

n=length(x);

y4(1)=1;

for i=1:n-1;

y4(i+2)=(8/(4-h)-1)*y4(i);

end

error=abs(y-y4).

3.分别用一级隐式中点公式、二级隐式中点公式解常微分方程初值问题

2xyy,0x1,yy(0)1 

取步长h0.1，计算结果保存于下表1

中，精确解为y(x)

解：fun函数：

function f=fun(x,y)

f=y-2*x/y;

1.应用一级隐式中点公式法，运用matlab软件，程序为：

format long

x=0:0.1:1.0;

y=sqrt(1+2*x);

h=0.1;

x=0:h:1.0;

y1=zeros(size(x));

y1=sqrt(1+2*x);

y1(1)=1;

N=length(x);

for n=1:N-1

error=1;

K1=feval(@fun,x(n),y1(n));

while(error>=1.0e-010)

K11=feval(@fun,x(n)+0.5*h,y1(n)+0.5*h*K1);

error=abs(K11-K1);

K1=K11;

end

y1(n+1)=y1(n)+h*K1;

end

e=abs(y-y1)

2.应用二级隐式中点公式法，运用matlab软件，程序为：

format long

x=0:0.1:1.0;

y=sqrt(1+2*x);

h=0.1;

x=0:h:1.0;

y2=zeros(size(x));

y2=sqrt(1+2*x);

y2(1)=1;

N=length(x);

for n=1:N-1

error=1;

K1=feval(@fun,x(n),y2(n));

K2=feval(@fun,x(n),y2(n));

while(error>=1.0e-010)

K11=feval(@fun,x(n)+(3-sqrt(3))/6*h,y2(n)+1/4*h*K1+(3-2*sqrt(3))/

12*h*K2);

K22=feval(@fun,x(n)+(3+sqrt(3))/6*h,y2(n)+1/4*h*K2+(3+2*sqrt(3))/

12*h*K1);

error=max(abs(K11-K1),abs(K22-K2));

K1=K11;

K2=K22;

end

y2(n+1)=y2(n)+0.5*h*(K1+K2);

end

e=abs(y-y2)

4.分别用Adans显式公式和Adans隐式公式解常微分方程初值问题

yyx22,0x1, y(0)1

取步长h0.1。

解：fun函数：

function f=fun(x,y);

f=-y+x.^2;

1.应用Adans显式法，运用matlab软件，程序为：

format long

h=0.1;

x=0:h:1.0;

y=zeros(size(x));

y1=x.^2-2*x+4-3*exp(-x);

y(1)=1;

N=length(x);

for n=1:N-1

K1=feval(@fun,x(n),y(n));

K2=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K1);

K3=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K2);

K4=feval(@fun,x(n)+h,y(n)+h*K3);

y(n+1)=y(n)+h/6*(K1+2*K2+2*K3+K4);

end

for n=5:N-1

M1=feval(@fun,x(n-1),y(n-1));

M2=feval(@fun,x(n-2),y(n-2));

M3=feval(@fun,x(n-3),y(n-3));

M4=feval(@fun,x(n-4),y(n-4));

y(n)=y(n-1)+h/24*(55*M1-59*M2+37*M3-9*M4);

end

y=abs(y-y1)

2.应用Adans隐式法，运用matlab软件，程序为：

format long

h=0.1;

x=0:h:1.0;

y=zeros(size(x));

y1=x.^2-2*x+4-3*exp(-x);

y(1)=1;

N=length(x);

for n=1:N-1

K1=feval(@fun,x(n),y(n));

K2=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K1);

K3=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K2);

K4=feval(@fun,x(n)+h,y(n)+h*K3);

y(n+1)=y(n)+h/6*(K1+2*K2+2*K3+K4);

end

for n=3:N-1

error=1

M2=feval(@fun,x(n),y(n));

M3=feval(@fun,x(n-1),y(n-1));

M4=feval(@fun,x(n-2),y(n-2));

temp1=y(n);

While(error>=1.0e-013)

M1=feval(@fun,x(n+1),temp1);

tempt2=y(n)+h/24*(9*M1+19*M2-5*M3+M4);

error=abs(temp1-temp2);

temp1=temp2;

end

y(n+1)=temp2;

end

y=abs(y-y1).

参考文献

[1] 戴嘉尊，邱建贤，微分方程数值解.南京：东南大学出版社，2002.

[2] 冯康. 数值计算方法.北京:国防工业出版社，1978.

[3] 周义仓，常微分方程及其应用.北京，科学出版社.2009.

[4] 李瑞遐，何志庆.微分方程数值解法.上海：华东理工大学出版社.2005.

[5] 徐树方，高平，张平文.数值线性代数.北京：北京大学出版社.2000.

[6] 余徳浩，汤华中.微分方程数值解法.北京：高等教育出版社.2003.

[7] 曾金平.微分计算方法.长沙：湖南大学出版社.2006.

[8] 谭浩强，C程序设计.北京：清华大学出版社.2009.

[9] 陆金甫，关治．偏微分方程数值解法.北京：清华大学出版社.1987.

[10] 王高雄，周之铭，朱思铭，王寿松.常微分方程.北京：高等教育出版社.

1987.

[11] 胡健伟，汤怀民.微分方程数值方法.北京:科学出版社.2000.

[12] 李信真，车刚明，欧阳洁，封建湖.计算方法.西安:西北工业大学出版社.

2000.

[13] 玛康.数值计算方法.杭州:浙江大学出版社.2003.

[14] 李荣华，玛果忱.微分方程数值方法(第三版).北京:高等教育出版社.1996.

[15] 刘会灯，朱飞.MATLAB编程基础与典型应用.北京:人民邮电出版社.2006.

[16] 苏金明，阮沈勇.MATLAB实用教程（第二版）.北京:电子工业出版社.2008.

[17] 刘卫国.MATLAB程序设计教程.北京:中国水利水电出版社.2005.

[18] 马知恩，周义仓.常微分方程定性与稳定性方法.北京:科学出版社.2007.

[19] C.W.吉尔.常微分方程初值问题的数值解法,费景高译.北京:科学出版社. 1978.

[20] 楼顺天，于卫，闰华梁.MATLAB程序设计语言西安:西安电了科技大学出版

社.1997.

[21] 李大侃，常微分方程数值解.浙江:浙江大学出版社.1994.

[22] 李立康，龄崇华，朱政华.微分方程数值解法.上海:复旦大学出版社.2003.

[23] 李庆扬，关治，白峰杉.数值计算原理.北京:清华大学出版社.2000.

[24] 李庆扬，王能超，易大义.数值分析.北京:清华大学出版社，施普林格

出版社.2001.

[25] 丁同仁，李承治.常微分方程.北京：高等教育出版社.1991.

[26] Hackbusch W. Iterative Solution of Large Sparse Systems of Equations.

New York:Springer-Verlag.1994

[27] Kincaid D,Cheney W. Numerical Analysis:Mathematics of Scientific

Computing.北京：机械工业出版社.2003.

[28] Kress R. Numerical Anaysis.New York:springer-Verlag.1998.

[29] Larson S,Thomee V. Partial Differential Equations with Numercial

Methods.Berlin:Springe-Verlag.2003.

[30] Quarteroni A,Valli A. Domain Decomposition Methods for Partial

Differential Equations.Oxford:Clarendon Press.1999.

[31] Lambert,J. D. Numerical Methods for Ordinary Differential Systems.

London: Wiley.1991.

[32] R. D.Richtmyer,K. W. Morton:Difference Methods of Initial value

Problems.1976.

[33] Butcher,J. C. The Numerical Analysis of Ordinary Differential

Equations.New York:John Wiley.1987.

[34] M. Braun. Differential Equations and Their Applications(Fourth edition).New York:Springer-Verlag.1993

[35] Sanz-Serna,J M and Ca1vo,M P. Numerical Hamiltonian Problems. London: Chapman & Hal.1994.

第一章 绪论

1.1 引言

常微分方程是现代数学的一个重要分支，是人们解决各种实际问题的有效工具。微分方程的理论和方法从17世纪末开始发展起来，很快成了研究自然现象的强有力工具，在17到18世纪，在力学、天文、科学技术、物理中，就已借助微分方程取得了巨大的成就。1864年Leverrer根据这个方程预见了海王星的存在，并确定出海王星在天空中的位置。现在，常微分方程在许多方面获得了日新月异的应用。这些应用也为常微分方程的进一步发展提供了新的问题，促使人们对微分方程进行更深入的研究，以便适应科学技术飞速发展的需要。

研究常微分方程常用数值解是数学工作者的一项基本的且重要的工作。在国内外众多数学家的不懈努力，使此学科基本上形成了一套完美的体系。微分方程的首要问题是如何求一个给定方程的通解或特解。到目前为止，人们已经对许多微分方程得出了求解的一般方法。由于在生产实际和科学研究中所遇到的微分方程问题比较复杂，使这些问题的解即使能求出解析表达式，也往往因计算量太大而难于求出，而对于一些典型的微分方程则可以运用基本方法求出其解析解，并可以根据初值问题的条件把其中的任意常数确定下来。

由于求通解存在许多困难，人们就开始研究带某种定解条件的特解。首先是Cauchy对微分方程初始解的存在惟一性进行了研究。目前解的存在惟一性、延拓性、大范围的存在性以及解对初始解和参数的延续性和可微性等理论问题都已发展成熟。与此同时，人们开始采取各种近似方法来求微分方程的特解，例如求微分方程数值解的Euler折线法、Runge-Kutta法等，可以求得若干个点上微分方程的近似解。最后，由于当代高科技的发展为数学的广泛应用和深入研究提供了更好的手段。用计算机结合Matlab软件求方程的精确解、近似解，对解的性态进行图示和定性、稳定性研究都十分方便有效。

本章先介绍常微分的一般概念、导出微分方程的一些典型例子及求解微分方程的思路分析。从而得到常微分方程的常用数值解法。

1.2 常微分方程的概念

1.常微分方程的定义

含有未知量的等式称为方程，它表达了未知量所必须满足的某些条件。一般说来，凡含有自变量、未知函数以及未知函数的导数或微分的方程称为微分方程。如果微分方程中的未知函数只依赖于一个自变量，则称为常微分方程；如果未知函数依赖于两个或多个的自变量，并且在方程中出现偏导数，则称为偏微分方程。

在在一个微分方程中，未知函数最高阶导数的阶数，称为方程的阶数。如果一个微分方程关于未知函数及其各阶导数都是线性的，则称它为线性微分方程，否则称之为非线性微分方程。

本论文主要介绍常微分方程，也简称微分方程。

以y错误！未找到引用源。为未知函数，x错误！未找到引用源。为自变量的一阶常微分方程的一般形式可表示为微分方程，

F(x,y,y,)0, (1.2.1) 将（1.2.1）中y,错误！未找到引用源。解出，则得到方程

y,f(x,y), (1.2.2)

或

M(x,y)dxN(x,y)dy0, (1.2.3)

也称（1.2.1）为一阶隐式微分方程，（1.2.2）为一阶显式微分方程，（1.2.3）为一阶微分方程的微分形式。

错误！未找到引用源。阶隐式方程的一般形式为

F(x,y,y,,,y(n))0, (1.2.4)

错误！未找到引用源。阶显式方程的一般形式

y(n)f(x,y,y,,,y(n1)), (1.2.5)

方程（1.2.4）中，如果函数错误！未找到引用源。对未知函数错误！未找到引用源。和它的各阶导数y,y,,,y(n1)错误！未找到引用源。都是一次的，则称其为线性常微分方程，否则，称其为非线性微分方程。

以y错误！未找到引用源。为未知函数，错误！未找到引用源。为自变量的n错误！未找到引用源。阶线性微分方程具有如下形式：

y(n)p1(x)y(n1)pn1(x)y,pn(x)yf(x), (1.2.6)

2.常微分方程的解

设函数y(x)是定义在区间错误！未找到引用源。上的错误！未找到引用源。阶可微导数。如果把错误！未找到引用源。代入方程（1.2.4）后能使其成为恒等式，即

Fx,x,'x,..,nx0,xI

则称错误！未找到引用源。是微分方程（1.2.4）在区间错误！未找到引用源。上的一个解。

例如，yekx是微分方程dyky0在,的一个解.ytanx是微分方程dx

y'1y2在区间,的一个解. 22

如果关系式Fx,y0决定的隐函数y(x)是方程（1.2.4）的解，则我们称Fx,y0是（1.2.4）的隐式解。例如一阶微分方程

xdxydy0

有隐式解

x2y2c0.

我们把含有n个相互独立的任意常数c1,c2,...cn的解

y(x,c1,c2,..,cn)

称为n阶微分方程（1.2.4）的通解。

在通解之中当一组任意常数确定时，所得到确定的解称为特解。例如，yc1cosxc2sinx是二阶线性方程y''y0的通解，而ysinx,ycosxsinx都是其特解，其中c1,c2是任意常数。一般地，方程的特解可由其通解中任意常

数取确定的常数导出，且方程的通解不一定表示方程的所有解。

3.常微分方程初值问题

为了确定微分方程的一个特解，我们可以给出这个微分方程的所满足的定解条件，常见的定解条件是初始条件，即方程（1.2.4）在某一点x0所满足的条件： yx0y0,y'x0y'0,...,ynx0yn0. （1.2.7） 微分方程（1.2.4）连同初始条件（1.2.7）一起称为初始值问题。

一阶常微分方程初值问题是求解函数y(x)，且满足 dyf(x,y),axb, dx （1.2.8）

y(a)y.

其中f(x,y)错误！未找到引用源。为已知函数，y为给定的初值。 定理1.2.1 假设f(x,y)错误！未找到引用源。在矩形区域

(x,y)x[a,b],y(,)内连续，且关于变元yLipschitz连续，即存在正常数L，使得对任意x[a,b]及，成立不等式 f(x,y1)f(x,y2)Ly1y2, (1.2.9) 其中常数错误！未找到引用源。称为Lipschitz常数,则处置问题存在惟一解yxC[a,b]，yx连续依赖于初值y.

由常微分方程的基本理论，我们有：

换句话说，若是yx如下初值问题 dyf(x,y),axb, dx (1.2.10)

y(a)y.

的解错误！未找到引用源。，则

limyxyx,x[a,b].错误！未找到引用源。 0

因此问题（1.2.8）是适定的。

定理1.2.2 当假定函数f(x,y)错误！未找到引用源。满足定理1.2.1中的

条件时，也即微分方程（1.2.8）的解y(x)错误！未找到引用源。是适定的。

综上，高阶常微分方程初值问题一般为 dnydydn1yf(x,y,,,n1),axb,dxndxdx i (1.2.11) dy(a)y,i0,1,,n1.i1dxi

其中y1,,yn错误！未找到引用源。为给定值。

引进新的变量函数

dk1y(x)yk(x),axb,k1,2,,n. (1.2.12) k1dx

则初值问题（1.2.11）化成了一阶常微分方程组初值问题

dy1dxy2,

,

dyn1axb (1.2.13) yn,dx

dyndxf(x,y1,,yn),

yi(a)yi,i1,2,,n.

通过求解（1.2.13）得到（1.2.11）的解y(x).

1.3 常微分方程常用数值解法的思路

一般说来，对于某些典型微分方程可以通过初等积分法求解，而对于复杂微分方程，需要得到解在若干个点上的近似值或者便于计算的近似表达式即可。本论文讲研究常微分方程常用数值解法及其matlab程序设计。

对于一阶常微分方程的初值问题为

dyf(x,y), (1.3.1) dxy(x0)y0.

求微分方程初值问题（1.3.1）的数值解，就是求解函yy(x)在一系列离

散点错误！未找到引用源。xi(i1,2,,n)错误！未找到引用源。上精确值y(x1),y(x2),,y(xn)错误！未找到引用源。的近似值y1,y2,,yn错误！未找到引用源。。

由于计算的复杂及繁冗，在计算常微分方程初值问题的数值解是一般使用Matlab软件进行编程计算，按如下步骤：

（1）引入点列xi错误！未找到引用源。，其中xixi1hi,i1,2,,h，错误！未找到引用源。称为步长。为了便于使用计算机进行编程计算，一般取步长为定值，即hih错误！未找到引用源。，

xix0ih,i0,1,， (2.2)

（2）使用常微分方程常用数值解方法，即输入由yi1错误！未找到引用源。计算出yi(i1,2,,n)错误！未找到引用源。的递推公式。

（3）利用（2）中的公式逐步求出近似解y1,y2,,yn.错误！未找到引用源。

第二章 单步法

所谓单步法是指这类方法在计算yn1时，只用到前一步的值xn1,xn,yn，然后逐步往下计算。单步法的精度不高，这是它的缺点。但我们只要给定初值就可以计算，所以还是比较方便的。以下我们介绍两种典型的单步法：Euler法和Runge-Kutta法。

2.1 Euler法

为了计算出初始值问题 dyf(x,y), dx （2.1.1）

y(x0)y0.

在[x0,x0b]这一区间上若干个点上解的近似解，我们将此区间n等分，令hb,xkx0kh,k1,2,...,n,yk为（2.1.1）的解，即xxk当时精确值y(xk)的n

近似解。我们的任务是要用尽可能简单的方法求出与y(xk)尽可能接近的值yk。Euler法又称Euler折线法，Euler折线法是数值解法最简单的一种。

Euler折线法的基本思想是利用微分中值定理对（2.1.1）的解yyx函数进行近似。由于yx是可微函数，故

yx0hyx0y'h,

其中，是介于x0和x0h之间的一个值。当h比较小时，y'和y'x0相差不大，故用x0代替上式中的，就得到了的近似值

y1y0fx0,y0h.

当y1知道以后，再用类似方法来求yx在点xk的近似值yk.

ykyk1fxk1,yk1h,k1,2,..,n. （2.1.2） 这样就得到了（2.1.1）的解在点x0,x1,...,xn的近似解，从几何上看，Euler折线法就是在局部范围内用切线上的值去替代解曲线上的值。

Euler公式在实际计算时较少采用，但由于它的结构简单，易于分析，在理论上具有非常重要的意义。

2.1.1 向前Euler法

Euler折线法的计算量小，但误差大，后来人们对Euler折线法进行了改进： 由方程（2.1.1），在节点xn处成立

y'xnfxn,yxn. （2.1.3） 将yxn1在xn处利用Taylor展开式展开，得

h2

y''n, yxn1y(xn)hy'xn2

其中将nxn,xn1.将式（2.1.3）代入上式得

h2

y''n. （2.1.4） yxn1y(xn)hyxn,yxn2

在式（2.1.4）中略去高阶项

h2

y''n （2.1.5） Exn,h 2

并用yn和yn1分别代替式（2.1.4）中的y(xn)和yxn1，可得Euler公式（也称

向前Euler公式：

yn1ynhfxn1,yn,n0,1,..,N1. （2.1.6）

2.1.2 向后Euler法

向后Euler法和向前Euler法相差不大，只是将yxn1在xn处利用Taylor展开式展开，得

h2

y''n, yxn1y(xn)hy'xn2

其中将nxn1,xn.将式（2.1.3）代入上式得

h2

y''n. （2.1.7） yxn1y(xn)hyxn,yxn 2

在式（2.1.7）中略去高阶项，并用yn和yn1分别代替式（2.1.7）中的y(xn)和yxn1，可得向后Euler公式：

yn1ynhfxn1,yn1,n0,1,..,N1. （2.1.8） 其中y0y(a)ya，局部截断误差为

h2

Exn,hy''n1O(h2). 2

因此，向后Euler公式同向前Euler公式，为一阶方法。

向前Euler公式（2.1.6）和向后Euler公式（2.1.8）相比，向前Euler公式关于yn1是显式格式。而向后Euler公式中，yn1隐含在方程中，称这样的格式为隐式格式。

2.1.3 中点差分公式法

若将展开式（2.1.4）和（2.1.7）多展开一项，得

h2h3

1y''ny"'(n), yxn1y(xn)hyxn,yxn26

h2h3

y''ny"'(n2). yxn1y(xn)hyxn,yxn26

1其中n,n2xn1,xn1.将上述两式相减并整理得

h3

yxn1y(xn1)2hyxn,yxny"'(n), 3

其中nxn1,xn1.略去高阶项，并用yn1,yn和yn1代替上式中y(xn1),yxn 和y(xn1)，即得中心差分公式

yn1yn12hfxn,yn,n0,1,..,N1. （2.1.9） 其局部截断误差为

h3

y"'nO(h3). Exn,h3

即中心差分公式为二阶方法。

2.1.4 梯形公式

如果将（2.1.1）中微分方程两端从xn到xn1积分，则得

b y(xn1)y(xn)xn1

xnfx,ydx.

如果用梯形求积公式计算上式中的积分，则有

hh3

f"(n,y(n)). yxn1y(xn)[fxn,yxnfxn1,yxn1212

略去高阶项可得梯形公式：

yn1ynh[fxn,ynfxn1,yn1],n0,1,..,N1. (2.1.10) 2

其中y0y(a)ya，局部截断误差为

h3

Exn,hy"'nO(h3). 12

显然，梯形公式为二阶隐式公式。

2.1.5 改进Euler法

一般来说，隐式格式比显示格式具有较好的数值稳定性。但是，由于在隐式格式里，yn1往往满足的是一个非线性方程，因而需要使用迭代法得到的一个yn1近似值，然后代入隐式格式作校正，并以这个校正值作为yxn1的近似值。印版地，可以选择适当的显示格式计算预测值，我们把这样的格式称为预估一校正公式。预估一校正公式往往既便于计算又具有较好的数值稳定性。

如，对于梯形公式，可用向前Euler公式计算出的预测值，再用梯形公式修正得到：

y0n1ynhf(xn,yn),



 h0

yn1ynf(xn,yn)f(xn1,yn1).

2



上述公式为改进Euler公式。

其总体截断误差为O(h3).在步长h取得足够小的时侯，我们用改进Euler方法解得数值与理论值非常接近。

2.2 Runge-Kutta法

Runge-Kutta法是一种间接使用Taylor展开式来获取高阶方法的数值方法，也是一种特殊的单步法，具有相当实用价值。其基本做法是：首先，用函数fx,y在xn,yn附近的m个点上的函数值的线性组合来代替yx的导数。然后，通过Taylor展开式，确定其中的组合系数。这样既可以避免计算高阶导数，又可以获得较高的方法阶。

设m为一个正整数（代表使用fx,y的函数值的个数）， cii1,2,...,m和

ai,biji2,3,...,m,1ji为一些特定参数，则数值方法



yyh(cKcK..cK),

n1122mmn1

K1f(xn,yn),

K2f(xna2h,ynhb21K1), (2.2.1)

......

m1



Kmf(xnamh,ynhbmjKj).

j1称为m级Runge-Kutta公式。

2.2.1 二阶Runge-Kutta法

显然，公式（2.2.1）是显式的。如果对于该代数方程组的一组解，对应的方法阶为p，则得到一个m级p阶显式Runge-Kutta公式。

当m2时，公式（2.2.1）变为

yn1ynh(c1K1c2K2),



K1f(xn,yn),

Kf(xah,yhbK).

n2n2112

其局部截断误差为

Exn,hO(h3). 经计算得参数方程组：



c1c21,

1

c2a2, （2.2.2）

2

1cb.2212

1

1.若取c1c2,a21,b211,则对应两二级二阶方法（又称变形Euler公式）：

2

h

yy(K1K2),nn1

2



K1f(xn,yn),

Kf(xh,yhK).

nn1

2

132

2.若取c1,c2,a2b21,则对应两级二阶方法（又称两级Huen公式）：

443

h

yy(K13K2),nn1

4



K1f(xn,yn),

22

K2f(xnh,ynhK1).

33

2.2.2 三阶Runge-Kutta法

1.三级Kutta公式（三级三阶格式）

h

yy(K14K2K3),nn1

6



K1f(xn,yn),



hhKf(x,yK),nn1

222

Kf(xh,yhK2hK).

nn123

2.三级Huen公式（三级三阶格式）

hyy(K13K3),nn1

4



K1f(xn,yn),

 hh

Kf(x,yK),nn12

33



Kf(x2h,y2hK).3nn233

这两个常用方法在实际使用中能达到较低的精确要求。

2.2.3 古典Runge-Kutta法

古典Runge-Kutta公式是精度高、最实用的一种常微分方程常用数值解法。



yn1

ynh(K12K22K3K4),

6

K1f(xn,yn), 

Kxhh

2f(n2,yn2K1),

K3

f(xnh2

,ynhK2),K4fxnh,ynhK3.

2.2.4 一级隐式中点公式

设m为一个正整数（代表使用fx,y的函数值的个数），ai,biji2,3,...,m,1ji为一些特定参数，则数值方法

yn1ynh(c1K1c2K2..cmKm),m

K1

f(xna1h,ynhb1jKj

),

j1



mK2f(x

na2h,ynhb2jKj), (2.2.2)

j1......

m1Kmf(xnamh,ynhbmjKj).

j1称为m级隐式Runge-Kutta公式。 当m1,p2时，有一级隐式中点公式为

yn1ynhK1,



Khh

1f(xn2,yn2K1).2.2.5 二级隐式中点公式

当m1,p2时，有二级隐式中点公式为

cii1,2,...,m和

h

yy(K1K2),nn1

2



3132

K1f(xnh,ynhK1hK2),

6412

33231Kf(xh,yhKhK2).2nn1

6124

相对显式公式而言，对于同样的m，隐式格式不仅可获取较高的精度，而且具有较好的数值稳定性。

第三章 多步法

所谓多步法指：这类方法在计算yn1时，除了用到前一步的值xn1,xn,yn，之外，还要用到

xnp,ynp(p1,2,,k;k0)， 这前面k步的值，这个算法的代表就是阿达姆斯(Adams)方法。

本章介绍线性多步法。线性多步法多使用于求解步骤较多的情况。一般地，

k步线性多步法公式可写成

yn1jyxnjhjfxnj,ynj. （3.1）

j0

j1

k1k1

当10时，式（3.1）为显式格式。否则为隐式格式。 1.显式格式

首先，将常微分方程（2.1.1）两端在区间[xnp,xn1]上求定积分，得 yxn1yxnp

xn1xnp

ft,ytdt.

如果在上面等式中的定积分用关于结点xn,xn1,...,xnq的数值积分近似，可得如下

显式格式

yn1ynphjfxnj,ynj.

j0q

其中

1xn1

jljtdt,j0,1,...,q,

hxnp

而ljtj0,1,...,q为函数fx,yx关于结点xn,xn1,...,xnq的Lagrange插值函数的基函数。

当q3时，四步Adams显式格式

yn1yn

h

[55fxn,yn59fxn1,yn137fxn2,yn29fxn3,yn3]. 24

25155

hynnO(h6). 720

其局部截断误差为

En1

2.隐式格式

如果在等式（3.1）中的定积分用关于结点xn1,xn,...,xn1q的数值积分近似，可得如下隐式格式

yn1ynphjfxn1j,yn1j.

j0q

其中

j

1xn1

ljtdt,j0,1,...,q, xnph

而ljtj0,1,...,q为函数fx,yx关于结点xn1,xn,...,xn1q的Lagrange插值函数的基函数。

当q3时，三步Adams隐式格式 yn1yn

h

[9fxn1,yn119fxn,yn5fxn1,yn1fxn2,yn2]. 24

1955

hynnO(h6). 720

其局部截断误差为

En1

Adams显式与隐式法相比较而言：同阶数情况下，隐式格式的局部截断误差

系数的绝对值比显式格式的小，隐式格式的稳定性相对显式格式的要大。

总结

本论文总结了单步法（向前Euler法，向后Euler法，中心差分法，梯形法，改进Euler法，Runge-Kutta法）及线性多步法（Adams法）几种比较常用的常 微分方程数值解法，并通过数值算例对同一问题的解的误差计算比较，总结各种算法的优缺点。一般而言，对Euler法而言，求解常微分方程初值问题数值时计算简单、程序简洁并且计算量小，但是，计算精度却很低。中心差分法及梯形法较好地提高计算精确度，同时又保持计算简单，程序简洁的优点。改进的Euler法比Euler法得到的计算结果要好很多。我们用改进Euler方法解得数值与理论值非常接近。对Runge-Kutta法来说，越高阶的方法计算精度越高，但相应的计算量也越大。古典Runge-Kutta公式是精度高、最实用的一种常微分方程常用数值解法。而线性多步法多使用于求解步骤较多的情况。

我们要把常微分方程常用数值解法掌握并灵活运用并加以推广应用，来解决我们的实际生活中所遇到问题，为我国的发展服务，充分体现了常微分方程常用数值解法的重大意义。作为一名数学工作者，我们有责任也有义务继续探索常微分方程组的数值解法，并且在不断探索其解法的前提下，开发出更好更优的软件来解决实际算法问题。但是，其现有理论也还远远不能满足需要，有待于进一步

研究，使这门学科的理论更加完善。

数值实例

1.分别用改进Euler公式和古典Runge-Kutta公式解常微分方程初值问题

2x

,0x1,yy y

y(0)1

取步长h0.1，计算结果保存于下表1

中，精确解为y(x) 解：fun函数：

function f=fun(x,y); f=y-2*x/y;

精确解：

x=0:0.1:1; y=sqrt(1+2*x) y =

Columns 1 through 9

1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125

Columns 10 through 11 1.6733 1.7321

1.应用改进Euler法，运用matlab软件，程序为： format long h=0.1; x=0:h:1.0; y=zeros(size(x)); y(1)=1;

y1=sqrt(1+2*x); N=length(x); for n=1:N-1

K1=feval(@fun,x(n),y(n));

K2=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K1); y(n+1)=y(n)+h*K1; end y

2.应用古典Runge-Kutta法，运用matlab软件，程序为： format long h=0.1; x=0:h:1.0; y=zeros(size(x)); y(1)=1;

y1=sqrt(1+2*x); N=length(x); for n=1:N-1

K1=feval(@fun,x(n),y(n));

K2=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K1); K3=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K2); K4=feval(@fun,x(n)+h,y(n)+h*K3);

y(n+1)=y(n)+h/6*(K1+2*K2+2*K3+K4); end y

表1

2.分别用向前Euler

公式、向后Euler公式、中心差分公式和梯形公式解常微分方程初值问题

yy,0x1, y(0)1

取步长h0.1。

解：M文件：

x=0:0.1:1;

y=-exp(x);

1.应用向前Euler法，运用matlab软件，程序为：

h=0.1;

x=0:h:1;

y=-exp(x);

y1=zeros(size(x));

n=length(x);

y1(1)=1;

for i=1:n-1;

y1(i+1)=(1+(h/2))*y1(i);

end

error=abs(y-y1)

2.应用向后Euler法，运用matlab软件，程序为：

h=0.1;

x=0:h:1;

y=-exp(x);

y2=zeros(size(x));

n=length(x);

y2(1)=1;

for i=1:n-1;

y2(i+1)=(1/(1-(h/2)))*y2(i);

end

error=abs(y-y2)

3.应用中心差分法，运用matlab软件，程序为：

h=0.1;

x=0:h:1;

y=-exp(x);

y3=zeros(size(x));

n=length(x);

y3(1)=1;

y3(2)=-exp(h);

for i=1:n-1;

y3(i+2)=y3(i)+h*y3(i+1);

end

error=abs(y-y3)

4.应用梯形法，运用matlab软件，程序为：

h=0.1;

x=0:h:1;

y=-exp(x);

y4=zeros(size(x));

n=length(x);

y4(1)=1;

for i=1:n-1;

y4(i+2)=(8/(4-h)-1)*y4(i);

end

error=abs(y-y4).

3.分别用一级隐式中点公式、二级隐式中点公式解常微分方程初值问题

2xyy,0x1,yy(0)1 

取步长h0.1，计算结果保存于下表1

中，精确解为y(x)

解：fun函数：

function f=fun(x,y)

f=y-2*x/y;

1.应用一级隐式中点公式法，运用matlab软件，程序为：

format long

x=0:0.1:1.0;

y=sqrt(1+2*x);

h=0.1;

x=0:h:1.0;

y1=zeros(size(x));

y1=sqrt(1+2*x);

y1(1)=1;

N=length(x);

for n=1:N-1

error=1;

K1=feval(@fun,x(n),y1(n));

while(error>=1.0e-010)

K11=feval(@fun,x(n)+0.5*h,y1(n)+0.5*h*K1);

error=abs(K11-K1);

K1=K11;

end

y1(n+1)=y1(n)+h*K1;

end

e=abs(y-y1)

2.应用二级隐式中点公式法，运用matlab软件，程序为：

format long

x=0:0.1:1.0;

y=sqrt(1+2*x);

h=0.1;

x=0:h:1.0;

y2=zeros(size(x));

y2=sqrt(1+2*x);

y2(1)=1;

N=length(x);

for n=1:N-1

error=1;

K1=feval(@fun,x(n),y2(n));

K2=feval(@fun,x(n),y2(n));

while(error>=1.0e-010)

K11=feval(@fun,x(n)+(3-sqrt(3))/6*h,y2(n)+1/4*h*K1+(3-2*sqrt(3))/

12*h*K2);

K22=feval(@fun,x(n)+(3+sqrt(3))/6*h,y2(n)+1/4*h*K2+(3+2*sqrt(3))/

12*h*K1);

error=max(abs(K11-K1),abs(K22-K2));

K1=K11;

K2=K22;

end

y2(n+1)=y2(n)+0.5*h*(K1+K2);

end

e=abs(y-y2)

4.分别用Adans显式公式和Adans隐式公式解常微分方程初值问题

yyx22,0x1, y(0)1

取步长h0.1。

解：fun函数：

function f=fun(x,y);

f=-y+x.^2;

1.应用Adans显式法，运用matlab软件，程序为：

format long

h=0.1;

x=0:h:1.0;

y=zeros(size(x));

y1=x.^2-2*x+4-3*exp(-x);

y(1)=1;

N=length(x);

for n=1:N-1

K1=feval(@fun,x(n),y(n));

K2=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K1);

K3=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K2);

K4=feval(@fun,x(n)+h,y(n)+h*K3);

y(n+1)=y(n)+h/6*(K1+2*K2+2*K3+K4);

end

for n=5:N-1

M1=feval(@fun,x(n-1),y(n-1));

M2=feval(@fun,x(n-2),y(n-2));

M3=feval(@fun,x(n-3),y(n-3));

M4=feval(@fun,x(n-4),y(n-4));

y(n)=y(n-1)+h/24*(55*M1-59*M2+37*M3-9*M4);

end

y=abs(y-y1)

2.应用Adans隐式法，运用matlab软件，程序为：

format long

h=0.1;

x=0:h:1.0;

y=zeros(size(x));

y1=x.^2-2*x+4-3*exp(-x);

y(1)=1;

N=length(x);

for n=1:N-1

K1=feval(@fun,x(n),y(n));

K2=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K1);

K3=feval(@fun,x(n)+0.5*h,y(n)+0.5*h*K2);

K4=feval(@fun,x(n)+h,y(n)+h*K3);

y(n+1)=y(n)+h/6*(K1+2*K2+2*K3+K4);

end

for n=3:N-1

error=1

M2=feval(@fun,x(n),y(n));

M3=feval(@fun,x(n-1),y(n-1));

M4=feval(@fun,x(n-2),y(n-2));

temp1=y(n);

While(error>=1.0e-013)

M1=feval(@fun,x(n+1),temp1);

tempt2=y(n)+h/24*(9*M1+19*M2-5*M3+M4);

error=abs(temp1-temp2);

temp1=temp2;

end

y(n+1)=temp2;

end

y=abs(y-y1).

参考文献

[1] 戴嘉尊，邱建贤，微分方程数值解.南京：东南大学出版社，2002.

[2] 冯康. 数值计算方法.北京:国防工业出版社，1978.

[3] 周义仓，常微分方程及其应用.北京，科学出版社.2009.

[4] 李瑞遐，何志庆.微分方程数值解法.上海：华东理工大学出版社.2005.

[5] 徐树方，高平，张平文.数值线性代数.北京：北京大学出版社.2000.

[6] 余徳浩，汤华中.微分方程数值解法.北京：高等教育出版社.2003.

[7] 曾金平.微分计算方法.长沙：湖南大学出版社.2006.

[8] 谭浩强，C程序设计.北京：清华大学出版社.2009.

[9] 陆金甫，关治．偏微分方程数值解法.北京：清华大学出版社.1987.

[10] 王高雄，周之铭，朱思铭，王寿松.常微分方程.北京：高等教育出版社.

1987.

[11] 胡健伟，汤怀民.微分方程数值方法.北京:科学出版社.2000.

[12] 李信真，车刚明，欧阳洁，封建湖.计算方法.西安:西北工业大学出版社.

2000.

[13] 玛康.数值计算方法.杭州:浙江大学出版社.2003.

[14] 李荣华，玛果忱.微分方程数值方法(第三版).北京:高等教育出版社.1996.

[15] 刘会灯，朱飞.MATLAB编程基础与典型应用.北京:人民邮电出版社.2006.

[16] 苏金明，阮沈勇.MATLAB实用教程（第二版）.北京:电子工业出版社.2008.

[17] 刘卫国.MATLAB程序设计教程.北京:中国水利水电出版社.2005.

[18] 马知恩，周义仓.常微分方程定性与稳定性方法.北京:科学出版社.2007.

[19] C.W.吉尔.常微分方程初值问题的数值解法,费景高译.北京:科学出版社. 1978.

[20] 楼顺天，于卫，闰华梁.MATLAB程序设计语言西安:西安电了科技大学出版

社.1997.

[21] 李大侃，常微分方程数值解.浙江:浙江大学出版社.1994.

[22] 李立康，龄崇华，朱政华.微分方程数值解法.上海:复旦大学出版社.2003.

[23] 李庆扬，关治，白峰杉.数值计算原理.北京:清华大学出版社.2000.

[24] 李庆扬，王能超，易大义.数值分析.北京:清华大学出版社，施普林格

出版社.2001.

[25] 丁同仁，李承治.常微分方程.北京：高等教育出版社.1991.

[26] Hackbusch W. Iterative Solution of Large Sparse Systems of Equations.

New York:Springer-Verlag.1994

[27] Kincaid D,Cheney W. Numerical Analysis:Mathematics of Scientific

Computing.北京：机械工业出版社.2003.

[28] Kress R. Numerical Anaysis.New York:springer-Verlag.1998.

[29] Larson S,Thomee V. Partial Differential Equations with Numercial

Methods.Berlin:Springe-Verlag.2003.

[30] Quarteroni A,Valli A. Domain Decomposition Methods for Partial

Differential Equations.Oxford:Clarendon Press.1999.

[31] Lambert,J. D. Numerical Methods for Ordinary Differential Systems.

London: Wiley.1991.

[32] R. D.Richtmyer,K. W. Morton:Difference Methods of Initial value

Problems.1976.

[33] Butcher,J. C. The Numerical Analysis of Ordinary Differential

Equations.New York:John Wiley.1987.

[34] M. Braun. Differential Equations and Their Applications(Fourth edition).New York:Springer-Verlag.1993

[35] Sanz-Serna,J M and Ca1vo,M P. Numerical Hamiltonian Problems. London: Chapman & Hal.1994.