第9章常微分方程初值问题数值解法课件

第第9章章常微分方程初值问题数值解法常微分方程初值问题数值解法1.Euler公式公式2.改进的欧拉公式改进的欧拉公式3.龙格龙格库塔法库塔法4.亚当斯法亚当斯法5.算法的稳定性及收敛性算法的稳定性及收敛性2024/7/1219.1 9.1 引言引言 包含自变量、未知函数及未知函数的导数或微包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中分的方程称为微分方程。在微分方程中,自变量的自变量的个数只有一个个数只有一个,称为常微分称为常微分方程。方程。自变量的个数为自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数程的阶数。如果未知函数y y及其各阶导数及其各阶导数都是一次的都是一次的,则称它是线性的则称它是线性的,否则称为非线性的。否则称为非线性的。2024/7/122 在高等数学中，对于常微分方程的求解，给出在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。线性方程的解法等。但能求解的常微分方程仍然是但能求解的常微分方程仍然是有限的，大多数的常微分方程是有限的，大多数的常微分方程是不能不能给出解析解。给出解析解。譬如譬如 这个一阶微分方程就不能用初等函这个一阶微分方程就不能用初等函数来数来表达表达它的解。它的解。2024/7/123再如，方程再如，方程 的解的解 ,虽然有表可查虽然有表可查,但对于表但对于表上没有给出上没有给出 的值的值,仍需插值方法来仍需插值方法来计算计算2024/7/124 从从实际问题当中归纳出来的微分方程，通常主实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分要依靠数值解法来解决。本章主要讨论一阶常微分方程初值问题方程初值问题(9.1)在区间在区间a x ba x b上的数值解法上的数值解法。可以证明可以证明,如果函数在带形区域如果函数在带形区域 R:R:axbaxb,-y y 内连续，且关于内连续，且关于y y满足李普希兹满足李普希兹(Lipschitz)(Lipschitz)条件，即存在条件，即存在常数常数 L(L(它与它与x,yx,y无关无关)使使对对R内任内任意意x x及两及两个个 都成立都成立,则方程则方程(9.1)的的解解 在在 a,b 上上存在且唯一存在且唯一。2024/7/125数值方法的基本思想数值方法的基本思想 对常微分方程初值问题对常微分方程初值问题(9.1)式的数值解法，就式的数值解法，就是要算出精确解是要算出精确解y(x)y(x)在区间在区间 a,b 上的一系列离散节上的一系列离散节点点 处的函数值处的函数值 的近似值的近似值相相邻两个节点的间距邻两个节点的间距 称为步长，步称为步长，步长可以相等，也可以不等。本章总是假定长可以相等，也可以不等。本章总是假定h为定数，为定数，称为称为定步长定步长，这时节点可表示为，这时节点可表示为 数数值解法需要把连续性的问题加以离散化，从值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。而求出离散节点的数值解。2024/7/126 对常微分方程数值解法的基本出发点就是离散对常微分方程数值解法的基本出发点就是离散化。其数值解法有化。其数值解法有两个基本特两个基本特点，点，11它它们都采用们都采用“步进式步进式”，即求解过程顺着节点排列的次序一步一，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信步地向前推进，描述这类算法，要求给出用已知信息息 计计算算 的递推公的递推公式。式。22建建立立这类递推公式的基本方法是在这些节点上用数值积这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问分、数值微分、泰勒展开等离散化方法，对初值问题题中的中的导数导数 进行不同的离散化进行不同的离散化处理处理。2024/7/127对于初值问题对于初值问题的数值解法，首先要解决的问题就是的数值解法，首先要解决的问题就是如何如何对微分方程对微分方程进行离散化，建立求数值解的递推公式。递推公式通进行离散化，建立求数值解的递推公式。递推公式通常有两类，一类是计算常有两类，一类是计算yi+1时只用到时只用到xi+1,xi和和yi，即前，即前一步的值，因此有了初值以后就可以逐步往下计算，一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为此类方法称为单步法单步法；其代表是；其代表是龙格龙格库塔库塔法。另一法。另一类是计算类是计算y yi+1i+1时，除用到时，除用到x xi+1i+1,x,xi i和和y yi i以外，还要以外，还要用到，用到，即即前面前面k步的值，此类方法称为步的值，此类方法称为多步法多步法；其代；其代表是表是亚亚当斯法。当斯法。2024/7/1289.2简单的数的数值方法与基本概念方法与基本概念9.2.1Euler公式公式 欧拉（欧拉（Euler）方法是解初值问题的最简）方法是解初值问题的最简单的数值方法。初值问题单的数值方法。初值问题的解的解y=y(x)y=y(x)为通为通过点过点 的一的一条积条积分曲线。分曲线。积分曲线上每积分曲线上每一点一点 的的切线的斜率切线的斜率 等于函数等于函数 在这点的值。在这点的值。2024/7/129Euler法的求解过程是法的求解过程是:从初从初始点始点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点的坐标。点的坐标。2024/7/1210同样同样,过过点点P1(x x1 1,y,y1 1),),作积分曲线作积分曲线y=y(x)y=y(x)的切线的切线交直线交直线x=xx=x2 2于于P2点点,切线切线 的的斜率斜率直线方程为直线方程为当当 时时,得得 2024/7/1211当当 时时,得得由此获得了由此获得了P P2 2的坐标。重复以上过程的坐标。重复以上过程,就可获得一系就可获得一系列的点列的点:P P1 1,P P1 1,P Pn n。对已求得点对已求得点以以 为为斜率作直线斜率作直线 取取2024/7/1212 从图形上看从图形上看,就获得了一条近似于曲线就获得了一条近似于曲线y=y(x)y=y(x)的折线的折线 。这样这样,从从x x0 0 出发出发逐个逐个算出算出2024/7/1213通常取通常取 (常常数数),),则则 微分方程的初值问题微分方程的初值问题Euler法的计算格法的计算格式为：式为：用折线近似于曲线得到的计算公式称为欧拉公式用折线近似于曲线得到的计算公式称为欧拉公式2024/7/1214 Euler格式格式还还可用可用数值微分数值微分、数值积数值积分分和和泰勒展泰勒展开开等方法得到。等方法得到。2024/7/1215 Euler格式格式用用泰泰勒展勒展开开方法得到。方法得到。2024/7/1216 Euler格式格式用用数值微数值微分分方法得到。方法得到。选择不同的计算方法计算上式的积分项选择不同的计算方法计算上式的积分项就就会得到不同的计算公式。会得到不同的计算公式。将方程将方程 的两端在区间的两端在区间 上积分得，上积分得，2024/7/1217 用左矩形方法计算积分项用左矩形方法计算积分项 代入代入(9.3)(9.3)式式,并用并用y yi i近似代替式中近似代替式中y(xy(xi i)即可得到即可得到向前欧拉（向前欧拉（EulerEuler）公式）公式 由于数值积分的矩形方法精度很低，所以由于数值积分的矩形方法精度很低，所以欧拉（欧拉（EulerEuler）公式当然很粗糙。）公式当然很粗糙。2024/7/1218例例9.1用欧拉法解初值问题用欧拉法解初值问题取步长取步长h=0.2,h=0.2,计算过程保留计算过程保留4 4位小数位小数 解解:当当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.46132024/7/1219clear;y=1,x=0,%初始化for n=1:10 y=1.1*y-0.2*x/y,x=x+0.1,endy=1 x=0 y=1.1000 x=0.1000y=1.1918 x=0.2000y=1.2774 x=0.3000y=1.3582 x=0.4000y=1.4351 x=0.5000y=1.5090 x=0.6000y=1.5803 x=0.7000y=1.6498 x=0.8000y=1.7178 x=0.9000y=1.7848 x=1.00002024/7/1220对对方程方程 的两端在的两端在区间区间 上积上积分得，分得，代入代入(9.4)(9.4)式式,即即可得到梯形公式可得到梯形公式 由由于于数数值值积积分分的的梯梯形形公公式式比比矩矩形形公公式式的的精精度度高高，因因此此梯梯形形公公式式（9.59.5）是是比比欧欧拉拉公公式式(9.2 9.2)精精度度高高的的一一个数值方法。个数值方法。为了提高精度为了提高精度,改改用梯形方法计算其积分项，即用梯形方法计算其积分项，即 9.2.2梯形公式梯形公式2024/7/1221(9.5)(9.2)和和(9.5)式粗看差不多式粗看差不多，(9.2)Euler公式梯形公式 但但(9.5)式的右端含有未知的式的右端含有未知的y yi+1i+1,它是一个关它是一个关于于y yi+1i+1的函数方程的函数方程,这类数值方法称为这类数值方法称为隐式方法隐式方法。相反地相反地,欧拉法是关于欧拉法是关于y yi+1i+1的一个直接的计算公式，的一个直接的计算公式，这类数值方法称为这类数值方法称为显式方法显式方法。2024/7/12229.2.3两步欧拉公式两步欧拉公式对方程对方程 的两端在区间上的两端在区间上 积分得积分得 (9.6)改用中矩形公式计算其积分项，即改用中矩形公式计算其积分项，即 代入上式代入上式,并用并用y yi i近似代替式中近似代替式中y(xy(xi i)即可得到两步欧即可得到两步欧拉公式拉公式 2024/7/1223 前面介绍过的数值方法前面介绍过的数值方法,无论是欧拉方无论是欧拉方法法,还是梯形方法，它们都是单步法还是梯形方法，它们都是单步法,其特点其特点是在计算是在计算y yi+1i+1时只用到前一步的信息时只用到前一步的信息y yi i;可是可是公式公式(9.7)中除了中除了y yi i外外,还用到更前一步的信还用到更前一步的信息息y yi-1i-1,即调用了前两步的信息即调用了前两步的信息,故称其为两故称其为两步欧拉公式步欧拉公式2024/7/1224定定义义9.1在在yi准准确确的的前前提提下下,即即时时,用用数数值值方方法法计计算算yi+1的的误误差差,称称为为该该数数值值方方法法计计算时算时yi+1的局部截断误的局部截断误差。差。9.2.4欧拉法的局部截断误差欧拉法的局部截断误差衡量求解公式好坏的一个主要标准是求解公式的衡量求解公式好坏的一个主要标准是求解公式的精度精度,因此引入局部截断误差和阶数的概念。因此引入局部截断误差和阶数的概念。2024/7/1225假定假定由欧由欧拉公式拉公式，则有，则有因此有因此有 局部截断误差局部截断误差2024/7/1226定义定义9.2如果数如果数值方法的局部截断误差为值方法的局部截断误差为,则则称这种数值方法的称这种数值方法的阶数是阶数是P。显然，步长显然，步长(hN结束。结束。2024/7/1235（2）改改进进欧欧拉拉法法的的流流程程图图 2024/7/1236例例9.2 9.2 用改进欧拉法解初值问题用改进欧拉法解初值问题 区间为区间为 0,10,1,取步长取步长h=0.1h=0.1 解解:改进欧拉法的具体形式改进欧拉法的具体形式 本题的精确本题的精确解为解为2024/7/1237clearx=0,yn=1%初始化for n=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%校正end2024/7/1238例例9.3对初值问题对初值问题 证明用梯形公式求得的近似解为证明用梯形公式求得的近似解为并证明当步长并证明当步长h h0 0时时,y,yn n收敛于精确解收敛于精确解整理成显式整理成显式 反复迭代反复迭代,得到得到 证明证明:解初值问题的梯形公式为解初值问题的梯形公式为2024/7/1239由于由于 ，有，有 证毕证毕 2024/7/12409.3 9.3 龙格龙格-库塔（库塔（Runge-KuttaRunge-Kutta）法）法只有对平均斜率提供一种算法便可得到一种计算格式只有对平均斜率提供一种算法便可得到一种计算格式9.3.1龙格龙格-库塔法的基本思想库塔法的基本思想2024/7/1241Euler公式可改写成公式可改写成改进的改进的Euler公式又可改写成公式又可改写成用左端点的斜率作为平均斜率得到的计算格式用左端点的斜率作为平均斜率得到的计算格式用左、右端点斜率的平均值作为平均斜率得到的计用左、右端点斜率的平均值作为平均斜率得到的计算格式算格式2024/7/1242 上述两组公式在形式上有一个共同点上述两组公式在形式上有一个共同点:都是用都是用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)的值，的值，它是二阶方法。它的局部截断误差为它是二阶方法。它的局部截断误差为 。2024/7/1243 于是可考于是可考虑在虑在 内多内多预报几个点的预报几个点的斜率值，然后将其加权平均作为平均斜率斜率值，然后将其加权平均作为平均斜率，构，构造时要求近似公式在造时要求近似公式在(x xi i,y yi i)处的处的TaylorTaylor展开式展开式与解与解y y(x x)在在x xi i处的处的TaylorTaylor展开式的前面几项重展开式的前面几项重合，从而使近似公式达到所需要的阶数合，从而使近似公式达到所需要的阶数。则。则可可构造出更高精度的计算格式，这就是龙格构造出更高精度的计算格式，这就是龙格库库塔（塔（Runge-KuttaRunge-Kutta）法的基本思想。）法的基本思想。2024/7/1244只有多计算几个点的函数值，平均高度就更精确只有多计算几个点的函数值，平均高度就更精确.龙格龙格-库塔法的基本思想库塔法的基本思想多预报几个点的斜率值，将其加权平均作为平均斜率多预报几个点的斜率值，将其加权平均作为平均斜率.2024/7/12459.3.2二阶龙格二阶龙格库塔法库塔法 在在 上取两点上取两点xi i和和 ,以该两点处的以该两点处的斜率值斜率值k1 1和和k2 2的加权平均的加权平均(或称为线性组合或称为线性组合)来求取平来求取平均斜率均斜率k*的近似值的近似值K，即，即式中式中:k1 1为为xi i点处的切线斜率值，点处的切线斜率值，k2 2为为 点处的切线斜率值点处的切线斜率值,比照改进的欧拉比照改进的欧拉法法,将将 视为视为 ，即可得，即可得 对常微分方程初值问题对常微分方程初值问题(9.1)(9.1)式的解式的解 y=y(x),),根据微分根据微分中值定理，存在点中值定理，存在点 ，使得，使得 2024/7/1246式中式中 K K可看作是可看作是y=y(x)y=y(x)在区间在区间 上的平均斜率。上的平均斜率。所以可得计算公式为：所以可得计算公式为：（9.14）将将 在在x=xx=xi i处进行二阶处进行二阶TaylorTaylor展开：展开：（9 9.15）也即也即 （9.13）2024/7/1247将以上结果将以上结果代入代入2024/7/1248对式对式(9.15)(9.15)和和(9.16)(9.16)进行比较系数后可知进行比较系数后可知,只要只要 成立成立,格式格式(9.14)(9.14)的局部截断误差就等于的局部截断误差就等于2024/7/1249式式(9.17)(9.17)中具有三个未知量中具有三个未知量,但只有两个方程但只有两个方程,因而有因而有无穷多解。若取无穷多解。若取 ,则则P P=1=1，这是无穷多解中，这是无穷多解中的一个解，将以上所解的值代入式的一个解，将以上所解的值代入式(9.14)(9.14)并改写可得并改写可得 不难发现，上面的格式就是改进的欧拉格式。不难发现，上面的格式就是改进的欧拉格式。凡满足条件式（凡满足条件式（9.179.17）有一簇形如上式的计算格式，）有一簇形如上式的计算格式，这些格式统称为二阶龙格这些格式统称为二阶龙格库塔格式。因此改进的库塔格式。因此改进的欧拉格式是众多的二阶龙格欧拉格式是众多的二阶龙格库塔法中的一种特殊库塔法中的一种特殊格式。格式。2024/7/1250若取若取 ,则则 ，此时二阶龙格，此时二阶龙格-库塔库塔法的计算公式为法的计算公式为 此计算公式称为变形的二阶龙格此计算公式称为变形的二阶龙格库塔法。式中库塔法。式中 为区间为区间 的中点。的中点。2024/7/12519.3.3三阶龙格三阶龙格-库塔法库塔法 为了进一步提高精度，设除为了进一步提高精度，设除 外再增加一点外再增加一点 并用三个点并用三个点 的的斜率斜率k1 1,k2 2,k3 3加权平均加权平均得出平均斜率得出平均斜率k*的近似值，这时计算格式具有形式的近似值，这时计算格式具有形式:（9.18）为了预报点为了预报点 的斜率值的斜率值k3 3,在区间在区间 内有两内有两个斜率值个斜率值k1 1和和k2 2可以用可以用,可将可将k1 1,k2 2加权平均得出加权平均得出 上的平均斜率上的平均斜率,从而得到从而得到 的预报值的预报值 2024/7/1252于是可得于是可得运用运用Taylor展开方法选择参数展开方法选择参数 ,可以使格式可以使格式(9.18)的局部截断误差为的局部截断误差为 ,即具有即具有三阶三阶精度精度.2024/7/1253局局部截断误差为部截断误差为 ,即具有三阶精度，这类格式即具有三阶精度，这类格式统称为统称为三阶龙格三阶龙格库塔方库塔方法法。（9.19）是其中的一种，称为库塔（是其中的一种，称为库塔（KuttaKutta）公式。）公式。2024/7/12549.3.4四阶龙格四阶龙格库塔法库塔法如如果果需需要要再再提提高高精精度度，用用类类似似上上述述的的处处理理方方法法，只只需需在在区区间间上上用用四四个个点点处处的的斜斜率率加加权权平平均均作作为为平平均均斜斜率率k*的的近近似似值值，构构成成一一系系列列四四阶阶龙龙格格库库塔塔公公式式。具有四阶精度，即局部截断误差是具有四阶精度，即局部截断误差是。由于推导复杂，这里从略，只介绍最常用的一种由于推导复杂，这里从略，只介绍最常用的一种四阶经四阶经典典RK公公式式。（9.20）2024/7/12559.3.5 9.3.5 四阶龙格四阶龙格库塔法算法实现库塔法算法实现(1)(1)计算步骤计算步骤 输入输入 ,h,Nh,N 使用龙格使用龙格库塔公式（库塔公式（9.209.20）计算出）计算出y y1 1 输出输出 ，并使，并使 转到转到 直至直至n n N N 结束。结束。2024/7/1256(2 2)四四阶阶龙龙格格库库塔塔算算法法流流程程图图2024/7/1257(3)(3)程序实现程序实现(四阶龙格四阶龙格-库塔法计库塔法计(4)(4)算常微分方程初值问题算常微分方程初值问题)例例9.4 9.4 取步长取步长h=0.2h=0.2，用经典格式求解初值问题，用经典格式求解初值问题 解解:由四阶龙格由四阶龙格-库塔公式可得库塔公式可得 2024/7/1258可同样进行其余可同样进行其余y yi i的计算。本例方程的解为的计算。本例方程的解为，从表中看到所求的数值解具有，从表中看到所求的数值解具有4位有效数字。位有效数字。龙龙格格库库塔法塔法的推导基于的推导基于Taylor展开方法，因而它展开方法，因而它要求所求的解具有较好的要求所求的解具有较好的光滑性。光滑性。如果解的光滑性如果解的光滑性差，那么，使用四阶龙格差，那么，使用四阶龙格库塔方法求得的数值解，库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法。在实际计算其精度可能反而不如改进的欧拉方法。在实际计算时，应当针对问题的具体特点选择合适的算法。时，应当针对问题的具体特点选择合适的算法。2024/7/12599.3.6变步长的龙格变步长的龙格-库塔法库塔法在在微微分分方方程程的的数数值值解解中中，选选择择适适当当的的步步长长是是非非常常重重要要的的。单单从从每每一一步步看看，步步长长越越小小，截截断断误误差差就就越越小小；但但随随着着步步长长的的缩缩小小，在在一一定定的的求求解解区区间间内内所所要要完完成成的的步步数数就就增增加加了了。这这样样会会引引起起计计算算量量的的增增大大，并并且且会会引引起起舍舍入入误误差差的的大大量量积积累累与与传传播播。因因此此微微分分方程数值解法也有选择步长的问题。方程数值解法也有选择步长的问题。以经典的四阶龙格以经典的四阶龙格-库塔法库塔法(9.20)为例。从节点为例。从节点x xi i出发，先以出发，先以h为步长求出一个近似值，记为为步长求出一个近似值，记为 ，由于局部截断误差为由于局部截断误差为 ，故有，故有当h h值不大时，式中的系数值不大时，式中的系数c c可近似地看作为常数。可近似地看作为常数。2024/7/1260然后将步长折半然后将步长折半,即以为即以为 步长步长,从节点从节点x xi i出发出发,跨跨两步到节点两步到节点x xi+1i+1,再求得一个近似值再求得一个近似值 ,每跨一步的每跨一步的截断误差是截断误差是 ,因此有因此有这样这样 由此可得由此可得 这表明以这表明以 作为作为 的近似值，其误差可用步的近似值，其误差可用步长折半前后两次计算结果的偏差长折半前后两次计算结果的偏差 来判断所选步长是否适当来判断所选步长是否适当2024/7/1261当要求的数值精度为当要求的数值精度为时：时：（1 1）如如果果，反反复复将将步步长长折折半半进进行行计计算算，直直至至为止为止,并取其最后一次步长的计算结果作为并取其最后一次步长的计算结果作为 （2 2）如如果果为止，并以上一次步长的计算结果作为为止，并以上一次步长的计算结果作为 。这这种种通通过过步步长长加加倍倍或或折折半半来来处处理理步步长长的的方方法法称称为为变变步步长长法法。表表面面上上看看，为为了了选选择择步步长长，每每一一步步都都要要反反复复判判断断，增增加加了了计计算算工工作作量量，但但在在方方程程的的解解y(x)y(x)变变化化剧剧烈烈的的情情况况下下，总总的的计计算算工工作作量量得得到到减减少少，结果还是合算的。结果还是合算的。2024/7/12629.4算法的稳定性及收敛性算法的稳定性及收敛性9.4.1稳定性稳定性稳稳定定性性在在微微分分方方程程的的数数值值解解法法中中是是一一个个非非常常重重要要的的问问题题。因因为为微微分分方方程程初初值值问问题题的的数数值值方方法法是是用用差差分分格格式式进进行行计计算算的的，而而在在差差分分方方程程的的求求解解过过程程中中，存存在在着着各各种种计计算算误误差差，这这些些计计算算误误差差如如舍舍入入误误差差等等引引起起的的扰扰动动，在在传传播播过过程程中中，可可能能会会大大量量积积累累，对对计计算算结结果果的的准准确确性性将将产产生生影影响响。这这就就涉涉及及到到算算法法稳定性问题。稳定性问题。2024/7/1263 当当在在某某节节点点上上x xi i的的y yi i值值有有大大小小为为的的扰扰动动时时，如如果果在在其其后后的的各各节节点点 上上的的值值y yi i产产生生的的偏偏差差都都不大于不大于，则称这种方法是稳定的。，则称这种方法是稳定的。稳定性不仅与算法有关，而且与方程中函数稳定性不仅与算法有关，而且与方程中函数f(x,y)f(x,y)也有关，讨论起来比较复杂。也有关，讨论起来比较复杂。为简单起见，为简单起见，通常只针对模型方程通常只针对模型方程来来讨讨论论。一一般般方方程程若若局局部部线线性性化化，也也可可化化为为上上述述形形式式。模模型型方方程程相相对对比比较较简简单单，若若一一个个数数值值方方法法对对模模型型方方程程是是稳稳定定的的，并并不不能能保保证证该该方方法法对对任任何何方方程程都都稳稳定定，但但若若某某方方法法对对模模型型方方程程都都不不稳稳定定，也也就就很很难难用于其他方程的求解。用于其他方程的求解。2024/7/1264先考察显式先考察显式EulerEuler方法的稳定性。模型方程方法的稳定性。模型方程的的EulerEuler公式为公式为 将上式反复递推后，可得将上式反复递推后，可得 或或式中式中2024/7/1265 要使要使y yi i有界，其充要条件为有界，其充要条件为 即即 由于由于 ，故有，故有 （9.269.26）可见，如欲保证算法的稳定，显式可见，如欲保证算法的稳定，显式EulerEuler格式的格式的步长步长h h的选取要受到式（的选取要受到式（9.269.26）的限制。）的限制。的绝对的绝对值越大，则限制的值越大，则限制的h h值就越小。值就越小。用隐式用隐式EulerEuler格式，对模型方程格式，对模型方程 的计算公式为，可化为的计算公式为，可化为2024/7/1266由于由于 ,则恒有则恒有 ,故恒有故恒有 因此，隐式因此，隐式EulerEuler格式是绝对稳定的（无条件稳格式是绝对稳定的（无条件稳定）的（对任何定）的（对任何h0h0）。）。9.4.2 9.4.2 收敛性收敛性 常常微微分分方方程程初初值值问问题题的的求求解解，是是将将微微分分方方程程转转化化为为差差分分方方程程来来求求解解，并并用用计计算算值值y yi i来来近近似似替替代代y(xy(xi i)，这这种种近近似似替替代代是是否否合合理理，还还须须看看分分割割区区间间 的的长长度度h h越越来来越越小小时时，即即 时时，是是否否成成立立。若若成立，则称该方法是收敛的，否则称为不收敛。成立，则称该方法是收敛的，否则称为不收敛。这里仍以这里仍以EulerEuler方法为例，来分析其收敛性。方法为例，来分析其收敛性。EulerEuler格式格式2024/7/1267设设表示取表示取时时,按按Euler公式的计算结果公式的计算结果,即即Euler方法局部截断误差为：方法局部截断误差为：设有常数设有常数 ，则，则 （9.27）总体截断误差总体截断误差 （9.28）又又 由于由于f(x,y)f(x,y)关于关于y y满足李普希茨条件满足李普希茨条件,即即 2024/7/1268代入代入(9.28)上式，有上式，有再利用式（再利用式（9.27），式（），式（9.28）即即上式反复递推后，可得上式反复递推后，可得（9.29）2024/7/1269设设 （T T为常数）为常数）因为因为 所以所以 把上式代入式（把上式代入式（9.299.29），得），得 若不计初值误差，即若不计初值误差，即 ，则有，则有 （9.30）式式(9.30)(9.30)说明说明,当时当时 ,即即 ,所以所以EulerEuler方法是收敛的，且其收敛速度为方法是收敛的，且其收敛速度为 ，即具，即具有一阶收敛速度。同时还说明有一阶收敛速度。同时还说明EulerEuler方法的整体截断方法的整体截断误差为误差为 ，因此算法的精度为一阶。，因此算法的精度为一阶。2024/7/1270龙龙格格-库库塔塔方方法法是是一一类类重重要要算算法法，但但这这类类算算法法在在每每一一步步都都需需要要先先预预报报几几个个点点上上的的斜斜率率值值，计计算算量量比比较大。考虑到计算较大。考虑到计算yi+1之前已得出一系列之前已得出一系列节点节点上上的的斜斜率率值值，能能否否利利用用这这些些已已知知的的斜率值斜率值来来减少计算量呢？减少计算量呢？这就是亚当姆斯（这就是亚当姆斯（Adams）方法的设计思想。）方法的设计思想。9.5亚当姆斯方法亚当姆斯方法9.5.1亚当姆斯格式亚当姆斯格式2024/7/1271 设用设用x xi i,x,xi+1i+1两点的斜率值加权平均作为区间两点的斜率值加权平均作为区间 上的平均斜率，有计算格式上的平均斜率，有计算格式（9.21）选取参数选取参数，使格式（，使格式（9.219.21）具有二阶精度。）具有二阶精度。2024/7/1272将将 在在x xi i处处Taylor展开展开 代入计算代入计算格式格式 化简化简,因此有因此有 与与y(xi+1)在在xi处的处的Taylor展开式展开式相比较相比较,需取需取 才使格式才使格式(9.21)具有二阶精度。这样导出的计算格式具有二阶精度。这样导出的计算格式称之为二阶亚当姆斯格称之为二阶亚当姆斯格式。式。2024/7/1273四四阶亚当姆斯格式。阶亚当姆斯格式。(9.22)上述几种亚当姆斯格式都是显式的，算法比较上述几种亚当姆斯格式都是显式的，算法比较简单，但用节点简单，但用节点的斜率值来预报的斜率值来预报区间区间上的平均斜率是个上的平均斜率是个外推过程，外推过程，效果不够效果不够理想。为了进一步改善精度，变外推为理想。为了进一步改善精度，变外推为内插内插，即，即增加节点增加节点xi+1的斜率值来得出的斜率值来得出上的平均斜率。上的平均斜率。譬如考察形如譬如考察形如类似地可以导出三阶亚当姆斯格式。类似地可以导出三阶亚当姆斯格式。2024/7/1274（9.23）的隐式格式的隐式格式,设设(9.23)右端的右端的Taylor展开有展开有 可见要使格式可见要使格式(9.23)(9.23)具有二阶精度具有二阶精度,需令需令 ,这样就可构造二阶隐式亚当姆斯格式这样就可构造二阶隐式亚当姆斯格式 其实是梯形格式。类似可导出三阶隐式亚当姆斯格式其实是梯形格式。类似可导出三阶隐式亚当姆斯格式 2024/7/1275和四阶隐式亚当姆斯格式和四阶隐式亚当姆斯格式（9.249.24）9.5.2亚当姆斯预报亚当姆斯预报-校正格式校正格式参照改进的欧拉格式的构造方法，以四阶亚当参照改进的欧拉格式的构造方法，以四阶亚当姆斯为例，将显式（姆斯为例，将显式（9.22）和隐式（）和隐式（9.24）相结合，）相结合，用显式公式做预报，再用隐式公式做校正，可构成用显式公式做预报，再用隐式公式做校正，可构成亚当姆斯预报亚当姆斯预报-校正格式校正格式（9.259.25）预报：预报：校正：校正：2024/7/1276 这种预报这种预报-校正格式是四步法，它在计校正格式是四步法，它在计算算y yi+1i+1时不但用到前一步的信息时不但用到前一步的信息 ，而且，而且要用到再前面三步的信息要用到再前面三步的信息 ，因此，因此它不能自行启动。在实际计算时，可借助于它不能自行启动。在实际计算时，可借助于某种单步法，譬如四阶龙格某种单步法，譬如四阶龙格库塔法提供开库塔法提供开始值始值 。2024/7/1277例例9.5取步长取步长h=0.1,h=0.1,用亚当姆斯预报用亚当姆斯预报-校正公式求解校正公式求解 初值问题初值问题的数值解。的数值解。解解:用四阶龙格用四阶龙格-库塔公式求库塔公式求出初值出初值 ,计算得：计算得：表中的表中的 ,y,yi i和和y(xy(xi i)分别为预报值、校正值和准确分别为预报值、校正值和准确解解(),(),以比较计算结果的精度。以比较计算结果的精度。再使用亚当姆斯预报再使用亚当姆斯预报-校正公式校正公式(9.25),(9.25),2024/7/12789.6一阶常微分方程组与高阶方程一阶常微分方程组与高阶方程我我们们已已介介绍绍了了一一阶阶常常微微分分方方程程初初值值问问题题的的各各种种数数值值解解法法，对对于于一一阶阶常常微微分分方方程程组组，可可类类似似得得到到各各种种解解法法，而而高高阶阶常常微微分分方方程程可可转转化化为为一一阶阶常常微微分分方方程组来求解。程组来求解。9.6.1一阶常微分方程组一阶常微分方程组对于一阶常微分方程组的初值问题对于一阶常微分方程组的初值问题（9.31）可以把单个方程可以把单个方程中的中的f和和y看作向量看作向量来处理，这样就可把前面介绍的各种差分算法推广来处理，这样就可把前面介绍的各种差分算法推广到求一阶方程组初值问题中来。到求一阶方程组初值问题中来。2024/7/1279设设 为节点上的近似解，为节点上的近似解，则有改进的则有改进的EulerEuler格式为格式为 预报：预报：校正：校正：（9.32）又，相应的四阶龙格又，相应的四阶龙格库塔格式（经典格式）为库塔格式（经典格式）为 （9.33）2024/7/1280式中式中（9.34）2024/7/1281把节点把节点xi上的上的yi和和zi值代入式值代入式(7.34),依次算出依次算出,然然后后把把它它们们代代入入式式(7.33),算算出节点出节点xi+1上的上的yi+1和和zi+1值。值。对于具有三个或三个以上方程的方程组的初值问对于具有三个或三个以上方程的方程组的初值问题题,也可用类似方法处理也可用类似方法处理,只是更复杂一些而已。此外只是更复杂一些而已。此外,多步方法也同样可以应用于求解方程组初值问题。多步方法也同样可以应用于求解方程组初值问题。例例7.6 7.6 用改进的用改进的EulerEuler法求解初值问题法求解初值问题 取步长取步长h=0.1h=0.1，保留六位小数。，保留六位小数。解解:改进的改进的EulerEuler法公式为法公式为2024/7/1282预报：预报：校正：校正：将将 及及h=0.1h=0.1代入上式代入上式,得得 2024/7/1283由初值由初值 ,计算得计算得2024/7/12849.6.2高阶方程组高阶方程组 高阶微分方程高阶微分方程(或方程组或方程组)的初值问题的初值问题,原则上都原则上都可以归结为一阶方程组来求解。例如可以归结为一阶方程组来求解。例如,有二阶微分方有二阶微分方程的初值问题程的初值问题（9.359.35）在引入新的变量在引入新的变量 后后,即化为一阶方程组初值问题即化为一阶方程组初值问题:（9.369.36）式（式（9.369.36）为一个一阶方程组的初值问题，对此可）为一个一阶方程组的初值问题，对此可用前面中介绍的方法来求解。例如应用四阶龙格用前面中介绍的方法来求解。例如应用四阶龙格-库塔公式得库塔公式得 2024/7/1285（9.379.37）（9.389.38）2024/7/1286消去消去 ，上式简化为：，上式简化为：（9.399.39）（9.409.40）上述方法同样可以用来处理三阶或更高阶的微分方上述方法同样可以用来处理三阶或更高阶的微分方程（或方程组）的初值问题程（或方程组）的初值问题 2024/7/1287例例9.7 9.7 求解下列二阶微分方程的初值问题求解下列二阶微分方程的初值问题 取步长取步长h=0.1 h=0.1 解解:先作变换：令先作变换：令 ，代入上式，得一阶方程组，代入上式，得一阶方程组 用用四四阶阶龙龙格格-库库塔塔方方法法求求解解,按按式式(9.37)(9.37)及及(9.38)(9.38)进行计算：进行计算：取步长取步长 ,时时 2024/7/12882024/7/1289然后计算然后计算 时的时的y y2 2和和z z2 2；依此类推，直到；依此类推，直到i=9i=9时的时的y y1010和和z z1010，即可得到，即可得到其数值解。其数值解。2024/7/1290本章小结本章小结 本本章章介介绍绍了了常常微微分分方方程程初初值值问问题题的的基基本本数数值值解解法法。包包括括单单步步法法和和多多步步法法。单单步步法法主主要要有有欧欧拉拉法法、改改进进欧欧拉拉法法和和龙龙格格库库塔塔方方法法。多多步步法法是是亚亚当当姆姆斯斯法法。它它们们都都是是基基于于把把一一个个连连续续的的定定解解问问题题离离散散化化为为一一个个差差分分方方程程来来求求解解，是是一一种种步步进进式式的的方方法法。用用多多步步法法求常微分方程的数值解可获得较高的精度。求常微分方程的数值解可获得较高的精度。实际应用时，选择合适的算法有一定的难度，实际应用时，选择合适的算法有一定的难度，既要考虑算法的简易性和计算量，又要考虑截断误既要考虑算法的简易性和计算量，又要考虑截断误差和收敛性、稳定性。差和收敛性、稳定性。2024/7/1291龙龙格格-库库塔塔法法较较为为常常用用，适适用用于于多多步步方方法法中中作作初初值值计计算算和和函函数数f(x,y)较较为为简简单单的的场场合合。四四阶阶标标准准龙龙格格库库塔塔法法精精度度高高，程程序序简简单单，易易于于改改变变步步长长，比比较较稳稳定定，也也是是一一个个常常用用的的方方法法，但但计计算算量量较较大大。当当函函数数f(x,y)较较为为复复杂杂，可可用用显显式式亚亚当当姆姆斯斯方方法法或或亚亚当当姆姆斯斯预预测测校校正正方方法法，不不仅仅计计算算量量较较小小，稳稳定定性性也也比比较好，但不易改变步长。较好，但不易改变步长。一一般般采采用用龙龙格格库库塔塔法法提提供供初初值值y1,y2,y3，然然后后用用亚亚当当姆姆斯斯外外推推公公式式求求得得预预测测值值，再再由由亚亚当当姆姆斯斯内内插插值值求求得得校校正正值值yi+1，如如此此求求得得的的值值近近似似程程度度好好且且节省计算量，是一种较好的方法。节省计算量，是一种较好的方法。2024/7/1292作业作业习题九习题九P3161、2、5(1)2024/7/1293Thank you very much!2024/7/1294