微分方程数值解

WORD格式整理

专业资料 值得拥有 微分方程数值解及其应用

绪论

自然界中的许多事物的运动和变化规律都可以用微分方程来描述，因此对工程和科学技术中的实际问题的研究中, 常常需要求解微分方程．但往往只有少数较简单和典型的微分方程可求出其解析解，在大多数情况下,只能用近似法求解，数值解法是一类重要的近似方法．本文主要讨论一阶常微分方程的初值问题的数值解法，探讨这些算法在处理来自生活实际问题中的应用，并结合MATLAB软件，动手编程予以解决．

１ 微分方程的初值问题[1]

1.1 预备知识

在对生活实际问题的研究中，通常需要考虑一阶微分方程的初值问题

00(,)()dyfxydxyxy （1）

这里,fxy是矩形区域R：00,xxayyb上的连续函数．

对初值问题（1）需要考虑以下问题：方程是否一定有解呢?若有解，有多少个解呢?下面给出相关的概念与定理．

定义1 Lipschitz条件[1][2]：矩形区域R：00,xxayyb上的连续函数,fxy若满足：存在常数0L，使得不等式1212,,fxyfxyLyy对所有12,,,xyxyR都成立，则称,fxy在R上关于y满足Lipschitz条件.

定理 1 解的存在唯一性定理[1][3]：设f在区域,,DxyaxbyR上连续，关于y满足Lipschitz条件，则对任意的00,,xabyR，常微分方程初值问题（1）当,xab时存在唯一的连续解yx． WORD格式整理

专业资料 值得拥有 该定理保证若一个函数,fxy关于y满足Lipschitz条件，它所对应的微分方程的初值问题就有唯一解.在解的存在唯一性得到保证的前提下，自然要考虑方程的求WORD格式整理

专业资料 值得拥有 解问题．求解微分方程虽然有多种解析方法，但根据工程和科学实践问题所得到的微分方程往往很复杂，在很多情况下不能或很难给出解析解，有时即使能求出形式解，也往往因形式过于复杂或计算量太大而不实用，因此从实际问题中归结出来的微分方程主要依靠数值解法．

定义 2 微分方程数值解：对初值问题（1）寻求数值解就是寻求解yx在一系列离散节点

0121nnxxxxx

上的近似解0121,,,,,,nnyyyyy，相邻两个节点的间距1nnnhxx称为步长.在一般情况下假定0,1,ihhi为常数，这时节点为0,0,1,2,nxxnhn．

要求微分方程数值解，首先要建立数值算法，即对初值问题（1）中的方程离散化，建立求解数值解法的递推公式．一类是计算1ny时只用到前一点的值ny，称为单步法；另一类是用到1ny前面k点的值11,,,nnnkyyy称为k步法.

对初值问题（1）式的单步法可用一般形式表示为

11(,,,)nnnnnyyhxyyh，

其中多元函数与,fxy有关，当含有1ny时，方法是隐式的；若中不含1ny，则为显式方法，所以显式单步法可表示为

1(,,)nnnnyyhxyh. （2）

设yx是初值问题（1）的准确解，称11(,,)nnnnnTyxyxhxyxh为显式单步法（2）的局部截断误差. 若存在最大正整数p,使显式单步法（2）式的局部截断误差满足11,,pnTyxhyxhxyhOh，则称（2）式有p阶精度.

1.2几种常用的数值解法及其分析、比较

1.2.1欧拉法与后退欧拉法 WORD格式整理

专业资料 值得拥有 1）欧拉法：欧拉曾简单地用差分代替微分，即利用公式

1()()()(,)nnnnnyxyxyxfxyh

将初值问题（1）离散化，则问题（1）可化为

1(,),nnnnyyhfxy0nxxnh， （3）

此方法称为欧拉法.

欧拉方法的几何意义在数值计算公式中体现了出来.在xy平面上，一阶微分方程的解yyx称作它的积分曲线.积分曲线上一点,xy的切线斜率等于函数,fxy，按函数,fxy在xy平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.

基于上述几何解释，从初始点000(,)Pxy出发，先依方向场在该点的方向上推进到1xx上一点1P，再从1P依方向场的方向推进到2xx上一点2P，循环前进便作出一条折线012PPP，因此欧拉方法又称为折线法.若初值0y已知，则由（3）式可逐步算出

10002111(,),(,),yyhfxyyyhfxy

为了分析计算公式的精确度，通常可用泰勒展开将1nyx在nx处展开，则有2''11,,.2nnnnnnnnhyxyxhyxyxhyxx

