[理学]第九章 常微分方程的数值解法.doc

第九章 常微分方程的数值解法 对于一阶常微分方程初值问题:例如微分方程:xy-2y=4x转变为y=2/x+4,这儿f(x,y)=2y/x+4,加上定解条件(初始条件):y(1)=-3,得一阶常微分方程的初始问题,微分方程的定解问题是求一个函数y=y(x)使得该函数满足微分方程并且符合初值条件。例如对于函数y(x)=x2-4x满足以上条件,因而是该初始问题的解。 但是,只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表示，而在实际工程技术、生产、科研上会出现大量的微分方程问题，需要得到其解，因此只能依赖于数值方法去获得微分方程的数值解。ab x0 x1 x2 . xn-1 xn微分方程的数值解：设方程问题的解y(x)的存在区间是a,b，初始点x0=a，将a,b内得一系列节点x0 , x1 ,.,xn。a= x0< x1<< xn =b,其中hk=xk+1-xk , 如是等距节点h=(b-a)/n , 对一切k，h= xk+1-xk ,k=0,1,2,.n-1。h称为步长，一般采用等距步长。y(x)的解析表达式不容易得到或根本无法得到，我们用数值方法求得y(x)在每个节点xk的值y(xk)的近似值，用yk表示，即yy(xk)，这样y0 , y1 ,.,yn称为微分方程的数值解。如果用数值方法，我们求的是f(xk)的近似值yk，如果要用差值方法或拟全方法再求y(xk)的近似函数。上例中微分方程的解y(x)=x2-4x设解的存在区间为1，3，n=10，h=(3-1)/10=0.2Kxky(xk)01-311.2-3.3621.4-3.6431.6-3.8441.8-3.9652.0-462.2-3.9672.4-3.8482.6-3.6492.8-3.36103.0-3§1 欧拉公式和改进欧拉方法一、欧拉公式： 1、方法构造的思想：微分方程初值问题 ，微分方程在(x0,y0)处成立，再用前差代替一阶导数，则。此时x0, y0均已知y(x0)= y0,则由此式可近似求出y(x1)的近似值y1，同样可利用x1处的微分方程可得，一般的利用在xn处的微分方程可得： 此式称为微分方程的差分格式，对于这个差分格式我们称为欧拉公式(或称欧拉格式Euler Scheme)。2、几何意义：对于微分方程y=2(x+1),其通解是y=(x+1)2+c,是一个曲线族(本问题是抛物线族)，如果加上定解条件y(0)=2，则是一个确定的解，y=(x+1)2+1,由y(0)=2，过该曲线上一点P0(x0,y0)作曲线的切线，其斜率，切线为YXP0P1Q1y1y(x1)h例：h=0.5,x1=x0+h切线与x=x1的交点为Q1 (x1,y1), x1-x0=h y1=y0+hf(x0，y0)这就是由欧拉方法求出的y1，而x=x1与曲线y=y(x)的交点是P1(x1, y(x1)，y(x1)与y1是有误差的。在本问题中h=0.5，则y(x1)=1.52，而y1=y0+hf(x0，y0)=3，误差为0.25。例：y=2x，y(0)=1，h=0.2，其解为y=x2+1，y(x1)=0.22+1=1.04y1=y0+hf(x0，y0)=1+h×0=1Q1点实际在曲线族的另一个抛物线y=x2+0.96上，过Q1点作切线x=x2交于Q2(x2, y2)。y(x2)=0.42+1=1.16 y2=y1+hf(x1，y1)=1+0.2×2×0.2=1.08Q2实际上在y=x2+0.92上，再通过Q2作曲线y=x2+0.92的切线，交x=x3于Q3(x3, y3)继续得到一系列的Q1，Q2， . ，Qn，得一折线P0Q1Q2. Qn。故欧拉法又称欧拉折线法。P4YP3P2P1Q4Q3Q1Q2h h h h x1 x2 x3 x4X例1：解：h=0.2 , xi=1+ih y1= y0+hf(x0，y0)=-1+0.2× y2= y1+hf(x1，y1)=-1+0.2× y3= y2+hf(x2，y2)= -0.9333+0.2×y4= y3+hf(x3，y3)=0.8+0.2×y5= y4+hf(x4，y4)=0.6+0.2×y6= y5+hf(x5，y5)=0.3333+0.2×精确解为：y=x2-2xxky(xk)ykek1.2-0.96-10.041.4-0.84-0.93330.09331.6-0.64-0.80.161.8-0.36-0.60.242.00-0.33330.33332.20.4400.44 可以看出误差随着计算在积累。格式：。特点：(1)、单步方法；(2)、显式格式；(3)、局部截断误差因而是一阶精度 。局部截断误差：当是精确解下，由按照欧拉方法计算出来的的误差称为局部截断误差。即，则称为局部截断误差。即是局部截断误差。换言之，局部截断误差是差分格式中均换成精确解时，所截断部分，即局部截断误差。例如：由，则欧拉格式是将项截断得：因此 局部截断误差是。自我练习例1：步长0.21.00000.96078940.03920.40.92000.85214370.06790.60.77280.69767630.07510.80.5873280.52729240.06001.00.3993830.36787940.03151.20.23962980.23692770.0027一、 改进欧拉方法：1、 格式：2、 特点：(1)、是单步方法；(2)、局部截断误差是因而是二阶精度，截断余项得改进欧拉格式，所以局部截断误差是；(3)、是隐式格式，无法从格式中直接求出必须要解方程。 3、用预测校正方法来求隐式格式中的。 例： 解：(1)、(2)、(3)、(4)、(5)、(6)、0.20.960.96078940.40.8509440.85214370.60.69709330.69767630.80.52867550.52729241.00.37218870.36787941.20.24415570.2369277 一般如果局部截断误差是,我们就说该方法具有p阶精度。,两者是有误差的,因而,因而书上的推导是错误的,其原因是局部截断误差的定义是错误的。可以证明:如果局部截断误差达到,满足某一条件,则整体误差因而称其精度达到一阶。4.数值例子:例1: 求的近似值。解:这儿,由欧拉公式得:,一直计算下去。精确解为:。从表中看出误差在逐步增加、积累,因是由不准确的来计算的,来看局部截断误差。局部截断误差,而误差是,显然很两者有明显的差异,前者是量级,而后者是量级。二.改进Euler方法:1. 格式的形成和构造:在微分方程初值问题,对其从到进行定积分得:将右端的定积分用梯形公式来进行近似计算。用和来分别代替和得差分格式:这就是改进欧拉方法。2. 显式格式和隐式格式:在欧拉式中每一步计算已知,直接用格式可以计算出,此类格式称为显式格式。而在改进欧拉方法中在每一步计算中是未知,待求的,未知量在中这是一个方程,如是非线形或超越函数,此方程是无法直接解出来(要依靠迭代法才能解出)。这类格式称为隐式格式。例: 用改进欧拉方法求解。解得:由于,是线形函数可以从隐式格式中解出。问题的精确解是欧拉方法误差可见改进欧拉方法误差比欧拉方法要好的多。yxyx欧拉方法相当于在中用一个底作高的矩阵来代替曲边梯形的面积。3.用预测校正方法求隐式格式的解:有时要从中要解出是相当困难的,因此我们用预测校正的方法计算,在已得到的情况下, 预测值: 校正值:此式相当于对隐式格式求时采用迭代的方法,用欧拉格式得到的作为初始值迭代公式迭代一次而已,此公式代入后得:如改写成平均的形式为:例2、。解：改进欧拉格式为 解yk+1出来比较困难，遇到的是一个二次方程， 关于yk+1是一元二次方程，求解比较难，我们用预测校正方法来求。 (2)求y2(3)求y3xkyky(xk)ek0.11.09596911.09544510.0004640.21.18409661.18321600.0008810.31.26620141.26491110.0012900.41.34336021.34164080.0017190.51.41640191.41421360.0021880.61.47595561.48323970.0027160.71.55251411.54919330.0033210.81.61647481.61245150.0040230.91.67816641.67332010.0047461.01.73786741.73205080.0058174、改进欧拉方法的截断误差假定y(xk)是精确的,由此一步得yk+1。两式相减得：改进欧拉方法的局部截断误差比欧拉方法高出一次是O(h3),也称之为具有二阶精度的。§14.2 龙格-库塔法一、龙格-库塔法的思想 考虑微分方程的初值问题:如果用均差代替导数用微分中值定理,存在0<<1 ,从而得：称为y(x)在区间xk, xk+1上的平均斜率,记作K,但中值是存在而未知的,对平均斜率K,提供一种近似的计算方法,就得到的一种近似公式，或称为微分方程的一种计算格式。 前面用f(xk,yk)作为平均斜率K的近似值就得到欧拉格式：。而平均斜率用就得到比欧拉格式高一阶精度的格式。而k1=f(xk, yk),k2=f(xk+1, yk+1)= f(xk+1, yk+hk1)用来作为平均斜率的近似值，就得到了改进欧拉格式的预测校正方法。 由此启发，能否在二维平面中xxk, xk+1, yyk, yk+1上多找一些f(x,y)在这些点上的函数值的平均数；来作为平均的近似值，由于自由度的增加，使得的p能够提高，从而达到提高精度的目标，这就是龙格库塔法的思想。Y(xk+1 , yk+1) yk+1(xk+1 ,)hf(xk) yk(xk , yk) hxk+1xk 在以上二维平面的区域上取N个点，得到公式，其中Km是二元函数f(x,y)在这N个点上的值。X K1=f(xk , yk),.其中是待定常数，适当选取可以提高精度，此方法在已计yk算好后，依次计算，最后计算yk+1均是直接带入公式计算，不必解方程，因而称为显式龙格-库塔法，当然也有隐式龙格-库塔法。二、龙格-库塔法1、 计算公式：用泰勒展式来得到这些参数应满足的方程，将f(x,y)在(xk , yk)上展开K1=f(xk , yk)=y(xk),而因y=f 以上均代入得比较两边的式子得:四个未知量,三个方程,有无穷多组解，但局部截断误差均为是是二阶精度的。 例:取,则,就是改进欧拉公式的预测-校正方法。如取，则，则,得二阶龙格库塔法为：此格式称为二阶中点格式。三、三阶龙格库塔方法:格式:参数有,共有八个参数,将f和均在点展开。利用,按展开式代入(高于的不计)得:以上展开式代入中比较各次得:(1)(2)(3)(4)(5)(6)(7)(8)八个方程真正独立的是六个方程。八个未知量,六个方程有无穷多组解,但其截断误差均为均是三阶精度。例如1:得:此方法称为Kutta方法。例如2:得:此方法称为Heun法。例如:求解微方分程初值问题：选,在区间上用三阶Kutta方法。其精确解为:0.21.9240961.9230770.0010190.41.7257011.7241380.0015630.61.4720031.47058820.0014150.81.2203811.2195120.0008691.01.0003711.000000.000371四、四阶龙格库塔方法: 显式四阶龙格库塔方法的一般形式是:同样用泰勒展开的方法,将f展开到,将展开到,四个,三个,六个共十三个未知待定参数,共有十一个独立的方程的二个自由度。十三个未知量，十一个方程有无穷多组解。这样就得到格式称为古典龙格库塔格式，书上称为标准龙格库塔方法，其格式为：得到的格式称为库塔法，格式为：得到的方法称为吉尔方法。总之，这三种方法的共同点是截断误差为阶数为四阶精度，均称为显式四阶龙格库塔方法，在从来计算出均要计算四个f的值一般一阶常微分方程初值问题均用四阶龙格库塔方法来计算，其精度均满足了实际问题的精度要求。数值例子：考虑微分方程初值问题：。解：其精确解为：，用古典龙格库塔方法，。(1)、求，此时(2)、求，此时以下计算用表格列出：误差1.11.200271.200271.21.401821.401821.31.605241.605241.41.810791.710791.52.018562.018561.62.228552.228551.72.440702.440701.82.654972.654971.92.871272.871272.03.089533.08953例：初值问题用四阶古典RungeKutta方法，。(1)、求，(2)、求，0.21.18322931.18321600.41.34168031.34164080.61.48328381.48323970.81.61251721.61245151.01.73214631.7320508§14.3 亚当斯方法(Adams)一、 单步方法和多步方法：前面讲的方法：欧拉方法、改进欧拉方法、龙格库塔方法均是单步方法，即在每一步要计算时，只要前面一个值已知的条件下秒可以计算出了， 特点：(1)、可以自成系统进行直接计算，因为初始条件只有一个已知，由可以计算，不必借助于其它方法，这种称为单步方法是自开始的。(2)、如果格式简单如欧拉方法，则只有一阶精度，如果提高精度，则计算很复杂，如RungeKutta方法。(3)、公式的构造推导也很复杂。多步方法：利用前面已知计算出来的，这前面计算好的个值来计算，这样自由度的增加来提高格式的精度，这样的方法称为多步方法，利用k个值计算，称为k步方法。多步方法的特点：(1)、因初始条件只有一个，运用多步方法设法开始，要借助高阶的单步方法来开始，例如，已知用单步的四阶RungeKutta方法计算，再计算，再由计算，用单步方法有后运用四阶的四步方法，由计算；由计算；由计算；一直下去均匀 可以用多步方法了，而且始终达到四阶精度。(2)、多步方法比较简单，只要在这四个点的函数值的线性组合，而且每步中后三个函数值下一步还可使用。二、 显式Adams方法：1、 方法构造的思想：Adams方法是多步方法中的某一类，而不就是多步方法，考虑微分方程初值问题，将微分方程在上积分，下面我们来推导四步显式Adams方法，即若已知来计算，简记，用的拉格朗日插值多项式代替f，截断部分，用等距步长，上面积分很简单，得到的方法就是显式Adams方法。以上可以看到该方法的局部截断误差是因而是四阶精度的。例如：解：取，首先用四阶RungeKutta方法来起步，计算出，下面不必用RungeKutta方法，而开始用四阶Adams方法。(1)、求(2)、求只要补算(3)、求只要补算现列表看用Adams方法求出的误差，精解为0.81.61142311.61245151.01.72984031.73205081.21.84066161.8439089三、 隐式Adams方法：与隐式不同我们用作为插值结点，关键不同在于也是插值结点，必带来从而导致是隐式格式。用插值多项式来代替积分中的得：截掉得近似公式：得：，从而得四阶隐式Adams方法。因，而是未知的，故这是隐式格式。隐式格式的解法用预测校正法：用显式格式作为预测值，再用隐式格式来校正。预测值：校正值：例：精确解为，取步长，用四阶多步方法，预测校正方法。首先用四阶龙格库塔方法求，(1)、求用显式四步格式得预测值：用隐式格式来校正：精确解为预测值的误差为：校正值的误差为：，精度大大提高。(2)、求只要补算预测值：校正值：精确解为：预测值的误差为：校正值的误差为：(3)、求只要补算预测值：校正值：精确解为：预测值的误差为：校正值的误差为：例：，用。解：先用RungeKutta方法求起步。(1)、求(2)、求(3)、求然后用四阶Adams方法，预测校正方法求。(4)、求预测值：校正值：(5)、求补充：预测值：校正值：(6)、求补充：预测值：校正值：精确解为：，列表如下：0.80.52664330.52729241.00.36682550.36787941.20.23598380.2369277自我练习例1：步长0.21.00000.96078940.03920.40.92000.85214370.06790.60.77280.69767630.07510.80.5873280.52729240.06001.00.3993830.36787940.03151.20.23962980.23692770.0027二、 改进欧拉方法：3、 格式：4、 特点：(1)、是单步方法；(2)、局部截断误差是因而是二阶精度，截断余项得改进欧拉格式，所以局部截断误差是；(3)、是隐式格式，无法从格式中直接求出必须要解方程。 3、用预测校正方法来求隐式格式中的。 例： 解：(1)、(2)、(3)、(4)、(5)、(6)、0.20.960.96078940.40.8509440.85214370.60.69709330.69767630.80.52867550.52729241.00.37218870.36787941.20.24415570.2369277三、RungeKutta法：二阶RungeKutta方法：系数：三阶RungeKutta方法：系数：局部截断误差，具有三阶精度。四阶RungeKutta方法：系数：截断误差，达到四阶精度。而根据RungeKutta方法的系数表为：01001自我练习例 3：取步长，用四阶RungeKutta方法求解实值问题。同样用泰勒展开的方法，将f展开到，将展开到,四个,三个，六个共十三个未知待定参数,共有十一个独立的方程的二个自由度。，十三个未知量，十一个方程有无穷多组解。这样就得到格式称为古典龙格库塔格式，书上称为标准龙格库塔方法，其格式为：得到的格式称为库塔法，格式为：得到的方法称为吉尔方法。总之,这三种方法的共同点是截断误差为阶数为四阶精度,均称为显式四阶龙格库塔方法，在从来计算出均要计算四个f的值一般一阶常微分方程初值问题均用四阶龙格库塔方法来计算，其精度均满足了实际问题的精度要求。数值例子：考虑微分方程初值问题：。解：其精确解为：，用古典龙格库塔方法,。(1)、求，此时(2)、求，此时以下计算用表格列出：误差1.11.200271.200273´10-72.7´10-41.21.401821.401825´10-71.06´10-31.31.605241.605246´10-72.08´10-31.41.810791.710796´10-73.22´10-31.52.018562.018567´10-74.42´10-31.62.228552.228557´10-75.66´10-31.72.440702.440708´10-76.90´10-31.82.654972.654978´10-78.16´10-31.92.871272.871278´10-79.40´10-32.03.089533.089538´10-710.65´10-3 例：初值问题用四阶古典RungeKutta方法，。(1)、求， (2)、求，0.21.18322931.18321601.3´10-50.41.34168031.34164084.0´10-50.61.48328381.48323974.4´10-50.81.61251721.61245156.6´10-51.01.73214631.73205089.6´10-5§14.3 亚当斯方法(Adams)一. 单步方法和多步方法：前面讲的方法：欧拉方法、改进欧拉方法、龙格库塔方法均是单步方法，即在每一步要计算时,只要前面一个值已知的条件下秒可以计算出了。特点：(1)、可以自成系统进行直接计算，因为初始条件只有一个已知，由可以计算，不必借助于其它方法,这种称为单步方法是自开始的。(2)、如果格式简单如欧拉方法，则只有一阶精度，如果提高精度，则计算很复杂.如RungeKutta方法。(3)、公式的构造推导也很复杂。多步方法：利用前面已知计算出来的，这前面计算好的个值来计算，这样自由度的增16加来提高格式的精度，这样的方法称为多步方法，利用k个值计算，称为k步方法。多步方法的特点：(1)、 因初始条件只有一个，运用多步方法设法开始，要借助高阶的单步方法来开始.例如,已知用单步的四阶RungeKutta方法计算，再计算，再由计算，用单步方法有后运用四阶的四步方法，由计算；由计算；由计算；一直下去均匀可以用多步方法了，而且始终达到四阶精度.(2)、多步方法比较简单，只要在这四个点的函数值的线性组合,而且每步中后三个函数值下一步还可使用。二. 显式Adams方法：1. 方法构造的思想：Adams方法是多步方法中的某一类，而不就是多步方法，考虑微分方程初值问题，将微分方程在上积分，下面我们来推导四步显式Adams方法，即若已知来计算，简记用的拉格朗日插值多项式代替f二.显式Adams方法:1.方法构造的思想:Adams方法是多步方法中的某一类,而不就是多步方法,考虑微分方程初值问题,将微分方程在上积分,下面我们来推导四步显式Adams方法,即若已知来计算,简记,用的拉格朗日插值多项式代替f,而截断部分,用等距步长,上面积分很简单,得到的方法就是显式Adams方法。以上可以看到该方法的局部截断误差是因而是四阶精度的。例如:解:取,首先用四阶RungeKutta方法来起步,计算出,下面不必用RungeKutta方法,而开始用四阶Adams方法。(1)求,。(2)求只要补算(3)求只要补算现列表看用Adams方法求出的误差,精解为0.81.61142311.61245151.01.72984031.73205081.21.84066161.8439089三.隐式Adams方法: 与隐式不同我们用作为插值结点,关键不同在于也是插值结点,必带来从而导致是隐式格式。用插值多项式来代替积分中的得:截掉得近似公式:得:，从而得四阶隐式Adams方法。因,而是未知的,故这是隐式格式。隐式格式的解法用预测校正法: 用显式格式作为预测值,再用隐式格式来校正。预测值:校正值:例:精确解为,取步长,用四阶多步方法,预测校正方法。首先用四阶龙格库塔方法求,(1)求,。用显式四步格式得预测值:用隐式格式来校正:精确解为预测值的误差为:校正值的误差为:,精度大大提高。(2)求只要补算预测值:校正值:精确解为:预测值的误差为:校正值的误差为:(3)求只要补算预测值:校正值:精确解为:预测值的误差为:校正值的误差为: 例：，用。解：先用RungeKutta方法求起步。(1)、求(2)、求(3)、求然后用四阶Adams方法，预测校正方法求。(4)、求预测值：校正值：(5)、求补充：预测值:校正值:(6)、求补充：预测值：校正值：精确解为：列表如下：0.80.52664330.52729241.00.36682550.36787941.20.23598380.2369277例：解：(1)、 (2)、 (3)、 (4)、(5)、 (6)、 0.20.960.96078940.40.8509440.85214370.60.69709330.69767630.80.52867550.52729241.00.37218870.36787941.20.24415570.2369277三、RungeKutta法： 二阶RungeKutta方法：系数： 三阶RungeKutta方法：系数：局部截断误差，具有三阶精度。四阶RungeKutta方法：系数：截断误差，达到四阶精度。而根据RungeKutta方法的系数表为：01001自我练习例 3：取步长，用四阶RungeKutta方法求解实值问题。,精确解: 解：(1)、求，,(2)、求，(3)、求，(4)、求，(5)、求，(6)、求， 试证任意参数t(0<t<1)下列RungeKutta方法的局部截断误差为。证明：要证要证右边等于：该格式的局部截断误差为。