chap常微分方程的数值解法实用

会计学1chap常微分方程的数值解法实用2第7章常微分方程（组）的数值解法

在科学研究和工程实践中会遇到很多微分方程，虽然从理论上可以证明其解的存在性，但其解的解析表达式往往是很难求解的，或者即使可以写出来，但也难于计算，此时，只能借助数值解来解决问题.常微分方程（组）定解问题是自然科学和工程技术领域中常见的数学模型.本章介绍求解此类问题的基本理论和数值解法。

第1页/共55页3定义7.0.1若存在常数L>0,使得对一切的

x∈[a,b]及y,,均有

则称f(x,y)在

D上关于

y

满足

Lipschitz条件,其中L

称为

Lipschitz常数.

我们首先考虑一阶常微分方程初值问题其中f(x,y)是区域上的实值函数.(7.0.1)第2页/共55页4

我们首先给出常微分方程初值问题解的存在惟一性定理。定理假设f(x,y)∈C(D),且关于

y

满足

Lipschitz条件,则一阶常微分方程初值问题(7.0.1)存在唯一解.

下面在此前提下,我们讨论上述初值问题(7.0.1)的数值解法。第3页/共55页5然后在节点上建立逼近于原初值问题的计算格式(或差分格式),由此计算出原问题的解y(x)

在节点x1,x2,...,

xN

处的近似值：

y1,y2,...,

yN，

称它们为常微分方程初值问题的数值解.

相邻两个节点的距离hn=xn+1-xn称为步长,

通常取定步长h>0,即节点xn=x0+nh,n=0,1,…

,

N.

其基本思想是在区间[a,b]上引入一系列节点第4页/共55页67.1初值问题计算格式的建立1.数值微分方法在等距节点下讨论问题.利用两点数值微分公式7.1.1计算格式的建立将上式代入初值问题(7.0.1)有(7.1.1)第5页/共55页7略去余项,并以数值解yn,

yn+1代替y(xn)及y(xn+1),则得差分方程上式称为Euler公式.利用此式可由初值y0

出发按“步进式”方法,逐步求得数值解y1,y2,...,

yN.由于计算yn+1时,只用到它前一步的结果yn,这类公式称为单步法.又因为其关于yn+1是显式形式,故称该Euler公式为显格式.第6页/共55页8如果利用下列数值微分公式由

类似的可导出上述公式称为后退的Euler公式,此公式为单步法公式.又因为它关于yn+1成隐式形式,所以该公式为隐式公式，简称隐格式.第7页/共55页9类似地，可导出上述公式称为Euler两步法公式.这因为，当计算yn+1时,要用到yn-1与yn

.显然它也是显格式.如果利用下列三点数值微分公式(7.1.3)

第8页/共55页10设y(x)∈C2[a,b],由Taylor公式有由于故上式即为

略去余项,并以yn,

yn+1代替y(xn)及y(xn+1),得到的差分方程正是Euler公式.(7.1.4)2.Taylor展开法第9页/共55页113.数值积分方法对，在区间[xn,xn+1]上积分，得则有对上式中的积分采用不同的数值积分公式可得到不同的差分方程.例如,对上式的积分采用左矩形公式,可得到Euler公式.第10页/共55页12若对此式的积分采用梯形公式,

则有若略去余项,以yn,

yn+1代替y(xn)及y(xn+1),得到的差分方程第11页/共55页13上式称为梯形公式.由于它关于yn+1

成隐式形式,故其为隐格式.隐格式求解比较困难,当yn

已知时,要求yn+1

，需解关于yn+1

的非线性方程.在实际应用时,上式常与Euler公式联合使用,构成如下计算格式:(7.1.6)第12页/共55页14隐式梯形公式的迭代格式(7.1.7)由上式可以得到一个序列:,k=0,1,…,关于此序列的收敛性,有如下的定理.第13页/共55页15定理7.1.1

设f(x,y)在区域

D上关于

y

满足

Lipschitz条件,即其中L为Lipschitz常数,当步长时,对任意的初值按格式（7.1.7）生成的序列收敛于梯形公式（7.1.6）的解

.第14页/共55页16

为了减少计算量,可采用预测-校正格式.方法是先用Euler公式求得一个初始近似值

称为预测值,再把带入梯形公式右端计算一次求得yn+1

称之为校正值,即预测:校正:上式称为预测-

校正公式或改进的

Euler

公式.上式也可写成如下形式:第15页/共55页17例7.1.1：利用Euler公式与改进的Euler公式求解初值问题（步长h=0.1）解：由步长h=0.1，知节点设数值解为利用Euler公式得第16页/共55页18计算结果见下表(见书P227表7.1)

此初值问题的解析解为,从上表可以看出,数值解yn与解析解y(xn)比较,yn精度较差.第17页/共55页19解此问题的改进的Euler公式为

同Euler公式比较,改进的Euler法显然精度提高了.由于误差大小是评价计算格式优劣的重要依据,故需要给出有关误差的概念.计算结果见下表(见书P228表7.2)第18页/共55页207.1.2

截断误差与方法的精度定义7.1.1

称误差

en+1=y(xn+1)-yn+1为数值方法在点xn+1的截断误差,又称整体截断误差.

设

yk=y(xk

)(k=0,1,...,n)，则

为数值方法在点xn+1的局部截断误差.第19页/共55页21整体截断误差en+1