在nnyyx的前提下，,,.nnnnnfxyfxyxyx可得欧拉法（3）的误差为

2211.22nnnnhhyxyyyx

容易看出，欧拉法（3）式具有一阶精度.

2）向后欧拉方法：如果对微分方程（1）从nx到1nx积分，得 WORD格式整理

专业资料 值得拥有 11,nnxnnxyxyxftytdx， （4）

如果（4）式右端积分用右矩形公式11,nnhfxyx近似，则得到另一个公式

111,nnnnyyhfxy， （5）

称为后退欧拉法.

值得一提的是:后退欧拉法与欧拉公式有着本质的区别，后者是关于1ny的直接计算公式，它是显式的，而（5）式的右端含有关于1ny的表达式，它是隐式的.在利用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化.具体迭代过程如下：

首先利用欧拉公式(0)1(,)nnnnyyhfxy给出迭代初值(0)1ny，把它代入（5）式的右端，使之转化为显式，直接计算得11(1)(0)1(,)nnnnyyhfxy.如此反复进行，得

(1)()111(,)kknnnnyyhfxy0,1,k ，

则得到后退欧拉法的迭代公式

(0)1(1)()111(,)(,)nnnnkknnnnyyhfxyyyhfxy，

可以看出，后退欧拉法具有一阶精度，且计算比较麻烦.

1.2.2梯形方法

为得到比欧拉法精确度高的计算公式，在等式（4）式右端积分中若用梯形求积公式近似，并用ny代替nyx，1ny代替1nyx，则得

111,,2nnnnnnhyyfxyfxy， （6）

称其为梯形方法.

梯形方法与后退欧拉法一样，都是隐式单步法，可用迭代法求解，其迭代公式为 WORD格式整理

专业资料 值得拥有 (0)1(1)()111(,),,2nnnnkknnnnnnyyhfxyhyyfxyfxy . （7）

为了分析梯形公式的收敛性，将（6）与（7）式相减，得

(1)()111111,,2kknnnnnnhyyfxyfxy，0,1,2,k

因为,fxy满足Lipschitz条件，于是有(1)()11112kknnnnhLyyyy，其中L为,fxy关于y的Lipschitz常数.如果选取h充分小，使得12hL，则当k时有(1)11knnyy，这说明迭代过程（7）式是收敛的[4].容易推导得出梯形法（7）式是二阶方法.

经分析，梯形方法虽然提高了精度，但是以增加计算量为代价的.从上述的迭代公式可以看出，每迭代一次都要重新计算,fxy的值，而且迭代又要进行若干次，计算相当的复杂.为此，有没有比较简便的计算方法呢？下面给出改进的欧拉方法.

1.2.3改进的欧拉方法

由前面的讨论可知，梯形法计算相对复杂，现对上面的梯形法进行简化，具体方法是只计算一两次就转入下一步的计算，先用欧拉公式（3）求得一个初步的近似解1ny，称为预测值，再利用公式（6）把它校正一次，这样建立的预测-校正系统通常称为改进的欧拉公式.具体公式如下

11,,,2nnnnnnnnhyyfxyfxyhfxy （8）改进的欧拉法与梯形法一样，是二阶方法.

1.2.4 Runge-Kutta方法

由前面讨论可知，从（4）式

11,nnxnnxyxyxfxyxdx

可以看出，若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必WORD格式整理

专业资料 值得拥有 然要增加求积积点，为此将（4）式的右端用求积公式表示为

11,,nnrxininixifxyxdxhcfxhyxh， （9）

一般来说，点数r越多，精度越高，上式右端相当于增量函数,,xyh，为得到便于计算的显式方法，将公式（9）表示为：

1,,,nnnnyyhxyh （10）

其中

1111,,,,,2,rnniiinniininijjjxyhcKKfxyKfxhyhKir （11）

这里,,iiijc均为常数. ic为加权因子，iK为第i段斜率，共有r段.我们把（10）和（11）称为r级显式Runge-Kutta法，简称为R-K方法.下面给出其中最经典最常用的一个公式：

11234121324322,6,,,22,,22,.nnnnnnnnnnhyyKKKKKfxyhhKfxyKhhKfxyKKfxhyhK （12）

Runge-Kutta方法作为一种重要的单步方法，具有很高的实用价值，它关于初值是稳定的，其解连续地依赖于初值，是一类便于应用的单步法，为了计算1ny，只用到前面一步的值ny即可，因此每步的步长可以独立取定.常用的Runge-Kutta方法精度较高，为了达到预定的精度，与欧拉方法与梯形法相比，步长h可取得大些，求解区间WORD格式整理