常微分方程初值问题数值解法.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,第9章 常微分方程初值问题数值解法,9.1 引言,包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中, 自变量的个数只有一个, 称为常微分方程.。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数y及其各阶导数,都是一次的,则称它是线性的,否则称为非线性的。,在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。,但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。 譬如,这个一阶微分方程就不能用初等函数及其积分来表达它的解。,再如，方程,的解 ,虽然有表可查,但对于表上没有给出 的值,仍需插值方法来计算,从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程初值问题,( 9.1 ),在区间a x b上的数值解法,。,可以证明,如果函数在带形区域 R=axb,-y内连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L(它与x,y无关)使,对,R,内任意两个 都成立,则方程,( 9.1 ),的解,在,a, b,上存在且唯一。,数值方法的基本思想,对常微分方程初值问题(,9.1),式的数值解法，就是要算出精确解y(x)在区间,a,b,上的一系列离散节点 处的函数值,的近似值,。相邻两个节点的间距 称为步长，步,长可以相等，也可以不等。本章总是假定,h,为定数，称为,定步长,，这时节点可表示为,数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。,对常微分方程数值解法的,基本出发点就是离散化。其数值解法有两个基本特点，它们都采用“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信息,计算 的递推公式,。建立这类递推公式的基本方法是在这些节点上用,数值积分、数值微分、泰勒展开等离散化方法,，对初值问题,中的导数 进行不同的离散化处理,。,对于初值问题,的数值解法，首先要解决的问题就是,如何,对微分方程进行离散化，建立求数值解的递推公式。递推公式通常有两类，一类是计算y,i+1,时只用到x,i+1, x,i,和y,i,，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为,单步法,；其代表是,龙格库塔,法。另一类是计算y,i+1,时，除用到x,i+1,x,i,和y,i,以外，还要用到,，即前面,k,步的值，此类方法称为,多步法,；其代表,是亚当斯法。,9.2,简单的数值方法与基本概念,9.2.1 Euler公式,欧拉（,Euler,）方法是解初值问题的最简单的数值方法。初值问题,的解y=y(x)代表通过点 的一条称之为微分方程的积分曲线。积分曲线上每一点,的切线的斜率 等于函数 在这点的值。,Euler,法的求解过程是:从初始点,P,0,(,即点(x,0,y,0,)出发,作积分曲线y=y(x)在,P,0,点上切线 (其斜率为,),与x=x,1,直线,相交于,P,1,点(即点(x,1,y,1,),得到y,1,作为y(x,1,)的近似值,如上图所示。过点(x,0,y,0,),以f(x,0,y,0,),为,斜率的切线方程为,当 时,得,这样就获得了P,1,点的坐标。,同样, 过,点,P,1,(,x,1,y,1,),作积分曲线y=y(x)的切线,交直线x=x,2,于,P,2,点,切线 的斜率 =,直线方程为,当 时,得,当 时,得,由此获得了P,2,的坐标。重复以上过程,就可获得一系列的点:,P,1,P,1,P,n,。,对已求得点,以 = 为斜率作直线,取,从图形上看,就获得了一条近似于曲线y=y(x),的折线 。,这样,从x,0,逐个算出,对应的数值解,通常取 (常数),则,Euler,法的计算格式,i,=0,1,n,( 9.2 ),还可用,数值微分,、,数值积分法,和,泰勒展开法,推导Euler格式。以数值积分为例进行推导。,将方程 的两端在区间 上积分得，,选择不同的计算方法计算上式的积分项,就会得到不同的计算公式。,(9.3),用左矩形方法计算积分项,代入(9.3)式,并用y,i,近似代替式中y(x,i,)即可得到向前欧拉（Euler）公式,由于数值积分的矩形方法精度很低，所以欧拉（Euler）公式当然很粗糙。,例,9.1,用欧拉法解初值问题,取步长h=0.2 ,计算过程保留4位小数,解: h=0.2, 欧拉迭代格式,当 k=0, x,1,=0.2时，已知x,0,=0，y,0,=1，有,y(0.2),y,1,=0.21(401)0.8,当 k=1, x,2,=0.4时，已知x,1,=0.2， y,1,=0.8，有,y(0.4), y,2,=0.20.8(40.20.8)0.6144,当 k=2, x,3,=0.6时，已知x,2,=0.4， y,2,=0.6144，有,y(0.6),y,3,=0.2,0.6144,(4-0.4,0.6144)=0.4613,clear; y=1, x=0, %初始化,for n=1:10,y=1.1*y-0.2*x/y, x=x+0.1,end,y = 1 x = 0,y = 1.1000 x = 0.1000,y = 1.1918 x = 0.2000,y = 1.2774 x = 0.3000,y = 1.3582 x = 0.4000,y = 1.4351 x = 0.5000,y = 1.5090 x = 0.6000,y = 1.5803 x = 0.7000,y = 1.6498 x = 0.8000,y = 1.7178 x = 0.9000,y = 1.7848 x = 1.0000,9.2.2 梯形公式,为了提高精度,对方程 的两端在区间上,积分得，,改用梯形方法计算其积分项，即,(9.4 ),代入(7.4)式,并用近似代替式中即可得到梯形公式,( 9.5 ),由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式（9.5）比欧拉公式( 9.2 )的精度高一个数值方法。,( 9.5 ),(,9.5),式的右端含有未知的y,i+1,它是一个关于y,i+1,的函数方程,这类数值方法称为,隐式方法,。相反地,欧拉法是关于y,i+1,的一个直接的计算公式， 这类数值方法称为,显式方法。,9.2.3 两步欧拉公式,对方程 的两端在区间上 积分得,( 9.6 ),改用中矩形公式计算其积分项，即,代入上式,并用y,i,近似代替式中y(x,i,)即可得到两步欧拉公式,( 9.7 ),前面介绍过的数值方法,无论是欧拉方法,还是梯形方法，它们都是单步法,其特点是在计算y,i+1,时只用到前一步的信息y,i,;可是公式(,7.7),中除了y,i,外,还用到更前一步的信息y,i-1,即调用了前两步的信息,故称其为两步欧拉公式,9.2.4 欧拉法的局部截断误差,衡量求解公式好坏的一个主要标准是求解公式的精度, 因此引入局部截断误差和阶数的概念。,定义9.1,在y,i,准确,的前提下, 即 时, 用数值方法计算y,i+1,的误差 , 称为该数值方法计算时y,i+1,的局部截断误差。,对于欧拉公式，假定 ，则有,而将真解y(x)在x,i,处按二阶泰勒展开,因此有,定义9.2 数值方法的局部截断误差为 ,则称这种数值方法的阶数是P。,步长(h N,结束。,（,2,）改进欧拉法的流程图,(,3),程序实现(改进欧拉法计算常微,分方程初值问题,),例9.2 用改进欧拉法解初值问题,区间为,0,1,取步长h=0.1,解: 改进欧拉法的具体形式,本题的精确解为 ,clear,x=0,yn=1 %初始化,for n=1:10,yp=yn+0.1*(yn-2*x/yn); %预测,x=x+0.1;,yc=yn+0.1*(yp-2*x/yp) ;,yn=(yp+yc)/2 %校正,end,例,9.3,对初值问题,证明用梯形公式求得的近似解为,并证明当步长h,0时,y,n,收敛于精确解,证明: 解初值问题的梯形公式为,整理成显式,反复迭代,得到,由于 ，有,证毕,9.3 龙格-库塔（Runge-Kutta）法,9.3.1 龙格-库塔(Runge-Kutta)法的基本思想,Euler公式可改写成,则,y,i+1,的表达式,y,(,x,i+1,)与的Taylor展开式的前两项完全相同,即局部截断误差为 。,改进的Euler公式又可改写成,上述两组公式在形式上有一个共同点:都是用f(x,y)在某些点上值的线性组合得出y(x,i+1,)的近似值y,i+1,而且增加计算的次数,f,(,x,y,)的次数,可提高截断误差的阶。如欧拉公式:每步计算一次,f,(,x,y,)的值,为一阶方法。改进欧拉公式需计算两次,f,(,x,y,)的值，它是二阶方法。它的局部截断误差为 。,于是可考虑用函数,f,(,x,y,)在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式在(,x,i,y,i,)处的,Taylor,展开式与解,y,(,x,)在,x,i,处的,Taylor,展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求偏导,又提高了计算方法精度的阶数。或者说,在,这一步内多,预报几个点的斜率值，然后将其加权平均作为平均斜率,，则可构造出更高精度的计算格式，这就是龙格库塔（,Runge-Kutta,）法的基本思想。,9.3.2 二阶龙格库塔法,在 上取两点,x,i,和 ,以该两点处的斜率值,k,1,和,k,2,的加权平均(或称为线性组合)来求取平均斜率,k,*,的近似值,K,，即,式中:,k,1,为,x,i,点处的切线斜率值，,k,2,为 点处的切线斜率值,比照改进的欧拉法,将 视为 ，即可得,对常微分方程初值问题(9.1)式的解,y,=,y,(,x,),根据微分中值定理，存在点 ，使得,式中,K可看作是y=y(x)在区间 上的平均斜率。所以可得计算公式为：,（,9.14,）,将y(x,i,)在x=x,i,处进行二阶Taylor展开：,（9,.15,）,也即,（,9.13,）,将,在,x=x,i,处进行一阶,Taylor,展开：,将以上结果代入（9.14）得：,（,9.16,）,对式(9.15)和(9.16)进行比较系数后可知,只要,（,9.17,）,成立,格式(9.14)的局部截断误差就等于,有2阶,精度,式(9.17)中具有三个未知量,但只有两个方程,因而有无穷多解。若取 ,则,p,=1，这是无穷多解中的一个解，将以上所解的值代入式(9.14)并改写可得,不难发现，上面的格式就是改进的欧拉格式。凡满足条件式（9.17）有一簇形如上式的计算格式，这些格式统称为二阶龙格库塔格式。因此改进的欧拉格式是众多的二阶龙格库塔法中的一种特殊格式。,若取 ,则 ，此时二阶龙格-库塔,法的计算公式为,此计算公式称为变形的二阶龙格库塔法。式中,为区间 的中点。,9.3.3,三阶龙格-库塔法,为了进一步提高精度，设除 外再增加一点,并用三个点 , , 的斜率,k,1,k,2,k,3,加权平均,得出平均斜率,k,*,的近似值，这时计算格式具有形式:,（,9.18,）,为了预报点 的斜率值,k,3,在区间 内有两,个斜率值,k,1,和,k,2,可以用,可将,k,1,k,2,加权平均得出,上的平均斜率,从而得到 的预报值,于是可得,运用,Taylor,展开方法选择参数 ,可以使格式(,9.18),的局部截断误差为 ,即具有三阶精度，这类格式统称为,三阶龙格库塔方法,。下列是其中的一种，称为,库塔（Kutta）公式。,（,9.19,）,9.3.4,四阶龙格库塔法,如果需要再提高精度，用类似上述的处理方法，只需在区间 上用四个点处的斜率加权平均作为平均斜率,k,*,的近似值，构成一系列四阶龙格库塔公式。具有四阶精度，即局部截断误差是 。,由于推导复杂，这里从略，只介绍最常用的一种,四阶经典龙格库塔公式,。,（,9.20,）,9.3.5 四阶龙格库塔法算法实现,(1)计算步骤, 输入 ,h,N, 使用龙格库塔公式（9.20）计算出y,1, 输出 ，并使,转到 直至,n,N,结束。,(2) 四阶龙格库塔算法流程图,程序实现,(,四阶龙格,-,库塔法计,算常微分方程初值问题,),例9.4 取步长h=0.2，用经典格式求解初值问题,解:,由四阶龙格-库塔公式可得,可同样进行其余y,i,的计算。本例方程的解为,，从表中看到所求的数值解具有,4,位有效数字。,龙格库塔方法的,推导基于,Taylor,展开方法，因而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么，使用四阶龙格库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法,。在实际计算时，应当针对问题的具体特点选择合适的算法。,9.3.6 变步长的龙格-库塔法,在微分方程的数值解中，选择适当的步长是非常重要的。单从每一步看，,步长越小，截断误差就越小；但随着步长的缩小，在一定的求解区间内所要完成的步数就增加了。这样会引起计算量的增大，并且会引起舍入误差的大量积累与传播。因此微分方程数值解法也有选择步长的问题。,以经典的四阶龙格-库塔法(,9.20),为例。从节点x,i,出发，先以,h,为步长求出一个近似值，记为 ，由于局部截断误差为 ，故有,当,h值不大时，式中的系数c可近似地看作为常数。,然后将步长折半,即以为 步长,从节点x,i,出发,跨两步到节点x,i+1,再求得一个近似值 ,每跨一步的截断误差是 ,因此有,这样,由此可得,这表明以 作为 的近似值，其误差可用步长折半前后两次计算结果的偏差,来判断所选步长是否适当,当要求的数值精度为时：,（1）如果，反复将步长折半进行计算，直至0）。,9.4.2 收敛性,常微分方程初值问题的求解，是将微分方程转化为差分方程来求解，并用计算值y,i,来近似替代y(x,i,)，,这种近似替代是否合理，还须看分割区间,的长度h越来越小时，即 时， 是否成立。若成立，则称该方法是收敛的，否则称为不收敛。,这里仍以Euler方法为例，来分析其收敛性。,Euler格式,设 表示取 时, 按Euler公式的计算结果, 即,Euler,方法局部截断误差为：,设有常数 ，则,（9.27）,总体截断误差,（9.28）,又,由于f(x,y)关于y满足李普希茨条件,即,代入(9.28)上式，有,再利用式（9.27），式（9.28）,即,上式反复递推后，可得,（,9.29,）,设 （T为常数）,因为,所以,把上式代入式（9.29），得,若不计初值误差，即 ，则有,（,9.30,）,式(9.30)说明,当时 , ,即 ,所以Euler方法是收敛的，且其收敛速度为 ，即具有一阶收敛速度。同时还说明Euler方法的整体截断误差为 ，因此算法的精度为一阶。,9.5 亚当姆斯方法,9.5.1 亚当姆斯格式,龙格-库塔方法是一类重要算法，但这类算法在每一步都需要先预报几个点上的斜率值，计算量比较大。考虑到计算y,i+1,之前已得出一系列节点上 的斜率值，能否利用这些已知值来减少计算量呢？,这就是亚当姆斯（Adams）方法的设计思想。,设用x,i,x,i+1,两点的斜率值加权平均作为区间,上的平均斜率，有计算格式,（,9.21,）,选取参数，使格式（9.21）具有二阶精度。,将 在x,i,处,Taylor,展开,代入计算格式(9.21)化简,并假设,因此有,与y(x,i+1,)在x,i,处的Taylor展开式,相比较, 需取,才使格式(9.21)具有二阶精度。这样导出的计算格式,称之为二阶亚当姆斯格式。类似地可以导出三阶亚当姆斯格式,。,和,四阶亚当姆斯格式。,(,9.22),这里和下面均记,上述几种亚当姆斯格式都是显式的，算法比较简单，但用节点 的斜率值来预报区间 上的平均斜率是个,外推过程，,效果不够理想。为了进一步改善精度，变外推为,内插,，即增加节点x,i+1,的斜率值来得出 上的平均斜率。譬如考察形如,（,9.23,）,的隐式格式,设(,9.23),右端的,Taylor,展开有,可见要使格式(9.23)具有二阶精度,需令 ,这样就可构造二阶隐式亚当姆斯格式,其实是梯形格式。类似可导出三阶隐式亚当姆斯格式,和四阶隐式亚当姆斯格式,（9.24）,9.5.2 亚当姆斯预报-校正格式,参照改进的欧拉格式的构造方法，以四阶亚当姆斯为例，将显式（9.22）和隐式（9.24）相结合，用显式公式做预报，再用隐式公式做校正，可构成亚当姆斯预报-校正格式,（9.25）,预报：,校正：,这种预报-校正格式是四步法，它在计算y,i+1,时不但用到前一步的信息 ，而且要用到再前面三步的信息 ，因此它不能自行启动。在实际计算时，可借助于某种单步法，譬如四阶龙格库塔法提供开始值 。,例,9.5,取步长h=0.1,用亚当姆斯预报-校正公式求解,初值问题,的数值解。,解:,用四阶龙格-库塔公式求出发值 ,计算得：,表中的 ,y,i,和y(x,i,)分别为预报值、校正值和准确解( ),以比较计算结果的精度。,再使用亚当姆斯预报-校正公式(9.25),9.6 一阶常微分方程组与高阶方程,我们已介绍了一阶常微分方程初值问题的各种数值解法，对于一阶常微分方程组，可类似得到各种解法，而高阶常微分方程可转化为一阶常微分方程组来求解。,9.6.1 一阶常微分方程组,对于一阶常微分方程组的初值问题,（,9.31,）,可以把单个方程 中的f 和y看作向量来处理，这样就可把前面介绍的各种差分算法推广到求一阶方程组初值问题中来。,设 为节点上的近似解，,则有改进的Euler格式为,预报：,校正：,（,9.32）,又，相应的四阶龙格库塔格式（经典格式）为,（,9.33）,式中,（,9.34,）,把节点x,i,上的y,i,和z,i,值代入式(7.34), 依次算出, 然后把它们代入式(7.33), 算出节点x,i+1,上的y,i+1,和z,i+1,值。,对于具有三个或三个以上方程的方程组的初值问题,也可用类似方法处理,只是更复杂一些而已。此外,多步方法也同样可以应用于求解方程组初值问题。,例7.6 用改进的Euler法求解初值问题,取步长h=0.1，保留六位小数。,解: 改进的Euler法公式为,预报：,校正：,将 及h=0.1代入上式,得,由初值 ,计算得,9.6.2 高阶方程组,高阶微分方程(或方程组)的初值问题,原则上都,可以归结为一阶方程组来求解。例如,有二阶微分方,程的初值问题,（9.35）,在引入新的变量 后,即化为一阶方程组初值问题:,（9.36）,式（9.36）为一个一阶方程组的初值问题，对此可用前面中介绍的方法来求解。例如应用四阶龙格-库塔公式得,（9.37）,（9.38）,消去 ，上式简化为：,（9.39）,（9.40）,上述方法同样可以用来处理三阶或更高阶的微分方程（或方程组）的初值问题,例9.7 求解下列二阶微分方程的初值问题,取步长h=0.1,解:先作变换：令 ，代入上式，得一阶方程组,用四阶龙格-库塔方法求解,按式(9.37)及(9.38)进行计算：,取步长 , , ,时,然后计算 时的,y,2,和z,2,；依此类推，直到i=9时的y,10,和z,10,，即可得到其数值解。,本章小结,本章介绍了常微分方程初值问题的基本数值解法。包括单步法和多步法。单步法主要有欧拉法、改进欧拉法和龙格库塔方法。多步法是亚当姆斯法。它们都是基于把一个连续的定解问题离散化为一个差分方程来求解，是一种步进式的方法。用多步法求常微分方程的数值解可获得较高的精度。,实际应用时，选择合适的算法有一定的难度，既要考虑算法的简易性和计算量，又要考虑截断误差和收敛性、稳定性。,龙格-库塔法较为常用，适用于多步方法中作初值计算和函数f(x,y)较为简单的场合。四阶标准龙格库塔法精度高，程序简单，易于改变步长，比较稳定，也是一个常用的方法，但计算量较大。当函数f(x,y)较为复杂，可用显式亚当姆斯方法或亚当姆斯预测校正方法，不仅计算量较小，稳定性也比较好，但不易改变步长。,一般采用龙格库塔法提供初值y,1, y,2, y,3,，然后用亚当姆斯外推公式求得预测值 ，再由亚当姆斯内插值求得校正值y,i+1,，如此求得的值近似程度好且节省计算量，是一种较好的方法。,作业,习题九 P,316,1、2、5(1),Thank you very much!,