第七章常微分方程数值解精选文档.ppt

第七章常微分方程数值解本讲稿第一页，共九十三页第七章第七章常微分方程的数值解法常微分方程的数值解法问题提出问题提出 倒葫芦形状容器壁上的刻度问题倒葫芦形状容器壁上的刻度问题.对于如图所示圆柱形状对于如图所示圆柱形状容器壁上的容积刻度容器壁上的容积刻度,可以利用圆柱体体积公式可以利用圆柱体体积公式其中直径其中直径D D为常数为常数.由于体积由于体积V V与相对于容器底部的任与相对于容器底部的任意高度意高度H H的函数关系明确的函数关系明确,因此在容器上可以方便地标因此在容器上可以方便地标出容器刻度出容器刻度,而对于几何形状不是规则的容器而对于几何形状不是规则的容器,比如倒比如倒葫芦形状容器壁上如何标出刻度呢葫芦形状容器壁上如何标出刻度呢?下表是经过测量下表是经过测量得到部分容器高度与直径的关系得到部分容器高度与直径的关系.本讲稿第二页，共九十三页x1o根据上表的数据根据上表的数据,可以拟合出倒葫可以拟合出倒葫芦形状容器的图芦形状容器的图,建立如图所示的建立如图所示的坐标轴后坐标轴后,问题即为如何根据任意问题即为如何根据任意高度高度x x标出容器体积标出容器体积V V的刻度的刻度,由微由微元思想分析可知元思想分析可知H 0 0.2 0.4 0.6 0.8 1.0D 0.04 0.11 0.26 0.56 1.04 1.17其中其中x x表示高度表示高度,直径直径D D是高度是高度x x的函数的函数,记为记为D(x),D(x),因此因此得到如下微分方程初值问题得到如下微分方程初值问题本讲稿第三页，共九十三页第七章第七章常微分方程的数值解法常微分方程的数值解法只要求解上述方程只要求解上述方程,就可求出体积就可求出体积V V与高度与高度x x之间的函数关之间的函数关系系,从而可标出容器壁上容积的刻度从而可标出容器壁上容积的刻度,但问题是但问题是(7.0)(7.0)中的中的函数函数D(x)D(x)无解析表达式无解析表达式,我们根本无法求出其解析解我们根本无法求出其解析解.怎样求解不能用解析方法求解的微分方程初制值怎样求解不能用解析方法求解的微分方程初制值问题问题,就是本章要讨论的重点就是本章要讨论的重点.(7.0)本讲稿第四页，共九十三页第七章第七章常微分方程的数值解法常微分方程的数值解法7.1 7.1 引言引言 包含自变量、未知函数及未知函数的导数或微分的包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中方程称为微分方程。在微分方程中,自变量的个数只有自变量的个数只有一个一个,称为常微分方程。自变量的个数为两个或两个以称为常微分方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数数y y及其各阶导数及其各阶导数都是一次的都是一次的,则称它是线性的则称它是线性的,否则称为非线性的。否则称为非线性的。本讲稿第五页，共九十三页 在高等数学中，对于常微分方程的求解，给出在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变量了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。但能求解的常微分方程仍然是有限的，方程的解法等。但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。大多数的常微分方程是不可能给出解析解。譬如譬如 这个一阶微分方程就不能用初等函数及其积分来表这个一阶微分方程就不能用初等函数及其积分来表达它的解。达它的解。本讲稿第六页，共九十三页再如，方程再如，方程 的解的解 ,虽然有表可查虽然有表可查,但对于表但对于表上没有给出上没有给出 的值的值,仍需插值方法来计仍需插值方法来计算算本讲稿第七页，共九十三页从实际问题当中归纳出来的微分方程，通常主要依靠数值从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程初值问题解法来解决。本章主要讨论一阶常微分方程初值问题(7.1)在区间在区间a x ba x b上的数值解法上的数值解法。可以证明可以证明,如果函数在带形区域如果函数在带形区域 R=axb,R=axb,-y y内连续，且关于内连续，且关于y y满足李普希兹满足李普希兹(Lipschitz)(Lipschitz)条件，条件，即存在常数即存在常数L(L(它与它与x,yx,y无关无关)使使对对R内任意两个内任意两个 都成立都成立,则方程则方程(7.1)的解的解 在在 a,b 上存在且唯一。上存在且唯一。本讲稿第八页，共九十三页常微分方程表示方法常微分方程表示方法本讲稿第九页，共九十三页在微分方程中在微分方程中,自变量的个数只有一个自变量的个数只有一个,称为称为常微分方程常微分方程.自变量的个数为两个或两个以上的微分方程自变量的个数为两个或两个以上的微分方程叫偏微分方程。叫偏微分方程。本讲稿第十页，共九十三页例如例如本讲稿第十一页，共九十三页当当x=0时时,y=1,可得可得c=1特解特解当当x=0时时,y=1,可得可得c=-1特解特解两边积分两边积分通解通解常微分方程解法回顾常微分方程解法回顾-本讲稿第十二页，共九十三页例：设降落伞从跳伞塔下落后，所受空气阻力与速度例：设降落伞从跳伞塔下落后，所受空气阻力与速度成正比，并设降落伞离开跳伞塔时（成正比，并设降落伞离开跳伞塔时（t=0）速度为）速度为零，求降落伞下落速度与时间的函数关系。零，求降落伞下落速度与时间的函数关系。解：设降落伞下落速度为解：设降落伞下落速度为v(t),降落伞在空中下落时，同时降落伞在空中下落时，同时受到重力受到重力P与阻力与阻力R的作用，重力大小为的作用，重力大小为mg，方向与，方向与v一致；一致；阻力大小为阻力大小为kv（k为比例系数），方向与为比例系数），方向与v相反，从而降落相反，从而降落伞所受外力为伞所受外力为F=mg-kv根据牛顿第二定律根据牛顿第二定律F=ma(a为加速度），得函数为加速度），得函数v(t)应满足的方程为应满足的方程为按题意，初始条件为按题意，初始条件为分离变量后得分离变量后得两端积分两端积分考虑到考虑到mg-kv0即即或或将初始条件将初始条件代如入上式得代如入上式得于是所求特解为于是所求特解为本讲稿第十三页，共九十三页本讲稿第十四页，共九十三页7.2数值方法的基本思想数值方法的基本思想 对常微分方程初值问题对常微分方程初值问题(7.1)式的数值解法，就是要式的数值解法，就是要算出精确解算出精确解y(x)y(x)在区间在区间 a,b 上的一系列离散节点上的一系列离散节点 处的函数值处的函数值 的近似值的近似值。相邻两个节点的间距。相邻两个节点的间距 称为步长，步称为步长，步长可以相等，也可以不等。本章总是假定长可以相等，也可以不等。本章总是假定h为定数，称为为定数，称为定步长定步长，这时节点可表示为，这时节点可表示为数值解法需要把连续性的问题加以离散化，从而求出数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。离散节点的数值解。本讲稿第十五页，共九十三页 对常微分方程数值解法的基本出发点就是离散化。对常微分方程数值解法的基本出发点就是离散化。其数值解法有两个基本特点，它们都采用其数值解法有两个基本特点，它们都采用“步进式步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信息描述这类算法，要求给出用已知信息 计算计算 的递推公式。建立这的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题数值微分、泰勒展开等离散化方法，对初值问题中的导数中的导数 进行不同的离散化处理进行不同的离散化处理。本讲稿第十六页，共九十三页对于初值问题对于初值问题的数值解法，首先要解决的问题就是如何对微分方程的数值解法，首先要解决的问题就是如何对微分方程进行离散化，建立求数值解的递推公式。递推公式通进行离散化，建立求数值解的递推公式。递推公式通常有两类，一类是计算常有两类，一类是计算yi+1时只用到时只用到xi+1,xi和和yi，即前一即前一步的值，因此有了初值以后就可以逐步往下计算，此类方步的值，因此有了初值以后就可以逐步往下计算，此类方法称为单步法；其代表是法称为单步法；其代表是龙格龙格库塔法。另一类是计算库塔法。另一类是计算y yi+1i+1时，除用到时，除用到x xi+1i+1,x,xi i和和y yi i以外，还要用到以外，还要用到，即前面，即前面k步的值，此类方法称为步的值，此类方法称为多步法多步法；其代表；其代表是亚当斯法。是亚当斯法。本讲稿第十七页，共九十三页7.3欧拉（欧拉（Euler）法）法7.3.1Euler公式公式 欧拉（欧拉（Euler）方法是解初值问题的最简单）方法是解初值问题的最简单的数值方法。初值问题的数值方法。初值问题的解的解y=y(x)y=y(x)代表通过点代表通过点 的一条称之为微的一条称之为微分方程的积分曲线。积分曲线上每一点分方程的积分曲线。积分曲线上每一点 的切线的斜率的切线的斜率 等于函数等于函数 在这在这点的值。点的值。本讲稿第十八页，共九十三页Euler法的求解过程是法的求解过程是:从初始点从初始点P0(即点即点(x(x0 0,y,y0 0)出发出发,作积分曲线作积分曲线y=y(x)y=y(x)在在P0点上切线点上切线 (其斜率为其斜率为 ),),与与x=xx=x1 1直线直线相交于相交于P1点点(即点即点(x(x1 1,y,y1 1),),得到得到y y1 1作为作为y(xy(x1 1)的近似值的近似值,如上如上图所示。过点图所示。过点(x(x0 0,y,y0 0),),以以f(xf(x0 0,y,y0 0)为为斜率的切线方程为斜率的切线方程为 当当 时时,得得 这样就获得了这样就获得了P P1 1点的坐标。点的坐标。本讲稿第十九页，共九十三页同样同样,过过点点P1(x x1 1,y,y1 1),),作积分曲线作积分曲线y=y(x)y=y(x)的切线的切线交直线交直线x=xx=x2 2于于P2点点,切线切线 的斜率的斜率 =直线方程为直线方程为当当 时时,得得 本讲稿第二十页，共九十三页当当 时时,得得由此获得了由此获得了P P2 2的坐标。重复以上过程的坐标。重复以上过程,就可获得一系列的就可获得一系列的点点:P P1 1,P P1 1,P Pn n。对已求得点对已求得点以以 =为斜率作直线为斜率作直线 取取本讲稿第二十一页，共九十三页 从图形上看从图形上看,就获得了一条近似于曲线就获得了一条近似于曲线y=y(x)y=y(x)的折线的折线 。这样这样,从从x x0 0逐个算出逐个算出对应的数值解对应的数值解 本讲稿第二十二页，共九十三页通常取通常取 (常数常数),),则则Euler法的计算格式法的计算格式 i=0,1,n(7.2)还可用数值微分、数值积分法和泰勒展开法推导还可用数值微分、数值积分法和泰勒展开法推导EulerEuler格式。格式。以数值积分为例进行推导。以数值积分为例进行推导。将方程将方程 的两端在区间的两端在区间 上积分得，上积分得，选择不同的计算方法计算上式的积分项选择不同的计算方法计算上式的积分项 ,就会得到不同的计算公式。就会得到不同的计算公式。(7.3)本讲稿第二十三页，共九十三页 公式推导本讲稿第二十四页，共九十三页 本讲稿第二十五页，共九十三页 用左矩形方法计算积分项用左矩形方法计算积分项 代入代入(7.3)(7.3)式式,并用并用y yi i近似代替式中近似代替式中y(xy(xi i)即可得到向前即可得到向前欧拉（欧拉（EulerEuler）公式）公式 由于数值积分的矩形方法精度很低，所以欧拉由于数值积分的矩形方法精度很低，所以欧拉（EulerEuler）公式当然很粗糙。）公式当然很粗糙。本讲稿第二十六页，共九十三页 本讲稿第二十七页，共九十三页例例7.1用欧拉法解初值问题用欧拉法解初值问题取步长取步长h=0.2,h=0.2,计算过程保留计算过程保留4 4位小数位小数 解解:h=0.2,: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本讲稿第二十八页，共九十三页例例7.2用用Euler法求初值问题法求初值问题y=xy2y(0)=0的数值解的数值解（取取h=0.1n=5)解：解：f(x,y)=x-y2;x0=y0=0;h=0.1由由Euler法的递推公式得：法的递推公式得：yn+1=yn+0.1(xny2n)yn=0n=0,1,2,3,4,5由上式计算所得数据由上式计算所得数据n012345xn00.10.20.30.40.5yn000.010.029990.059990.0995本讲稿第二十九页，共九十三页二、隐式二、隐式Euler(尤拉尤拉)格式格式本讲稿第三十页，共九十三页7.3.2梯形公式梯形公式为了提高精度为了提高精度,对方程对方程 的两端在区间上的两端在区间上 积分得，积分得，改用梯形方法计算其积分项，即改用梯形方法计算其积分项，即 (7.4)代入代入(7.4)(7.4)式式,并用近似代替式中即可得到梯形公式并用近似代替式中即可得到梯形公式 (7.5)由由于于数数值值积积分分的的梯梯形形公公式式比比矩矩形形公公式式的的精精度度高高，因因此此梯梯形公式（形公式（7.57.5）比欧拉公式）比欧拉公式(7.2)(7.2)的精度高一个数值方法。的精度高一个数值方法。本讲稿第三十一页，共九十三页(7.5)(7.5)式的右端含有未知的式的右端含有未知的y yi+1i+1,它是一个关于它是一个关于y yi+1i+1的函的函数方程数方程,这类数值方法称为这类数值方法称为隐式方法隐式方法。相反地。相反地,欧拉法是欧拉法是关于关于y yi+1i+1的一个直接的计算公式，的一个直接的计算公式，这类数值方法称为这类数值方法称为显式方法。显式方法。本讲稿第三十二页，共九十三页例例7.3用梯形公式求下面初值问题的解用梯形公式求下面初值问题的解在在x=0.01上的值上的值y(0.01)y=yy(0)=1解：取解：取h=0.01,由梯形公式得由梯形公式得y=exe0.01=1.010050167本讲稿第三十三页，共九十三页7.3.3两步欧拉公式两步欧拉公式对方程对方程 的两端在区间上的两端在区间上 积分得积分得 (7.6)改用中矩形公式计算其积分项，即改用中矩形公式计算其积分项，即 代入上式代入上式,并用并用y yi i近似代替式中近似代替式中y(xy(xi i)即可得到两步欧拉即可得到两步欧拉公式公式 (7.7)本讲稿第三十四页，共九十三页 前面介绍过的数值方法前面介绍过的数值方法,无论是欧拉方法无论是欧拉方法,还是梯形方法，它们都是单步法还是梯形方法，它们都是单步法,其特点是在计其特点是在计算算y yi+1i+1时只用到前一步的信息时只用到前一步的信息y yi i;可是公式可是公式(7.7)中除了中除了y yi i外外,还用到更前一步的信息还用到更前一步的信息y yi-1i-1,即调用了前两步的信息即调用了前两步的信息,故称其为两步欧拉公式故称其为两步欧拉公式本讲稿第三十五页，共九十三页7.3.4欧拉法的局部截断误差欧拉法的局部截断误差衡衡量量求求解解公公式式好好坏坏的的一一个个主主要要标标准准是是求求解解公公式式的的精精度度,因此引入局部截断误差和阶数的概念。因此引入局部截断误差和阶数的概念。定定义义7.1在在yi准准确确的的前前提提下下,即即时时,用用数数值值方方法法计计算算yi+1的的误误差差,称称为为该该数数值值方方法法计计算算时时yi+1的的局局部截断误差。部截断误差。对于欧拉公式，假定对于欧拉公式，假定，则有，则有而将真解而将真解y(x)在在xi处按二阶泰勒展开处按二阶泰勒展开 因此有因此有 本讲稿第三十六页，共九十三页定定义义7.2数数值值方方法法的的局局部部截截断断误误差差为为,则则称称这这种种数数值值方方法法的的阶阶数数是是P。步步长长(hN结束。结束。本讲稿第四十一页，共九十三页（2）改改进进欧欧拉拉法法的的流流程程图图 本讲稿第四十二页，共九十三页(3)3)程序实现程序实现(见附录见附录AA-15改进欧拉法计算常微改进欧拉法计算常微 分方程初值问题分方程初值问题)例例7.2 7.2 用改进欧拉法解初值问题用改进欧拉法解初值问题 区间为区间为 0,10,1,取步长取步长h=0.1h=0.1 解解:改进欧拉法的具体形式改进欧拉法的具体形式 本题的精确解为本题的精确解为 ,计算见计算见P P158158列表所示列表所示 本讲稿第四十三页，共九十三页例例7.3对初值问题对初值问题 证明用梯形公式求得的近似解为证明用梯形公式求得的近似解为并证明当步长并证明当步长h h0 0时时,y,yn n收敛于精确解收敛于精确解证明证明:解初值问题的梯形公式为解初值问题的梯形公式为 整理成显式整理成显式 反复迭代反复迭代,得到得到 本讲稿第四十四页，共九十三页由于由于 ，有，有 证毕证毕 本讲稿第四十五页，共九十三页7.4 7.4 龙格龙格-库塔（库塔（Runge-KuttaRunge-Kutta）法）法7.4.1 7.4.1 龙格龙格-库塔库塔(Runge-Kutta)(Runge-Kutta)法的基本思想法的基本思想 Euler Euler公式可改写成公式可改写成 则则yi+1i+1的的表表达达式式与与y(xi+1i+1)的的TaylorTaylor展展开开式式的的前前两两项项完完全全相相同同,即局部截断误差为即局部截断误差为 。改进的改进的EulerEuler公式又可改写成公式又可改写成 本讲稿第四十六页，共九十三页 上述两组公式在形式上有一个共同点上述两组公式在形式上有一个共同点:都是用都是用f(x,y)f(x,y)在某些点上值的线性组合得出在某些点上值的线性组合得出y(xy(xi+1i+1)的近似值的近似值y yi+1i+1,而且增加计算的次数而且增加计算的次数f f(x x,y y)的次数的次数,可提高截断误差可提高截断误差的阶。如欧拉公式的阶。如欧拉公式:每步计算一次每步计算一次f f(x x,y y)的值的值,为一阶方为一阶方法。改进欧拉公式需计算两次法。改进欧拉公式需计算两次f f(x x,y y)的值，它是二阶的值，它是二阶方法。它的局部截断误差为方法。它的局部截断误差为 。本讲稿第四十七页，共九十三页 于是可考虑用函数于是可考虑用函数f(x,y)在若干点上的函数在若干点上的函数值的线性组合来构造近似公式，构造时要求近值的线性组合来构造近似公式，构造时要求近似公式在似公式在(xi i,yi i)处的处的Taylor展开式与解展开式与解y(x)在在xi i处处的的Taylor展开式的前面几项重合，从而使近似公式达展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求偏导到所需要的阶数。既避免求偏导,又提高了计算方法又提高了计算方法精度的阶数。或者说精度的阶数。或者说,在在 这一步内多预报几个点的斜率值，然后将这一步内多预报几个点的斜率值，然后将其加权平均作为平均斜率，则可构造出更高精其加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格度的计算格式，这就是龙格库塔（库塔（Runge-Kutta）法的基本思想。）法的基本思想。本讲稿第四十八页，共九十三页7.4.2二阶龙格二阶龙格库塔法库塔法 在在 上取两点上取两点xi i和和 ,以该两点处的斜率以该两点处的斜率值值k1 1和和k2 2的加权平均的加权平均(或称为线性组合或称为线性组合)来求取平均斜率来求取平均斜率k*的的近似值近似值K，即，即式中式中:k1 1为为xi i点处的切线斜率值，点处的切线斜率值，k2 2为为 点处的切线斜率值点处的切线斜率值,比照改进的欧拉法比照改进的欧拉法,将将 视为视为 ，即可得，即可得 对常微分方程初值问题对常微分方程初值问题(7.1)(7.1)式的解式的解 y=y(x),),根据微分中根据微分中值定理，存在点值定理，存在点 ，使得，使得 本讲稿第四十九页，共九十三页式中式中 K K可看作是可看作是y=y(x)y=y(x)在区间在区间 上的平均斜率。所以可上的平均斜率。所以可得计算公式为：得计算公式为：（7.14）将将y(xy(xi i)在在x=xx=xi i处进行二阶处进行二阶TaylorTaylor展开：展开：（7.15）也即也即 （7.13）本讲稿第五十页，共九十三页将将 在在x=xx=xi i处进行一阶处进行一阶TaylorTaylor展开：展开：将以上结果代入（将以上结果代入（7.147.14）得：）得：（7.16）对式对式(7.15)(7.15)和和(7.16)(7.16)进行比较系数后可知进行比较系数后可知,只要只要 （7.17）成立成立,格式格式(7.14)(7.14)的局部截断误差就等于的局部截断误差就等于有有2 2阶阶精度精度本讲稿第五十一页，共九十三页式式(7.17)(7.17)中具有三个未知量中具有三个未知量,但只有两个方程但只有两个方程,因而有因而有无穷多解。若取无穷多解。若取 ,则则p p=1=1，这是无穷多解中的一，这是无穷多解中的一个解，将以上所解的值代入式个解，将以上所解的值代入式(7.14)(7.14)并改写可得并改写可得 不难发现，上面的格式就是改进的欧拉格式。凡满不难发现，上面的格式就是改进的欧拉格式。凡满足条件式（足条件式（7.177.17）有一簇形如上式的计算格式，这些格）有一簇形如上式的计算格式，这些格式统称为二阶龙格式统称为二阶龙格库塔格式。因此改进的欧拉格式是众库塔格式。因此改进的欧拉格式是众多的二阶龙格多的二阶龙格库塔法中的一种特殊格式。库塔法中的一种特殊格式。本讲稿第五十二页，共九十三页若取若取 ,则则 ，此时二阶龙格，此时二阶龙格-库塔库塔法的计算公式为法的计算公式为 此计算公式称为变形的二阶龙格此计算公式称为变形的二阶龙格库塔法。式中库塔法。式中 为区间为区间 的中点。的中点。本讲稿第五十三页，共九十三页7.4.3三阶龙格三阶龙格-库塔法库塔法 为了进一步提高精度，设除为了进一步提高精度，设除 外再增加一点外再增加一点 并用三个点并用三个点 ,的斜率的斜率k1 1,k2 2,k3 3加权平均加权平均得出平均斜率得出平均斜率k*的近似值，这时计算格式具有形式的近似值，这时计算格式具有形式:（7.18）为了预报点为了预报点 的斜率值的斜率值k3 3,在区间在区间 内有两内有两个斜率值个斜率值k1 1和和k2 2可以用可以用,可将可将k1 1,k2 2加权平均得出加权平均得出 上的平均斜率上的平均斜率,从而得到从而得到 的预报值的预报值 本讲稿第五十四页，共九十三页于是可得于是可得运用运用Taylor展开方法选择参数展开方法选择参数 ,可以使格式可以使格式(7.18)的局部截断误差为的局部截断误差为 ,即具有三阶精度，这类即具有三阶精度，这类格式统称为格式统称为三阶龙格三阶龙格库塔方法库塔方法。下列是其中的一种，。下列是其中的一种，称为称为库塔（库塔（Kutta）公式。）公式。（7.19）本讲稿第五十五页，共九十三页7.4.4四阶龙格四阶龙格库塔法库塔法如如果果需需要要再再提提高高精精度度，用用类类似似上上述述的的处处理理方方法法，只只需需在在区区间间上上用用四四个个点点处处的的斜斜率率加加权权平平均均作作为为平平均均斜斜率率k*的的近近似似值值，构构成成一一系系列列四四阶阶龙龙格格库库塔塔公公式式。具具有有四阶精度，即局部截断误差是四阶精度，即局部截断误差是。由于推导复杂，这里从略，只介绍最常用的一种由于推导复杂，这里从略，只介绍最常用的一种四阶经四阶经典龙格典龙格库塔公式库塔公式。（7.20）本讲稿第五十六页，共九十三页7.4.5 7.4.5 四阶龙格四阶龙格库塔法算法实现库塔法算法实现(1)(1)计算步骤计算步骤 输入输入 ,h,Nh,N 使用龙格使用龙格库塔公式（库塔公式（7.207.20）计算出）计算出y y1 1 输出输出 ，并使，并使 转到转到 直至直至n n N N 结束。结束。本讲稿第五十七页，共九十三页(2 2)四四阶阶龙龙格格库库塔塔算算法法流流程程图图本讲稿第五十八页，共九十三页(3)(3)程序实现程序实现(见附录见附录AA-16四阶龙格四阶龙格-库塔法计库塔法计 算常微分方程初值问题算常微分方程初值问题)例例7.4 7.4 取步长取步长h=0.2h=0.2，用经典格式求解初值问题，用经典格式求解初值问题 解解:由四阶龙格由四阶龙格-库塔公式可得库塔公式可得 本讲稿第五十九页，共九十三页可同样进行其余可同样进行其余y yi i的计算。本例方程的解为的计算。本例方程的解为，数值解，数值解y yi i与准确解与准确解y(xy(xi i)的对照表见教材的对照表见教材P P163163所示，从表所示，从表中看到所求的数值解具有中看到所求的数值解具有4位有效数字。位有效数字。龙格龙格库塔方法的推导基于库塔方法的推导基于Taylor展开方法，因展开方法，因而它要求所求的解具有较好的光滑性。如果解的光滑而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么，使用四阶龙格性差，那么，使用四阶龙格库塔方法求得的数值解，库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法。在实际计算时，其精度可能反而不如改进的欧拉方法。在实际计算时，应当针对问题的具体特点选择合适的算法。应当针对问题的具体特点选择合适的算法。本讲稿第六十页，共九十三页7.4.6变步长的龙格变步长的龙格-库塔法库塔法在在微微分分方方程程的的数数值值解解中中，选选择择适适当当的的步步长长是是非非常常重重要要的的。单单从从每每一一步步看看，步步长长越越小小，截截断断误误差差就就越越小小；但但随随着着步步长长的的缩缩小小，在在一一定定的的求求解解区区间间内内所所要要完完成成的的步步数数就就增增加加了了。这这样样会会引引起起计计算算量量的的增增大大，并并且且会会引引起起舍舍入入误误差差的的大大量量积积累累与与传传播播。因因此此微微分分方方程程数数值值解解法也有选择步长的问题。法也有选择步长的问题。以经典的四阶龙格以经典的四阶龙格-库塔法库塔法(7.20)为例。从节点为例。从节点x xi i出发出发，先以，先以h为步长求出一个近似值，记为为步长求出一个近似值，记为 ，由于局部截，由于局部截断误差为断误差为 ，故有，故有当h h值不大时，式中的系数值不大时，式中的系数c c可近似地看作为常数。可近似地看作为常数。本讲稿第六十一页，共九十三页然后将步长折半然后将步长折半,即以为即以为 步长步长,从节点从节点x xi i出发出发,跨两步到跨两步到节点节点x xi+1i+1,再求得一个近似值再求得一个近似值 ,每跨一步的截断误差是每跨一步的截断误差是 ,因此有因此有这样这样 由此可得由此可得 这表明以这表明以 作为作为 的近似值，其误差可用步长的近似值，其误差可用步长折半前后两次计算结果的偏差折半前后两次计算结果的偏差 来判断所选步长是否适当来判断所选步长是否适当本讲稿第六十二页，共九十三页当要求的数值精度为当要求的数值精度为时：时：（1 1）如如果果，反反复复将将步步长长折折半半进进行行计计算算，直直至至为止为止,并取其最后一次步长的计算结果作为并取其最后一次步长的计算结果作为 （2 2）如如果果为为止止，并以上一次步长的计算结果作为并以上一次步长的计算结果作为 。这这种种通通过过步步长长加加倍倍或或折折半半来来处处理理步步长长的的方方法法称称为为变变步步长长法法。表表面面上上看看，为为了了选选择择步步长长，每每一一步步都都要要反反复复判判断断，增增加加了了计计算算工工作作量量，但但在在方方程程的的解解y(x)y(x)变变化化剧剧烈烈的的情情况况下下，总总的计算工作量得到减少，结果还是合算的。的计算工作量得到减少，结果还是合算的。本讲稿第六十三页，共九十三页7.5亚当姆斯方法亚当姆斯方法7.5.1亚当姆斯格式亚当姆斯格式龙龙格格-库库塔塔方方法法是是一一类类重重要要算算法法，但但这这类类算算法法在在每每一一步步都都需需要要先先预预报报几几个个点点上上的的斜斜率率值值，计计算算量量比比较较大大。考考虑虑到到计计算算yi+1之之前前已已得得出出一一系系列列节节点点上上的的斜斜率率值值，能能否否利利用用这些已知值来减少计算量呢？这些已知值来减少计算量呢？这这就就是是亚亚当当姆姆斯斯（Adams）方方法法的的设设计计思思想。想。本讲稿第六十四页，共九十三页 设用设用x xi i,x,xi+1i+1两点的斜率值加权平均作为区间两点的斜率值加权平均作为区间 上的平均斜率，有计算格式上的平均斜率，有计算格式（7.21）选取参数选取参数，使格式（，使格式（7.217.21）具有二阶精度。）具有二阶精度。本讲稿第六十五页，共九十三页将将 在在x xi i处处Taylor展开展开 代入计算格式代入计算格式(7.21)(7.21)化简化简,并假设并假设,因此有因此有 与与y(xi+1)在在xi处的处的Taylor展开式展开式相比较相比较,需取需取 才使格式才使格式(7.21)具有二阶精度。这样导出的计算格式具有二阶精度。这样导出的计算格式称之为二阶亚当姆斯格式。类似地可以导出三阶亚当姆斯称之为二阶亚当姆斯格式。类似地可以导出三阶亚当姆斯格式格式。本讲稿第六十六页，共九十三页和和四阶亚当姆斯格式。四阶亚当姆斯格式。(7.22)这里和下面均记这里和下面均记上述几种亚当姆斯格式都是显式的，算法比较简单，上述几种亚当姆斯格式都是显式的，算法比较简单，但用节点但用节点的斜率值来预报区间的斜率值来预报区间上的上的平均斜率是个外推过程，效果不够理想。为了进一步改善平均斜率是个外推过程，效果不够理想。为了进一步改善精度，变外推为内插，即增加节点精度，变外推为内插，即增加节点xi+1的斜率值来得出的斜率值来得出上的平均斜率。譬如考察形如上的平均斜率。譬如考察形如本讲稿第六十七页，共九十三页（7.23）的隐式格式的隐式格式,设设(7.23)右端的右端的Taylor展开有展开有 可见要使格式可见要使格式(7.23)(7.23)具有二阶精度具有二阶精度,需令需令 ,这样就可构造二阶隐式亚当姆斯格式这样就可构造二阶隐式亚当姆斯格式 其实是梯形格式。类似可导出三阶隐式亚当姆斯格式其实是梯形格式。类似可导出三阶隐式亚当姆斯格式 本讲稿第六十八页，共九十三页和四阶隐式亚当姆斯格式和四阶隐式亚当姆斯格式（7.247.24）7.5.2亚当姆斯预报亚当姆斯预报-校正格式校正格式参照改进的欧拉格式的构造方法，以四阶亚当姆参照改进的欧拉格式的构造方法，以四阶亚当姆斯为例，将显式（斯为例，将显式（7.22）和隐式（）和隐式（7.24）相结合，用）相结合，用显式公式做预报，再用隐式公式做校正，可构成亚当显式公式做预报，再用隐式公式做校正，可构成亚当姆斯预报姆斯预报-校正格式校正格式（7.257.25）预报：预报：校正：校正：本讲稿第六十九页，共九十三页 这种预报这种预报-校正格式是四步法，它在计算校正格式是四步法，它在计算y yi+1i+1时不但用到前一步的信息时不但用到前一步的信息 ，而且要用，而且要用到再前面三步的信息到再前面三步的信息 ，因此它不能，因此它不能自行启动。在实际计算时，可借助于某种单步法，自行启动。在实际计算时，可借助于某种单步法，譬如四阶龙格譬如四阶龙格库塔法提供开始值库塔法提供开始值 。本讲稿第七十页，共九十三页例例7.5取步长取步长h=0.1,h=0.1,用亚当姆斯预报用亚当姆斯预报-校正公式求解校正公式求解 初值问题初值问题的数值解。的数值解。解解:用四阶龙格用四阶龙格-库塔公式求出发值库塔公式求出发值 ,计算得：计算得：表中的表中的 ,y ,yi i和和y(xy(xi i)分别为预报值、校正值和准确解分别为预报值、校正值和准确解(),),以比较计算结果的精度。以比较计算结果的精度。再使用亚当姆斯预报再使用亚当姆斯预报-校正公式校正公式(7.25),(7.25),见教材见教材P P166166列表列表算得其余的计算结果算得其余的计算结果本讲稿第七十一页，共九十三页7.6算法的稳定性及收敛性算法的稳定性及收敛性7.6.1稳定性稳定性稳稳定定性性在在微微分分方方程程的的数数值值解解法法中中是是一一个个非非常常重重要要的的问问题题。因因为为微微分分方方程程初初值值问问题题的的数数值值方方法法是是用用差差分分格格式式进进行行计计算算的的，而而在在差差分分方方程程的的求求解解过过程程中中，存存在在着着各各种种计计算算误误差差，这这些些计计算算误误差差如如舍舍入入误误差差等等引引起起的的扰扰动动，在在传传播播过过程程中中，可可能能会会大大量量积积累累，对对计计算算结结果果的的准准确确性性将将产生影响。这就涉及到算法稳定性问题。产生影响。这就涉及到算法稳定性问题。本讲稿第七十二页，共九十三页 当当在在某某节节点点上上x xi i的的y yi i值值有有大大小小为为的的扰扰动动时时，如如果果在在其其后后的的各各节节点点 上上的的值值y yi i产产生生的的偏偏差差都都不不大大于于，则称这种方法是稳定的。，则称这种方法是稳定的。稳定性不仅与算法有关，而且与方程中函数稳定性不仅与算法有关，而且与方程中函数f(x,y)f(x,y)也有也有关，讨论起来比较复杂。关，讨论起来比较复杂。为简单起见，通常只针对模型为简单起见，通常只针对模型方程方程来来讨讨论论。一一般般方方程程若若局局部部线线性性化化，也也可可化化为为上上述述形形式式。模模型型方方程程相相对对比比较较简简单单，若若一一个个数数值值方方法法对对模模型型方方程程是是稳稳定定的的，并并不不能能保保证证该该方方法法对对任任何何方方程程都都稳稳定定，但但若若某某方方法法对对模模型型方方程程都都不不稳稳定定，也也就就很很难难用用于于其其他他方方程程的的求求解。解。本讲稿第七十三页，共九十三页先考察显式先考察显式EulerEuler方法的稳定性。模型方程方法的稳定性。模型方程的的EulerEuler公式为公式为 将上式反复递推后，可得将上式反复递推后，可得 或或式中式中本讲稿第七十四页，共九十三页 要使要使y yi i有界，其充要条件为有界，其充要条件为 即即 由于由于 ，故有，故有 （7.267.26）可见，如欲保证算法的稳定，显式可见，如欲保证算法的稳定，显式EulerEuler格式的步格式的步长长h h的选取要受到式（的选取要受到式（7.267.26）的限制。）的限制。的绝对值越大，的绝对值越大，则限制的则限制的h h值就越小。值就越小。用隐式用隐式EulerEuler格式，对模型方程格式，对模型方程 的计算公式为，可化为的计算公式为，可化为本讲稿第七十五页，共九十三页由于由于 ,则恒有则恒有 ,故恒有故恒有 因此，隐式因此，隐式EulerEuler格式是绝对稳定的（无条件稳定）格式是绝对稳定的（无条件稳定）的（对任何的（对任何h0h0）。）。7.6.2 7.6.2 收敛性收敛性 常常微微分分方方程程初初值值问问题题的的求求解解，是是将将微微分分方方程程转转化化为为差差分分方方程程来来求求解解，并并用用计计算算值值y yi i来来近近似似替替代代y(xy(xi i)，这这种种近近似似替替代代是是否否合合理理，还还须须看看分分割割区区