是在没有引进舍入误差的情况下,纯粹因为不准确的计算格式造成的,故又称为方法误差.它不仅与x=xn+1这一步的计算有关,而且和xn,xn-1,...,x1

这几步的计算都有关系.

局部截断误差是假设xn之前各数值解没有误差,仅由xn到xn+1这一步计算由计算格式引起的误差.第20页/共55页22如Euler公式在点xn+1的整体截断误差en+1=y(xn+1)-yn+1局部截断误差第21页/共55页23定义7.1.2若某数值方法的局部截断误差为

则称该方法具有P

阶精度,或称其为P

阶方法.可以证明:Euler方法的局部截断误差其具有一阶精度.梯形方法的局部截断误差其具有二阶精度.改进的Euler方法的局部截断误差具有二阶精度.第22页/共55页247.2Runge-Kutta方法

继续讨论前面的Taylor展开法。设y(x)∈C2[a,b],由Taylor公式有

由故上式即为

略去余项,并以yn,

yn+1代替y(xn)及y(xn+1),得到Euler公式.第23页/共55页25进一步假设设y(x)∈Cp+1[a,b],由Taylor公式有其中

由故式(7.2.1)即为第24页/共55页267.2Runge-Kutta方法略去余项,并以数值解yn,

yn+1代替(7.2.3)中的解析解y(xn)及y(xn+1),可得到一个差分方程，即其中余项可写成注：这里第25页/共55页27在(7.2.3)中略去余项,用yn,

yn+1代替y(xn)及y(xn+1)其中称(7.2.4)式为求解常微分方程初值问题数值解Taylor的格式.第26页/共55页28由于局部截断误差

可知它是一个p阶方法。当p=1时，上式正是Euler公式。但当p>2时，需要计算f(x,y(x))的高阶导数，特别是对于复杂函数f(x,y(x))的求导，这无疑是大大增加计算量，这是它最大的缺点。因此高阶的Taylor方法是不实用的。德国数学家C.Runge及M.W.Kutta提出了一种改进策略，得到了至今还被作为高精度的单步法广泛使用龙格-库塔法(Runge-Kuttamethod)。

第27页/共55页297.2.1Runge-Kutta方法的基本思想

Runge-Kutta方法是利用f

在某些点处函数值的线性组合替代（7.2.4）步长h后面括号中的因子来构造差分方程,从而避免了高阶导数的计算,这就是Runge-Kutta方法的基本思想.

用f

在某些点处函数值的线性组合替代这一部分第28页/共55页30其一般形式为:其中r是上式中调用

f

的个数,r称为级数，为待定参数,适当确定这些参数,可使之具有尽可能高的精度.如局部截断误差满足第29页/共55页31

7.2.2

二阶Runge-Kutta

方法考虑r=2的情况,此时有利用二元函数的一阶Taylor公式，即全微分公式

希望适当选择参数使上式的局部截断误差为即为二阶方法.,下面将yn+1与y(xn+1)作比较第30页/共55页32从而有将上式代入再由第31页/共55页33得到在下面要将yn+1与y(xn+1)作比较，使它们的局部截断误差满足为此考虑y(xn+1)。第32页/共55页34再根据y(xn+1)在点xn的一元3阶Taylor展开式由刚才已得到的让它们满足第33页/共55页35即由左式含有四个未知元三个方程,因此解不唯一.参数满足左式的一族公式统称二阶Runge-Kutta公式.

可得参数应满足下列方程组:第34页/共55页36

取上式称为中点公式.

取上式称为Heun公式.

可见,二阶Runge-Kutta公式,每计算一步需要两次调用f

的函数值.取得这正是改进的Euler公式.

第35页/共55页37

7.2.3

四阶Runge-Kutta方法当r=4时,类似地可导出四阶Runge-Kutta公式,这种公式也有一族,其中常用地有:标准(经典)的Runge-Kutta方法第36页/共55页38

Gill公式

Gill公式是标准的

Runge-Kutta公式的改进形式,这种算法可节省存储单元,并能控制舍入误差的增长.四阶Runge-Kutta公式,每一步计算需四次调用f的函数值,计算量较大,但其局部截断误差可达O(h5),精度较高.第37页/共55页39例7.2.1用标准的四阶Rung-Kutta法解初值问题,取步长h=0.2.解:

解此问题的计算公式为xn0.20.40.60.81.0ynyn-y(xn)1.18321.34171.48331.61251.73210.00000.00000.00000.00010.0001计算结果如下:显然在计算量大致相同的情况下,标准的Runge-Kutta方法比改进的

Euler方法精确度更高.

（参见p227和p228的表）第38页/共55页407.5一阶常微分方程组与高阶方程初值问题的数值解法

7.5.1一阶常微分方程组初值问题

第39页/共55页41写成向量形式：

其中注意:在形式上(7.5.2)与(7.0.1)一样,所以可以把求解常微分方程初值问题的各种数值方法推广到方程组上来.第40页/共55页42利用向量值函数的微积分理论,很容易推导出一阶常微分方程组初值问题的数值解法.

如Euler公式其中

第41页/共55页43(7.5.3)的分量形式为

或

第42页/共55页44四阶标准的Runge-Kutta公式

设则(7.5.5)的分量形式为第43页/共55页45四阶标准的Runge-Kutta公式的分量形式

其中第44页/共55页46例7.5.1：试写