第八章常微分方程的数值解法详解

第八章常微分方程的数值解法实际问题中，经常需要求解常微分方程。只有一些特殊类型的常微分方程有解析解，实际中的大多数常微分方程无法求得解析解。在多数的工程应用中，往往只需要获得若干点上的函数值，即求数值解。包含自变量、未知函数和未知函数的导数或微分的方程称为微分方程。在求解微分方程时，必须附加某种定解条件。微分方程和定解条件一起组成定解问题。定解条件通常有两种，一种是给出积分曲线在初始时刻的性态，这类条件称为初始条件，相应的定解问题就是初值问题；另一种是给出了积分曲线首末两端的性态，这类条件称为边界条件，相应的定解问题称为边值问题。第八章常微分方程的数值解法未知函数为一元函数的微分方程称为常微分方程，未知函数为多元函数，从而有多元函数偏导数的方程称为偏微分方程。微分方程中各阶导数的最高阶数称为微分方程的阶。本章讨论一阶常微分方程的初值问题。由常微分方程基本理论知道：当在矩形区域上连续，并且关于变量满足李普希兹（Lipschitz）条件，即对任意的，都有其中为常数，则一阶常微分方程初值问题的解存在且唯一。),(yxfdycbxayxD,|),(yDyxDyx),(,),(212121),(),(yyLyxfyxfL第八章常微分方程的数值解法求解常微分方程有各种各样的解析方法，但是解析方法只能求解一些特殊类型的方程，求解实际应用中得到的微分方程主要利用数值解法。对于一阶常微分方程初值问题求其数值解，就是计算出解函数在离散节点处的近似值。)(,)(),(000xxyxyyxfdxdy)(xybxxxan10nyyy,,10第八章常微分方程的数值解法相邻两个节点的距离称为步长，假设节点是等距的，即为常数，这时对应的函数值为这样把一个连续问题转化为一个离散型问题的过程称为离散化过程。离散化过程是把连续的微分方程初值问题转化为一个离散的差分方程问题，然后将求得的差分方程初值问题的解，作为微分方程的解在处的值的近似值。这种离散变量法又称为差分法，对应的离散方程称为差分方程或差分公式。iiixxh1ihnihixxi,,2,1,0,0niihxyxyi,,1,0),()(0iyxyixxixy§1欧拉（Euler）方法欧拉公式1.数值差分法微分方程数值求解的第一步就是设法消除导数项，称为离散化。由于差分是微分的近似运算，实现离散化的基本途径是用差商代替导数。对于初值问题用向前差商代替上式中的导数项，即)()(),,(000xxyxyyxfyhxyxyxyiii)()()(1§1欧拉（Euler）方法得到假设，则得到显式欧拉公式若用向后差商近似导数得到假设，则得隐式欧拉公式))(,()()(1iiiixyxhfxyxynixyyyxhfyyiiii,,2,1,0,)(),(001)(iixyyhxyxyxyiii)()()(11))(,()()(111iiiixyxhfxyxy)(iixyy1,,2,1,0),,(111niyxhfyyiiii§1欧拉（Euler）方法类似地，利用中心差商代替导数可得两步欧拉公式计算时，需要用到前两步信息。hxyxyxyiii2111,,2,1),,(211niyxhfyyiiii1iy1,iiyy§1欧拉（Euler）方法xy0x1x2x3x0P1P2P3PnPQnxO2.几何方法§1欧拉（Euler）方法如图所示，在几何上，的精确解是一条过点的积分曲线，即图中曲线。方程的显式欧拉公式近似解对应着折线上的一系列点的纵坐标。首先过点作斜率为的直线取，得；再过点作斜率为的直线，令得。依此类推，得到折线上的一系列点，对应的纵坐标就是显式欧拉公式的数值解。),(yxfy)(xyy),(000yxPQP0nPP0),(000yxP),(00yxf))(,(0000xxyxfyy1xx),(0001yxhfyy),(111yxP),(11yxf))(,(1111xxyxfyy2xx),(1112yxhfyynyyy,,10§1欧拉（Euler）方法一般地，已求得点，过此点以为斜率做直线当时，得到取。这样就可以从初始条件开始，依次计算出对应的数值解。欧拉公式的几何意义就是用一条初始点重合的折线来近似表示曲线。iiyx,iiyxf,iiiixxyxfyy,1ixxiiiiiixxyxfyy11,11iiyxy00yxynxxx,,,21nyyy,,,21xyy§1欧拉（Euler）方法3.数值积分法把方程的两端从到积分，得到近似地计算出上式中的积分项，就能够获得的近似值；而选用不同的计算方法计算此积分项，就会得到不同的差分公式。用左矩形公式计算积分项得到就是欧拉公式。yxfy,ix1ix111,,1iiiiiixxiixxxxdxxyxfxyxydxxyxfdxy1ixyiixxxyxhfdxxyxfii,,1iiiixyxhfxyxy,1§1欧拉（Euler）方法4.泰勒展开法对在处做二阶泰勒展开做线性近似，得到用近似值代替，则上式写为1ixyix121,!2iiiiiiiixxyhxyhxyhxyxyiiiiiixyxhfxyxyhxyxy,1iyixyiiiiyxhfyy,1§1欧拉（Euler）方法欧拉预估-校正方法对方程两边从到积分得利用梯形公式计算积分项有假设，则得梯形公式),(yxfyix1ix1),()()(1iixxiidxyxfxyxy))(,())(,(2)()(111iiiiiixyxfxyxfhxyxy)(iixyy1,,1,0,),(),(2111niyxfyxfhyyiiiiii§1欧拉（Euler）方法特点：隐式，不便计算；处理方法：用欧拉公式先对进行预报，记为，称为预报值，再将预报值代入梯形公式右端，得到校正的，称为欧拉预估-校正公式或改进的欧拉公式等价形式：1iy1iy1iy)],(),([2),(1111iiiiiiiiiiyxfyxfhyyyxhfyy校正预估1,,1,0,)(21),(),(11niyyyyxhfyyyxhfyycpipiiciiip§1欧拉（Euler）方法例1利用欧拉公式和预估-校正公式求初值问题在区间上的数值解（取），并与精确解进行比较。解：依题意有欧拉公式：1)0(,2yyxyy]1,0[1.0h12xy1.0,1)0(,2),(hyyxyyxf)9,,1,0()2(1.0),(1iyxyyyxhfyyiiiiiiii§1欧拉（Euler）方法预估-校正公式：)9,,1,0()22(05.0)],(),([2)2(1.0),(1111111iyxyyxyyyxfyxfhyyyxyyyxhfyyiiiiiiiiiiiiiiiiiiiii§1欧拉（Euler）方法结果列表如下：iixiyiy)(ixy欧拉公式预估-校正公式精确值节点号0011110.11.11.0959091.09544520.21.1918181.1840971.18321630.31.2774381.2662011.26491140.41.3582131.3433601.34164150.51.4351331.4164021.41421460.61.5089661.4859561.48324070.71.5803381.5525141.54919380.81.6497831.6164751.61245290.91.7177791.6781661.673320101.01.7847711.7378671.732051§1欧拉（Euler）方法局部截断误差：假定为准确值，即，用某种数值方法计算的误差称为该数值方法计算的局部截断误差。阶方法：若某一数值方法的局部截断误差为，为正整数，则称该方法为阶方法，或说该方法具有阶精度。pppp)(1pihORiy)(iixyy1iy11)(iiiyxyR1iy§1欧拉（Euler）方法已知则欧拉方法的误差将在点泰勒展开有局部截断误差：误差为，一阶方法。))(,()(xyxfxy)())(,())(,()(xyyxyxfxxyxfxy)(1ixyix)(!2)()()(!2)()()(221yhxyhxyxyhxyhxyxyiiiiii)(2)(2)(22211hOyhxyhyxyiii)(2hO§1欧拉（Euler）方法欧拉预估-校正公式将在点泰勒展开有则)],(),([2111iiiiiiyxfyxfhyy),(11iiyxfix)()()()()())(,())(,())(,(),(2211hOxyhxyhOxyxyxfhxyxfhxyxfyxfiiiiiyiixiiii)()(2)()()]()()([2)(23221hOxyhxyhxyhOxyhxyhxyhyyiiiiiiii§1欧拉（Euler）方法将在点处做泰勒展开有局部截断误差：误差为，是二阶方法。小结：数值方法误差分析常用泰勒展开，展开的阶数视截断误差的阶数而定。)(1ixyix)()(!2)()()(321hOxyhxyhxyxyiiii)()(311hOyxyii)(3hO§1欧拉（Euler）方法例：证明下面求解公式是一个3阶方法。)32,32()3,3(),()3(423121311khyhxfkkhyhxfkyxfkkkhyyiiiiiiii§1欧拉（Euler）方法证明：))(,()(xyxfxy)())(,())(,()(xyyxyxfxxyxfxyyxyxfxyyxyxfxyyxxyxfxyxxyxfxyyxyxfyxxyxfxyyxyxfxyxyyxxyxfxxyxfxy))(,()())(,()())(,()(2))(,()())(,())(,()())(,()()())(,())(,()(222222222222§1欧拉（Euler）方法由已知，可以分别记为局部截断误差为321,,kkk))(,(1iixyxfk)3)(,3(12khxyhxfkii)32)(,32(23khxyhxfkii)3(4)()()(311111kkhxyxyyxyRiiiii§1欧拉（Euler）方法又)())(,(1iiixyxyxfk)())(,()()(18)(3)()())(,()())(,()(2))(,(18))(,()())(,(3)()())(,(9))(,(332))(,(9!21))(,(3))(,(3))(,()3)(,3(323222222232221221222112hOyxyxfxyxyhxyhxyhOyxyxfxyyxxyxfxyxxyxfhyxyxfxyxxyxfhxyhOyxyxfkhyxxyxfkhhxxyxfhyxyxfkhxxyxfhxyxfkhxy