第3章-常微分方程的差分方法

数值计算方法（数值分析）,第三章：常微分方程的差分方法,包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中,自变量的个数只有一个,称为常微分方程.。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数y及其各阶导数都是一次的,则称它是线性的,否则称为非线性的。,引言,在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。譬如,这个一阶微分方程就不能用初等函数及其积分来表达它的解。,再如，方程,的解,虽然有表可查,但对于表上没有给出的值,仍需插值方法来计算,从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程初值问题,在区间axb上的数值解法。,可以证明,如果函数在带形区域R=axb,-y内连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L(它与x,y无关)使,对R内任意两个都成立,则方程的解在a,b上存在且唯一。,数值方法的基本思想对常微分方程初值问题数值解法，就是要算出精确解y(x)在区间a,b上的一系列离散节点处的函数值的近似值。相邻两个节点的间距称为步长，步长可以相等，也可以不等。本章总是假定h为定数，称为定步长，这时节点可表示为数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。,对常微分方程数值解法的基本出发点就是离散化。其数值解法有两个基本特点，它们都采用“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信息计算的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题中的导数进行不同的离散化处理。,对于初值问题的数值解法，首先要解决的问题就是如何对微分方程进行离散化，建立求数值解的递推公式。递推公式通常有两类，一类是计算yi+1时只用到xi+1,xi和yi，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为单步法；其代表是龙格库塔法。另一类是计算yi+1时，除用到xi+1,xi和yi以外，还要用到，即前面k步的值，此类方法称为多步法；其代表是亚当姆斯法。,Euler公式欧拉（Euler）方法是解初值问题的最简单的数值方法。初值问题的解y=y(x)代表通过点的一条称之为微分方程的积分曲线。积分曲线上每一点的切线的斜率等于函数在这点的值。,简单数值方法,Euler法的求解过程是:从初始点P0(即点(x0,y0)出发,作积分曲线y=y(x)在P0点上切线(其斜率为),与x=x1直线,相交于P1点(即点(x1,y1),得到y1作为y(x1)的近似值,如上图所示。过点(x0,y0),以f(x0,y0)为斜率的切线方程为当时,得,这样就获得了P1点的坐标。,同样,过点P1(x1,y1),作积分曲线y=y(x)的切线交直线x=x2于P2点,切线的斜率=直线方程为,当时,得,当时,得,由此获得了P2的坐标。重复以上过程,就可获得一系列的点:P1,P1,Pn。对已求得点以=为斜率作直线,取,从图形上看,就获得了一条近似于曲线y=y(x)的折线。,这样,从x0逐个算出对应的数值解,通常取(常数),则Euler法的计算格式,i=0,1,n(3.2),还可用数值微分、数值积分法和泰勒展开法推导Euler格式。以数值积分为例进行推导。将方程的两端在区间上积分得，,选择不同的计算方法计算上式的积分项,就会得到不同的计算公式。,(3.3),用左矩形方法计算积分项,代入(3.3)式,并用yi近似代替式中y(xi)即可得到向前欧拉（Euler）公式,由于数值积分的矩形方法精度很低，所以欧拉（Euler）公式当然很粗糙。,例1用欧拉法解初值问题,取步长h=0.2,计算过程保留4位小数,解:h=0.2,欧拉迭代格式,当k=0,x1=0.2时，已知x0=0，y0=1，有y(0.2)y1=0.21(401)0.8当k=1,x2=0.4时，已知x1=0.2，y1=0.8，有y(0.4)y2=0.20.8(40.20.8)0.6144当k=2,x3=0.6时，已知x2=0.4，y2=0.6144，有y(0.6)y3=0.20.6144(4-0.40.6144)=0.4613,clear;y=1,x=0,%初始化forn=1:10y=1.1*y-0.2*x/y,x=x+0.1,end,y=1x=0y=1.1000 x=0.1000y=1.1918x=0.2000y=1.2774x=0.3000y=1.3582x=0.4000y=1.4351x=0.5000y=1.5090 x=0.6000y=1.5803x=0.7000y=1.6498x=0.8000y=1.7178x=0.9000y=1.7848x=1.0000,梯形公式为了提高精度,对方程的两端在区间上积分得，改用梯形方法计算其积分项，即,(3.4),代入(3.4)式,并用近似代替式中即可得到梯形公式,(3.5),梯形公式（3.5）比欧拉公式(3.2)的精度高,(3.5),(3.5)式的右端含有未知的yi+1,它是一个关于yi+1的函数方程,这类数值方法称为隐式方法。相反地,欧拉法是关于yi+1的一个直接的计算公式，这类数值方法称为显式方法。,两步欧拉公式对方程的两端在区间上积分得,(3.6),改用中矩形公式计算其积分项，即,代入上式,并用yi近似代替式中y(xi)即可得到两步欧拉公式,(3.7),前面介绍过的数值方法,无论是欧拉方法,还是梯形方法，它们都是单步法,其特点是在计算yi+1时只用到前一步的信息yi;可是公式(3.7)中除了yi外,还用到更前一步的信息yi-1,即调用了前两步的信息,故称其为两步欧拉公式,欧拉法的局部截断误差衡量求解公式好坏的一个主要标准是求解公式的精度,因此引入局部截断误差和阶数的概念。定义1在yi准确的前提下,即时,用数值方法计算yi+1的误差,称为该数值方法计算时yi+1的局部截断误差。对于欧拉公式，假定，则有,而将真解y(x)在xi处按二阶泰勒展开,因此有,定义2数值方法的局部截断误差为,则称这种数值方法的阶数是P。步长(hN结束。,(3)程序实现(改进欧拉法计算常微分方程初值问题),例9.2用改进欧拉法解初值问题,区间为0,1,取步长h=0.1,解:改进欧拉法的具体形式,本题的精确解为,clearx=0,yn=1%初始化forn=1:10yp=yn+0.1*(yn-2*x/yn);%预测x=x+0.1;yc=yn+0.1*(yp-2*x/yp);yn=(yp+yc)/2%校正end,例3对初值问题,证明用梯形公式求得的近似解为,并证明当步长h0时,yn收敛于精确解证明:解初值问题的梯形公式为,整理成显式,反复迭代,得到,由于，有,证毕,作业题！,1、列出下列初值问题的欧拉格式：（1）0=x=0.4,y(0)=1,取h=0.2（2）1=x=1.2y(1)=1,取h=0.12、用梯形方法做1=，反复将步长折半进行计算，直至为止，并以上一次步长的计算结果作为。这种通过步长加倍或折半来处理步长的方法称为变步长法。表面上看，为了选择步长，每一步都要反复判断，增加了计算工作量，但在方程的解y(x)变化剧烈的情况下，总的计算工作量得到减少，结果还是合算的。,作业题！,3、取h=0.1用欧拉方法，取h=0.2用改进欧拉方法，取h=0.4用四阶龙格-库塔方法求解初值问题：,4、求解微分方程初值问题：,（1）取h=0.1，用欧拉方法求解。(5分)（2）取h=0.4，用经典四阶龙格-库塔方法求解。(6分)（3）精确解y（0.4）=1.51277，分别确定上述两种方法计算结果有几位有效数字。(4分),亚当姆斯格式龙格-库塔方法是一类重要算法，但这类算法在每一步都需要先预报几个点上的斜率值，计算量比较大。考虑到计算yi+1之前已得出一系列节点上的斜率值，能否利用这些已知值来减少计算量呢？这就是亚当姆斯（Adams）方法的设计思想。,亚当姆斯方法,设用xi,xi-1两点的斜率值加权平均作为区间上的平均斜率，有计算格式,（3.21）,选取参数，使格式（3.21）具有二阶精度。,亚当姆斯方法,将在xi处Taylor展开,代入计算格式(3.21)化简,并假设,因此有,与y(xi+1)在xi处的Taylor展开式,相比较,需取,才使格式(3.21)具有二阶精度。这样导出的计算格式,称之为二阶亚当姆斯格式。类似地可以导出三阶亚当姆斯格式。,和四阶亚当姆斯格式。,(3.22),这里和下面均记,上述几种亚当姆斯格式都是显式的，算法比较简单，但用节点的斜率值来预报区间上的平均斜率是个外推过程，效果不够理想。为了进一步改善精度，变外推为内插，即增加节点xi+1的斜率值来得出上的平均斜率。譬如考察形如,（3.23）,的隐式格式,设(3.23)右端的Taylor展开有,可见要使格式(3.23)具有二阶精度,需令,这样就可构造二阶隐式亚当姆斯格式,其实是梯形格式。,类似可导出三阶隐式亚当姆斯格式,和四阶隐式亚当姆斯格式,（3.24）,参照改进的欧拉格式的构造方法，以四阶亚当姆斯为例，将显式（3.22）和隐式（3.24）相结合，用显式公式做预报，再用隐式公式做校正，可构成亚当姆斯预报-校正格式,（3.25）,预报：,校正：,亚当姆斯预报校正-格式,这种预报-校正格式是四步法，它在计算yi+1时不但用到前一步的信息，而且要用到再前面三步的信息，因此它不能自行启动。在实际计算时，可借助于某种单步法，譬如四阶龙格库塔法提供开始值。,例5取步长h=0.1,用亚当姆斯预报-校正公式求解初值问题,的数值解。,解:用四阶龙格-库塔公式求出发值,计算得：,表中的,yi和y(xi)分别为预报值、校正值和准确解(),以比较计算结果的精度。,再使用亚当姆斯预报-校正公式(3.25),稳定性在微分方程的数值解法中是一个非常重要的问题。因为微分方程初值问题的数值方法是用差分格式进行计算的，而在差分方程的求解过程中，存在着各种计算误差，这些计算误差如舍入误差等引起的扰动，在传播过程中，可能会大量积累，对计算结果的准确性将产生影响。这就涉及到算法稳定性问题。,算法稳定性,当在某节点上xi的yi值有大小为的扰动时，如果在其后的各节点上的值yi产生的偏差都不大于，则称这种方法是稳定的。稳定性不仅与算法有关，而且与方程中函数f(x,y)也有关，讨论起来比较复杂。为简单起见，通常只针对模型方程,来讨论。一般方程若局部线性化，也可化为上述形式。模型方程相对比较简单，若一个数值方法对模型方程是稳定的，并不能保证该方法对任何方程都稳定，但若某方法对模型方程都不稳定，也就很难用于其他方程的求解。,先考察显式Euler方法的稳定性。模型方程的Euler公式为,将上式反复递推后，可得,或,式中,要使yi有界，其充要条件为,即,由于，故有,（3.26）,可见，如欲保证算法的稳定，显式Euler格式的步长h的选取要受到式（3.26）的限制。的绝对值越大，则限制的h值就越小。,用隐式Euler格式，对模型方程的计算公式为，可化为,由于,则恒有,故恒有,因此，隐式Euler格式是绝对稳定的（无条件稳定）的（对任何h0）。,常微分方程初值问题的求解，是将微分方程转化为差分方程来求解，并用计算值yi来近似替代y(xi)，这种近似替代是否合理，还须看分割区间的长度h越来越小时，即时，是否成立。若成立，则称该方法是收敛的，否则称为不收敛。这里仍以Euler方法为例，来分析其收敛性。Euler格式,算法收敛性,设表示取时,按Euler公式的计算结果,即,Euler方法局部截断误差为：,设有常数，则,（3.27）,总体截断误差,（3.28）,又,由于f(x,y)关于y满足李普希茨条件,即,代入(3.28)上式，有,再利用式（3.27），式（3.28）,即,上式反复递推后，可得,（3.29）,设（T为常数）,因为,所以,把上式代入式（3.29），得,若不计初值误差，即，则有,（3.30）,式(3.30)说明,当时,即,所以Euler方法是收敛的，且其收敛速度为，即具有一阶收敛速度。同时还说明Euler方法的整体截断误差为，因此算法的精度为一阶。,一阶常微分方程组与高阶方程我们已介绍了一阶常微分方程初值问题的各种数值解法，对于一阶常微分方程组，可类似得到各种解法，而高阶常微分方程可转化为一阶常微分方程组来求解。一阶常微分方程组对于一阶常微分方程组的初值问题,（3.31）,可以把单个方程中的f和y看作向量来处理，这样就可把前面介绍的各种差分算法推广到求一阶方程组初值问题中来。,式中,（3.34）,设为节点上的近似解，则有改进的Euler格式为,预报：,校正：,（3.32）,又，相应的四阶龙格库塔格式（经典格式）为,（3.33）,把节点xi上的yi和zi值代入式(7.34),依次算出,然后把它们代入式(7.33),算出节点xi+1上的yi+1和zi+1值。对于具有三个或三个以上方程的方程组的初值问题,也可用类似方法处理,只是更复杂一些而已。此外,多步方法也同样可以应用于求解方程组初值问题。,例7.6用改进的Euler法求解初值问题,取步长h=0.1，保留六位小数。,解:改进的Euler法公式为,预报：,校正：,将及h=0.1代入上式,得,由初值,计算得,高阶方程组高阶微分方程(或方程组)的初值问题,原则上都可以归结为一阶方程组来求解。例如,有二阶微分方程的初值问题,（3.35）,在引入新的变量后,即化为一阶方程组初值问题:,（3.36）,式（3.36）为一个一阶方程组的初值问题，对此可用前面中介绍的方法来求解。例如应用四阶龙格-库塔公式得,（3.37）,（3.38）,消去，上式简化为：,（3.39）,（3.40）,上述方法同样可以用来处