常微分方程数值解法

第五章常微分方程数值解法§1引言在实际工作中常常会遇到求解常微分方程的定解问题。虽然我们已经学过一些常微分方程的定解问题的解法，但是那只是对一些特殊类型的方程给出了解析方法。而实际问题中归纳出来的常微分方程往往不能用解析方法来求解。有时虽然能求出其解析式子，但真要求它在某一点的值时，还会遇到一些麻烦。见下例。例已知初值问题0)0(21yxyy求其解。易知，dxeexyxtx022)(。若要求它在某一点的值，还要求积分，而这个积分是不能求出的，需要用数值积分来求解。因此不如直接用数值解法求其解了。以下我们就研究常微分方程的数值解法。本章主要考虑下列一阶方程的定解问题（又称初值问题）。00)(),(yxyyxfy（1-1）这里假设函数),(yxf满足Lipschitz（李卜西兹）条件。即),(yxf满足：yyLyxfyxf),(),(其中，L为某一常数。这是保证问题（1-1）有唯一解。下面就研究关于问题（1-1）的数值解法。所谓数值解法，就是直接寻求解函数)(xy在一系列离散点121nnxxxx上的近似值1y，2y，…，ny，1ny，…。其中，iiihxx1称为步长，今后如果不特殊说明，我们总是假定是等步长的。即有hxxii1，此时节点nhxxn0，,2,1,0n。由于00)(yxy为已知，所以自然设想利用这个已知信息求出)(1xy的近似值1y，然后由1y求得)(2xy的近似值2y，如此继续下去，这就是初值问题数值解法的一般思想。称为“步进法”。§2Euler方法（折线法）2-1Euler公式这一方法是初值问题数值解中最简单的一个方法，其精度不高。所以实际计算中很少应用。但在某种程度上反映了数值方法的基本思想，且在此基础上得到的某些改进的方法目前还在使用，因此我们有必要介绍一下这种方法。Euler公式的推导方法很多，在此我们用Tarlor展开的方法。∵)()(1hxyxynn)(!2)()(2nnnyhhxyxy其中，),(1nnnxx，由（1-1）式知，)(!2)](,[)()(21nnnnnyhxyxhfxyxy,2,1,0n当h充分小时，略去)(!22nyh，即得)](,[)()(1nnnnxyxhfxyxy取近似，并写成等式得，),(1nnnnyxhfyy,2,1,0n（2-1）此即称为Euler公式。由于略去的是)(!22nyh，可见这就是误差。但这只是在计算第n步时产生的误差，并非是从一开始到现在所产生的误差。因此我们称这个误差为局部截断误差。即在假定)(nnxyy的假设下，计算)(1nxy时产生的误差11)(nnyxy，称为局部截断误差。由此知Euler公式的局部截断误差约为)(!22nxyh。2-2后退的Euler公式后退的Euler公式的推导同样也有许多方法。在此，我们用差商代替导数的方法。因为nnnnnxxxyxyxy111)()()(假设)(nnxyy，又有（1-1）知))(,()(111nnnxyxfxy，取近似整理得，),(111nnnnyxhfyy,2,1,0n（2-2）此即称为后退的（隐式）Euler公式。公式（2-2）在计算1ny时要用1ny。这样的公式称为隐式公式，而象公式（2-1）的式子称为显式公式。对于隐式公式，在计算时要先给1ny提供一个初值，然后再用给定的公式开始计算，通常称为迭代法。对（2-2）式的迭代公式如下，),(),()(11)1(1)0(1knnnknnnnnyxhfyyyxhfyy,1,0k（2-3）直到)(1)1(1knknyy为止。以下研究式（2-2）的局部截断误差。∵),(111nnnnyxhfyy)(1nnxyhy又hyxyhxyxynnnn)()()()(1代入上式得，21)()(hyxyhyynnnn又有)(!2))(,()()(21nnnnnyhxyxhfxyxy下式减上式，且注意到已假设)(nnxyy，故得，11)(nnyxy=)(22nyh∴后退的Euler公式的局部截断误差约为)(!22nxyh。2-3梯形公式今考察以上两个公式，可见它们的局部截断误差相差一个负号，即有),(1nnnnyxhfyy)()(!232hoxyhn),(111nnnnyxhfyy)()(!232hoxyhn以上两式相加除以2，则可以消去)(!22nxyh，从而得到)],(),([2111nnnnnnyxfyxfhyy)(3ho略去)(3ho，则得一公式，)],(),([2111nnnnnnyxfyxfhyy（2-4）此称为梯形公式。其几何意义如图。注意到（2-4）式为隐式公式，同样需要迭代法来求解。其迭代公式为，)],(),([2),()(11)1(1)0(1knnnnnknnnnnyxfyxfhyyyxhfyy,1,0k（2-5）以下分析迭代公式（2-5）的收敛性。（2-4）减去（2-5）得),(),(2)(1111)1(11knnnnknnyxfyxfhyy由于已假设),(yxf满足Lipschitz（李卜西兹）条件。即),(yxf满足：yyLyxfyxf),(),(所以有)(11)1(112knnknnyyLhyy，取12Lh，则有)(11)1(11knnknnyyyy从而收敛（带条件收敛）。2-4改进的Euler公式由（2-5）式知，用梯形公式计算时要反复求函数值，因而工作量是很大的。而改进的Euler公式是指在用梯形公式迭代时只迭代一次便取作1ny，这就是，预测：),(1nnnnyxhfyy校正：)],(),([2111nnnnnnyxfyxfhyy（2-6）以上称为预测-校正系统。也称为改进的Euler公式为便于上机常采用以下形式，)(21),(),(11cpnpnncnnnpyyyyxhfyyyxhfyy（2-7）例（见教材）在以上推导后退的Euler公式时，我们用差商代替了导数。但若改用中心差商代替导数就有)(2)()(11nnnxyhxyxy取近似并整理得，),(2)()(11nnnnyxhfxyxy从而得，),(211nnnnyxhfyy用此公式作为预测值，用梯形公式作为校正值，便得到以下的预测-校正系统，预测：),(211nnnnyxhfyy校正：)],(),([2111nnnnnnyxfyxfhyy（2-8）（2-8）与（2-6）相比，突出的特点是它们有相同的精度。不过（2-8）式是两步公式，而前面几式是单步公式。两步以上的公式称为多步法，多步法不能自开始（自起步），需要借助其他单步法才能开始进行。而单步法可以自开始。§3Runge-Kutta方法在§1中，我们曾经用Taylor展开的方法推导了Euler公式。那时我们略去了)(2ho项，自然想到如果多取几项，也就是略去)(1pho项（p取的更大些）应该得到精度更高的公式。这就是Taylor级数法。它的一般公式如下，)(21!!2pnpnnnnyphyhyhyy（3-1）（3-1）式的局部截断误差为11)(nnyxy=)()!1()1(1nppyph),(1nnnxx一般地有以下定义，定义若一种方法的局部截断误差为)(1pho，则称该方法具有p阶精度。由（3-1）知，当p=1时，（3-1）式即为Euler公式。显然Euler公式具有1阶精度。同理，后退的Euler公式也具有1阶精度；而梯形公式具有2阶精度。由（3-1）知，p越大精度越高。但是它要用到f(x,y)的高阶导数，而求f(x,y)的高阶导数是不容易求出的。因此，Tarlor级数法在实际问题中很少应用，它的作用在于启发我们去探索更好的方法。3-1Runge-Kutta方法的基本思想根据微分中值定理，存在10，使得)()()(1hxyhxyxynnn由方程（1-1）知，))(,()(hxyhxfhxynnn∴有hxyxynn)()(1))(,(hxyhxfnn（3-2）令))(,(*hxyhxfKnn称为区间],[1nnxx上的平均斜率。由此可见，只要对这个平均斜率提供一种算法，由（3-2）式便可以得到一种计算公式。若取0，即取),(nnyx点的斜率作为整个区间],[1nnxx上的平均斜率，这时),(*nnyxfK，则得Euler公式),(1nnnnyxhfyy,2,1,0n若取1，即取),(11nnyx点的斜率作为整个区间],[1nnxx上的平均斜率，这时),(11*nnyxfK，则得后退的Euler公式),(111nnnnyxhfyy,2,1,0n若取),(1nnyxfK，)),(,(12nnnnyxhfyxfK，令)(2121*KKK即取),(nnyx和),(11nnyx两点斜率的算术平均作为整个区间],[1nnxx上的平均斜率，则得下列公式，),(),()(21121211hKyxfKyxfKKKhyynnnnnn（3-3）这显然是改进的Euler公式（即梯形公式）。我们知道，改进的Euler公式是具有2阶精度，而Euler公式和后退的Euler公式只有1阶精度。可见用两个点斜率的算术平均作为平均斜率比只取一个点的斜率作为平均斜率精度要高些。由此得，如果设法在],[1nnxx内多预测几个点的斜率值，然后将它们加权平均作为平均斜率*K，则有可能构造出具有更高精度的计算公式。这就是Runge-Kutta方法的基本思想。3-2二阶Runge-Kutta方法公式（3-3）是一个特殊的二阶Runge-Kutta方法。它用了nx，1nx这两个点上的斜率值。我们现在在],[1nnxx内取nx，pnx10p将这两个点上的斜率值1K，2K的线性组合作为平均斜率*K。显然，),(1nnyxfK，而),(2pnpnyxfK，对于其中的pny我们用Euler公式来计算得，),(nnnpnyxphfyy（3-4）为此有，),(),()(12122111phKyxfKyxfKKKhyynpnnnnn（3-5）其中含有三个待定参数p,,21。确定这三个待定参数的原则是使得公式（3-5）具有二阶精度。为此将2K在),(nnyx点作二元函数Taylor展开，得nyxnnfffphyxfK)(),(2将1K，2K代入（3-5）的第一式得)(22111Kkhyynn])([221nyxnnnfffphffhynyxnnfffphhfy)()(2221而二阶Taylor级数法公式为nnnnyhyhyy!221比较二者系数得到，三个待定参数p,,21应满足211221p（3-6）满足（3-6）形如（3-5）的所有公式均称为二阶Runge-Kutta公式。特别当2121，p=1时就是改进的Euler公式。3-3三阶Runge-Kutta方法为提高精度，我们现在在],[1nnxx内取nx，pnx，qnx10qp三个点，用这三个点上的斜率值1K，2K，3K的线性组合作为平均斜率*K。显然，),(1nnyxfK，),(12phKyxfKnpn，而3K的预测值我们用),(3qnqnyxfK))(,(2