重庆工商大学数学建模算法讲义第15章 常微分方程的解法

-179-第十五章常微分方程的解法建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以肯定得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如22xydxdy+=，于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。§1常微分方程的离散化下面主要讨论一阶常微分方程的初值问题，其一般形式是⎪⎩⎪⎨⎧=≤≤=0)(),(yaybxayxfdxdy(1)在下面的讨论中，我们总假定函数),(yxf连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L，使得|||),(),(|yyLyxfyxf−≤−这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。所谓数值解法，就是求问题（1）的解)(xy在若干点bxxxxaN==L210处的近似值),,2,1(NnynL=的方法，),,2,1(NnynL=称为问题（1）的数值解，nnnxxh−=+1称为由nx到1+nx的步长。今后如无特别说明，我们总取步长为常量h。建立数值解法，首先要将微分方程离散化，一般采用以下几种方法：（i）用差商近似导数若用向前差商hxyxynn)()(1−+代替)('nxy代入（1）中的微分方程，则得),1,0())(,()()(1L=≈−+nxyxfhxyxynnnn化简得))(,()()(1nnnnxyxhfxyxy+≈+如果用)(nxy的近似值ny代入上式右端，所得结果作为)(1+nxy的近似值，记为1+ny，则有),1,0(),(1L=+=+nyxhfyynnnn（2）这样，问题（1）的近似解可通过求解下述问题⎩⎨⎧==+=+)(),1,0(),(01ayynyxhfyynnnnL（3）得到，按式（3）由初值0y可逐次算出L,,21yy。式（3）是个离散化的问题，称为差分方程初值问题。-180-需要说明的是，用不同的差商近似导数，将得到不同的计算公式。（ii）用数值积分方法将问题（1）的解表成积分形式，用数值积分方法离散化。例如，对微分方程两端积分，得∫+==−+1),1,0())(,()()(1nnxxnnndxxyxfxyxyL(4)右边的积分用矩形公式或梯形公式计算。（iii）Taylor多项式近似将函数)(xy在nx处展开，取一次Taylor多项式近似，则得))(,()()(')()(1nnnnnnxyxhfxyxhyxyxy+=+≈+再将)(nxy的近似值ny代入上式右端，所得结果作为)(1+nxy的近似值1+ny，得到离散化的计算公式),(1nnnnyxhfyy+=+以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。其中的Taylor展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。§2欧拉（Euler）方法2.1Euler方法Euler方法就是用差分方程初值问题（3）的解来近似微分方程初值问题（1）的解，即由公式（3）依次算出)(nxy的近似值),2,1(L=nyn。这组公式求问题（1）的数值解称为向前Euler公式。如果在微分方程离散化时，用向后差商代替导数，即hxyxyxynnn)()()('11−≈++，则得计算公式⎩⎨⎧==+=+++)(),1,0(),(0111ayynyxhfyynnnnL（5）用这组公式求问题（1）的数值解称为向后Euler公式。向后Euler法与Euler法形式上相似，但实际计算时却复杂得多。向前Euler公式是显式的，可直接求解。向后Euler公式的右端含有1+ny，因此是隐式公式，一般要用迭代法求解，迭代公式通常为⎪⎩⎪⎨⎧=+=+=+++++),2,1,0(),(),()(11)1(1)0(1Lkyxhfyyyxhfyyknnnknnnnn（6）2.2Euler方法的误差估计对于向前Euler公式（3）我们看到，当L,2,1=n时公式右端的ny都是近似的，所以用它计算的1+ny会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的所谓局部截断误差。假定用（3）式时右端的ny没有误差，即)(nnxyy=，那么由此算出))(,()(1nnnnxyxhfxyy+=+（7）-181-局部截断误差指的是，按（7）式计算由nx到1+nx这一步的计算值1+ny与精确值)(1+nxy之差11)(++−nnyxy。为了估计它，由Taylor展开得到的精确值)(1+nxy是)()(''2)(')()(321hOxyhxhyxyxynnnn+++=+（8）（7）、（8）两式相减（注意到),('yxfy=）得)()()(''2)(23211hOhOxyhyxynnn≈+=−++（9）即局部截断误差是2h阶的，而数值算法的精度定义为：若一种算法的局部截断误差为)(1+phO，则称该算法具有p阶精度。显然p越大，方法的精度越高。式（9）说明，向前Euler方法是一阶方法，因此它的精度不高。§3改进的Euler方法3.1梯形公式利用数值积分方法将微分方程离散化时，若用梯形公式计算式(4)中之右端积分，即))](,())(,([2))(,(111+++≈∫+nnnnxxxyxfxyxfhdxxyxfnn并用1,+nnyy代替)(),(1+nnxyxy，则得计算公式)],(),([2111+++++=nnnnnnyxfyxfhyy这就是求解初值问题（1）的梯形公式。直观上容易看出，用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为⎪⎪⎩⎪⎪⎨⎧=++=+=+++++),2,1,0()],(),([2),()(11)1(1)0(1Lkyxfyxfhyyyxhfyyknnnnnknnnnn(10)由于函数),(yxf关于y满足Lipschitz条件，容易看出||2||)1(1)(1)(1)1(1−+++++−≤−knknknknyyhLyy其中L为Lipschitz常数。因此，当120hL时，迭代收敛。但这样做计算量较大。如果实际计算时精度要求不太高，用公式（10）求解时，每步可以只迭代一次，由此导出一种新的方法—改进Euler法。3.2改进Euler法按式（5）计算问题（1）的数值解时，如果每步只迭代一次，相当于将Euler公式与梯形公式结合使用：先用Euler公式求1+ny的一个初步近似值1+ny，称为预测值，然后用梯形公式校正求得近似值1+ny，即-182-⎪⎩⎪⎨⎧++=+=++++校正预测]),(),([2),(1111nnnnnnnnnnyxfyxfhyyyxhfyy（11）式（11）称为由Euler公式和梯形公式得到的预测—校正系统，也叫改进Euler法。为便于编制程序上机，式（11）常改写成⎪⎪⎩⎪⎪⎨⎧+=++=+=+)(21),(),(1qpnpnnqnnnpyyyyhxhfyyyxhfyy（12）改进Euler法是二阶方法。§4龙格—库塔（Runge—Kutta）方法回到Euler方法的基本思想—用差商代替导数—上来。实际上，按照微分中值定理应有10),(')()(1+=−+θθhxyhxyxynnn注意到方程),('yxfy=就有))(,()()(1hxyhxhfxyxynnnnθθ+++=+（13）不妨记))(,(hxyhxfKnnθθ++=，称为区间],[1+nnxx上的平均斜率。可见给出一种斜率K，（13）式就对应地导出一种算法。向前Euler公式简单地取),(nnyxf为K，精度自然很低。改进的Euler公式可理解为K取),(nnyxf，),(11++nnyxf的平均值，其中),(1nnnnyxhfyy+=+，这种处理提高了精度。如上分析启示我们，在区间],[1+nnxx内多取几个点，将它们的斜率加权平均作为K，就有可能构造出精度更高的计算公式。这就是龙格—库塔方法的基本思想。4.1首先不妨在区间],[1+nnxx内仍取2个点，仿照（13）式用以下形式试一下⎪⎩⎪⎨⎧++==++=+1,0),,(),()(12122111βαβαλλhkyhxfkyxfkkkhyynnnnnn（14）其中βαλλ,,,21为待定系数，看看如何确定它们使（14）式的精度尽量高。为此我们分析局部截断误差11)(++−nnyxy，因为)(nnxyy=，所以（14）可以化为-183-⎪⎪⎪⎩⎪⎪⎪⎨⎧+++=++===++=+)())(,())(,())(,())(,()('))(,()()(2112122111hOxyxfhkxyxhfxyxfhkxyhxfkxyxyxfkkkhxyynnynnxnnnnnnnnnβαβαλλ（15）其中2k在点))(,(nnxyx作了Taylor展开。（15）式又可表为)()()(')()(322211hOfffhxhyxyyyxnnn+++++=+αβαλλλ注意到)()(''2)(')()(321hOxyhxhyxyxynnnn+++=+中fy='，yxfffy+=''，可见为使误差)()(311hOyxynn=−++，只须令121=+λλ，212=αλ，1=αβ（16）待定系数满足（16）的（15）式称为2阶龙格—库塔公式。由于（16）式有4个未知数而只有3个方程，所以解不唯一。不难发现，若令2121==λλ，1==βα，即为改进的Euler公式。可以证明，在],[1+nnxx内只取2点的龙格—库塔公式精度昀高为2阶。4.24阶龙格—库塔公式要进一步提高精度，必须取更多的点，如取4点构造如下形式的公式：⎪⎪⎪⎩⎪⎪⎪⎨⎧++++=+++=++==++++=+),(),(),(),()(3625143423122311121443322111hkhkhkyhxfkhkhkyhxfkhkyhxfkyxfkkkkkhyynnnnnnnnnnβββαββαβαλλλλ（17）其中待定系数iiiβαλ,,共13个，经过与推导2阶龙格—库塔公式类似、但更复杂的计算，得到使局部误差)()(511hOyxynn=−++的11个方程。取既满足这些方程、又较简单的一组iiiβαλ,,，可得-184-⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=++=++==+++=+),()2,2()2,2(),()22(6342312143211hkyhxfkhkyhxfkhkyhxfkyxfkkkkkhynnnnnnnnn（18）这就是常用的4阶龙格—库塔方法（简称RK方法）。§5线性多步法以上所介绍的各种数值解法都是单步法，这是因为它们在计算1+ny时，都只用到前一步的值ny，单步法的一般形式是)1,,1,0(),,(1−=+=+NnhyxhyynnnnLϕ（19）其中),,(hyxϕ称为增量函数，例如Euler方法的增量函数为),(yxf，改进Euler法的增量函数为))],(,(),([21),,(yxhfyhxfyxfhyx+++=ϕ如何通过较多地利用前面的已知信息，如rnnnyyy−−,,,1L，来构造高精度的算法计算1+ny，这就是多步法的基本思想。经常使用的是线性多步法。让我们试验一下1=r，即利用1,−nnyy计算1+ny的情况。从用数值积分方法离散化方程的（4）式∫+=−+1))(,()()(1nnxxnndxxyxfxyxy出发，记nnnfyxf=),(，111),(−−−=nnnfyxf，式中被积函数))(,(xyxf用二节点),(11−−nnfx，),(nnfx的插值公式得到（因)nxx≥，所以是外插。])()[(1))(,(111111−−−−−−−−−=−−+−−=nnnnnnnnnnnnfxxfxxhxxxxfxxxxfxyxf（20）此式在区间],[1+nnxx上积分可得1223))(,(1−−=∫+nnxxfhfhdxxyxfnn于是得到)3(211−+−+=nnnnffhyy（21）注意到插值公式（20）的误差项含因子))((1nnxxxx−−