W Y第八章第八章常微分方程常微分方程数值解法数值解法1 W Y第八章目录第八章目录§1 欧拉（欧拉（Euler）方法）方法 1.1 Euler法及其简单改进法及其简单改进 1.2 改进的改进的Euler法法§2 龙格库塔（龙格库塔（Runge-kutta）方法）方法 2.1 龙格龙格-库塔方法的基本思想库塔方法的基本思想 2.2 二阶龙格二阶龙格-库塔公式库塔公式 2.3 高阶高阶R-K公式公式 2.4 变步长变步长R-K法法§3 线性多步法线性多步法§4 一阶方程组与高阶方程初值问题一阶方程组与高阶方程初值问题§5 收敛性与稳定性收敛性与稳定性 2 W Y第八章第八章 序序 许多科学技术问题，例如天文学中的星体运动，空间许多科学技术问题，例如天文学中的星体运动，空间许多科学技术问题，例如天文学中的星体运动，空间许多科学技术问题，例如天文学中的星体运动，空间技术中的物体飞行，自动控制中的系统分析，力学中的振技术中的物体飞行，自动控制中的系统分析，力学中的振技术中的物体飞行，自动控制中的系统分析，力学中的振技术中的物体飞行，自动控制中的系统分析，力学中的振动，工程问题中的电路分析等，都可归结为常微分方程的动，工程问题中的电路分析等，都可归结为常微分方程的动，工程问题中的电路分析等，都可归结为常微分方程的动，工程问题中的电路分析等，都可归结为常微分方程的初值问题。
初值问题初值问题初值问题 所谓初值问题，是函数及其必要的导数在积分的起始所谓初值问题，是函数及其必要的导数在积分的起始所谓初值问题，是函数及其必要的导数在积分的起始所谓初值问题，是函数及其必要的导数在积分的起始点为已知的一类问题，一般形式为：点为已知的一类问题，一般形式为：点为已知的一类问题，一般形式为：点为已知的一类问题，一般形式为： 我们先介绍我们先介绍我们先介绍我们先介绍简单的一阶问题：简单的一阶问题：简单的一阶问题：简单的一阶问题：3 W Y第八章第八章 序序由由由由常微分方程常微分方程常微分方程常微分方程的理论可知：上述问题的解唯一存在的理论可知：上述问题的解唯一存在的理论可知：上述问题的解唯一存在的理论可知：上述问题的解唯一存在 常微分方程常微分方程常微分方程常微分方程求解求什么？应求一满足初值问题（求解求什么？应求一满足初值问题（求解求什么？应求一满足初值问题（求解求什么？应求一满足初值问题（8—18—1））））的解函数的解函数的解函数的解函数y y = = y y( (x x) ) ，如对下列微分方程：，如对下列微分方程：，如对下列微分方程：，如对下列微分方程：《高等数学》中，《高等数学》中，《高等数学》中，《高等数学》中，微分方程求解，如对一阶微分方程：微分方程求解，如对一阶微分方程：微分方程求解，如对一阶微分方程：微分方程求解，如对一阶微分方程：y y =f(x,y) =f(x,y)是求解解函数是求解解函数是求解解函数是求解解函数y y = = y y( (x x) ) ，使满足上述方程。
但能，使满足上述方程但能，使满足上述方程但能，使满足上述方程但能够够够够求出准确的解析函数求出准确的解析函数求出准确的解析函数求出准确的解析函数y(x)y(x)的微分方程是很少的，《高数》的微分方程是很少的，《高数》的微分方程是很少的，《高数》的微分方程是很少的，《高数》中研究微分方程的求解，是中研究微分方程的求解，是中研究微分方程的求解，是中研究微分方程的求解，是分门别类讨论分门别类讨论分门别类讨论分门别类讨论，对不同类型的，对不同类型的，对不同类型的，对不同类型的微分方程，求解方法不一样，因此，要求解微分方程，微分方程，求解方法不一样，因此，要求解微分方程，微分方程，求解方法不一样，因此，要求解微分方程，微分方程，求解方法不一样，因此，要求解微分方程，首首首首先必须认清类型先必须认清类型先必须认清类型先必须认清类型4 W Y微分方程微分方程 数值解数值解 而而而而常微分方程常微分方程常微分方程常微分方程 初值问题的数值解法，是要寻求解函数初值问题的数值解法，是要寻求解函数初值问题的数值解法，是要寻求解函数初值问题的数值解法，是要寻求解函数y y( (x x) )在一系列点在一系列点在一系列点在一系列点y(xy(xi i ) ) （离散点）（离散点）（离散点）（离散点）: : 上上上上 y(x y(xi i ) )的的的的近似值近似值近似值近似值y yi i（（（（ i= i=1,2,…,n1,2,…,n），并且还可由这些），并且还可由这些），并且还可由这些），并且还可由这些（（（（x xi i, ,y yi i) )（（（（i=i=1,2,…,n1,2,…,n）构造插值函数作为近似函数。
上述离散点相）构造插值函数作为近似函数上述离散点相）构造插值函数作为近似函数上述离散点相）构造插值函数作为近似函数上述离散点相邻两点间的距离邻两点间的距离邻两点间的距离邻两点间的距离h hi i=x=xi i-1-1- -x xi i 称为步长，若称为步长，若称为步长，若称为步长，若h hi i 都相等为一定数都相等为一定数都相等为一定数都相等为一定数h, h, 则称为定步长，否则为变步长则称为定步长，否则为变步长则称为定步长，否则为变步长则称为定步长，否则为变步长 由于在实际问题和科学研究中遇到的微分方程往往很由于在实际问题和科学研究中遇到的微分方程往往很由于在实际问题和科学研究中遇到的微分方程往往很由于在实际问题和科学研究中遇到的微分方程往往很复杂，复杂，复杂，复杂，绝大多数很难，甚至不可能求出解析绝大多数很难，甚至不可能求出解析绝大多数很难，甚至不可能求出解析绝大多数很难，甚至不可能求出解析函数函数函数函数y y( (x x) )，因，因，因，因此只能考虑求其数值解此只能考虑求其数值解此只能考虑求其数值解此只能考虑求其数值解 本章重点讨论如下本章重点讨论如下本章重点讨论如下本章重点讨论如下一阶微分方程：一阶微分方程：一阶微分方程：一阶微分方程：在此基础上介绍一阶微分方程组与在此基础上介绍一阶微分方程组与在此基础上介绍一阶微分方程组与在此基础上介绍一阶微分方程组与高阶微分方程的数值解法。
高阶微分方程的数值解法高阶微分方程的数值解法高阶微分方程的数值解法5 W Y§1 欧拉（欧拉（Euler）法）法 以以Euler法及其改进方法为例法及其改进方法为例,说明说明常微分方程常微分方程初值问题数值解法的一般概初值问题数值解法的一般概念，念，Euler法很简单，准确度也不高，法很简单，准确度也不高，介绍此方法的目的，是由于对它的分析介绍此方法的目的，是由于对它的分析讨论能够比较清楚地显示出方法的一些讨论能够比较清楚地显示出方法的一些特点，而这些特点及基本方法反映了其特点，而这些特点及基本方法反映了其它方法的特点它方法的特点 Euler法用于求法用于求解一阶微分方解一阶微分方程初值问题：程初值问题：6 W Y1.1 Euler法及其简单改进法及其简单改进 EulerEuler公公公公式为：式为：式为：式为： 由由由由x x0 0出发出发出发出发x x1 1, ,x x2 2,…,,…,x xN N，而利用此式可算出对应的，而利用此式可算出对应的，而利用此式可算出对应的，而利用此式可算出对应的 y y1 1, ,y y2 2,…,,…,y yN N，式（，式（，式（，式（8-28-2）称为）称为）称为）称为差分方程差分方程差分方程差分方程（序列（序列（序列（序列{ {y yn n} }满足的满足的满足的满足的方方方方 程）。
程） 下面是下面是下面是下面是EulerEuler公式的推导公式的推导公式的推导公式的推导 ：：：：一一、从几何意义出发：、从几何意义出发：、从几何意义出发：、从几何意义出发：y y = =f f ( (x x, ,y y) )的解函数的解函数的解函数的解函数y=y=y y ( (x x) ) 在在在在xoyxoy平平平平面上是一族解曲线，面上是一族解曲线，面上是一族解曲线，面上是一族解曲线， 而初值问题则是其中一条积分曲线，而初值问题则是其中一条积分曲线，而初值问题则是其中一条积分曲线，而初值问题则是其中一条积分曲线，假定假定假定假定y y = = y y( (x x) )的曲线如图的曲线如图的曲线如图的曲线如图8-18-1从给定的点从给定的点从给定的点从给定的点P P0 0( (x x0 0, ,y y0 0) ) 出发，以出发，以出发，以出发，以P P0 0为切点，作切线，切线斜率为曲线为切点，作切线，切线斜率为曲线为切点，作切线，切线斜率为曲线为切点，作切线，切线斜率为曲线y y( (x x) )的切线斜率的切线斜率的切线斜率的切线斜率 y y = =f f ( (x x0 0, ,y y0 0) )，因此可，因此可，因此可，因此可 得切线：（点斜式）得切线：（点斜式）得切线：（点斜式）得切线：（点斜式） P1 P2y(x)P0 x2x1x07 W YEulerEuler公式的推导公式的推导公式的推导公式的推导( (续续续续1 1））））几何意义几何意义几何意义几何意义：用折线近似解曲线：用折线近似解曲线：用折线近似解曲线：用折线近似解曲线y y = = y y( (x x) )，折线不会偏离太远，折线不会偏离太远，折线不会偏离太远，折线不会偏离太远 ，因为每项以，因为每项以，因为每项以，因为每项以f f ( (x x, , y y) )（斜率）修正。
斜率）修正斜率）修正斜率）修正 切线与切线与切线与切线与x x = = x x1 1交于交于交于交于P P1 1( (x x1 1, ,y y1 1) )，在，在，在，在[ [x x0 0, ,x x1 1] ]上以切线上以切线上以切线上以切线 近似曲线，近似曲线，近似曲线，近似曲线， 8 W YEuler公式公式的推导（续的推导（续2））二二二二、利用、利用、利用、利用TaylorTaylor级数：将级数：将级数：将级数：将y y( (x x) )在在在在x xn n处展开：处展开：处展开：处展开： 9 W YEuler公式的推导（续公式的推导（续3））公式（公式（公式（公式（8-38-3）称为后退）称为后退）称为后退）称为后退EulerEuler公式公式公式公式 所谓局部载断误差是指以所谓局部载断误差是指以所谓局部载断误差是指以所谓局部载断误差是指以y yn n代替代替代替代替y y( (x xn n) )而得到而得到而得到而得到y y( (x xn n+1+1) )的近的近的近的近似值似值似值似值y yn n+1+1的误差（只估计这一步的误差）。
的误差（只估计这一步的误差）的误差（只估计这一步的误差）的误差（只估计这一步的误差）三、利用数值微分公式：利用两点公式三、利用数值微分公式：利用两点公式并且并且并且并且EulerEuler公式的公式的公式的公式的局部截断误差为局部截断误差为局部截断误差为局部截断误差为: :后退后退后退后退EulerEuler公式的公式的公式的公式的局部截断误差为局部截断误差为局部截断误差为局部截断误差为: : 10 W YEuler公式的推导（续公式的推导（续4））11 W YEuler公式的推导（续公式的推导（续5））四四四四、、、、利用数值积分公式：在利用数值积分公式：在利用数值积分公式：在利用数值积分公式：在[ [x xn n, ,x xn n+1+1] ]上对上对上对上对y y ( (x x)=)=f f ( (x x, ,y y( (x x)) )) 积分积分积分积分 对右端积分项采用不同的数值积分公式，便可得到各种对右端积分项采用不同的数值积分公式，便可得到各种对右端积分项采用不同的数值积分公式，便可得到各种对右端积分项采用不同的数值积分公式，便可得到各种不同的求解不同的求解不同的求解不同的求解dEdE初值问题的计算公式。
初值问题的计算公式初值问题的计算公式初值问题的计算公式 如，以矩形面积代替曲边梯形面积如，以矩形面积代替曲边梯形面积如，以矩形面积代替曲边梯形面积如，以矩形面积代替曲边梯形面积 1 1））））以左矩形面积代替曲边梯形面积如图以左矩形面积代替曲边梯形面积如图以左矩形面积代替曲边梯形面积如图以左矩形面积代替曲边梯形面积如图8-28-2，亦即以，亦即以，亦即以，亦即以 yf(x, y)xnxxn+1图8-2 12 W Yyf(x, y)xnxxn+1图8-3yf(x, y)xnxxn+1图8-43 3）以梯形公式（面积）代替曲边形如图）以梯形公式（面积）代替曲边形如图）以梯形公式（面积）代替曲边形如图）以梯形公式（面积）代替曲边形如图8-48-4则有则有则有则有 式（式（式（式（8-58-5）称为求）称为求）称为求）称为求dEdE初值问题的梯形公式，梯形公式看作初值问题的梯形公式，梯形公式看作初值问题的梯形公式，梯形公式看作初值问题的梯形公式，梯形公式看作是以（是以（是以（是以（x xn n, ,y yn n））））( (x xn n+1+1, ,y yn n+1+1) )构造的插值多项式代替被积函数构造的插值多项式代替被积函数构造的插值多项式代替被积函数构造的插值多项式代替被积函数得到的，而得到的，而得到的，而得到的，而EulerEuler公式公式公式公式则是以左端点函数值近似被积函数则是以左端点函数值近似被积函数则是以左端点函数值近似被积函数则是以左端点函数值近似被积函数而得到，还可以用多个点做插值多项式近似被积函数构造而得到，还可以用多个点做插值多项式近似被积函数构造而得到，还可以用多个点做插值多项式近似被积函数构造而得到，还可以用多个点做插值多项式近似被积函数构造另一些精度更高的解微分方程的数值公式，梯形公式比另一些精度更高的解微分方程的数值公式，梯形公式比另一些精度更高的解微分方程的数值公式，梯形公式比另一些精度更高的解微分方程的数值公式，梯形公式比EulerEuler公式公式公式公式更准确一些，误差更小。
更准确一些，误差更小更准确一些，误差更小更准确一些，误差更小 Euler公式的推导（续公式的推导（续6））2 2））））以右矩形面积代替曲边梯形：如图以右矩形面积代替曲边梯形：如图以右矩形面积代替曲边梯形：如图以右矩形面积代替曲边梯形：如图8-38-3亦即以亦即以亦即以亦即以 13 W YEuler公式注释公式注释注注注注1 1 ：：：：EulerEuler公式为显式，右矩形，梯形公式为隐式；公式为显式，右矩形，梯形公式为隐式；公式为显式，右矩形，梯形公式为隐式；公式为显式，右矩形，梯形公式为隐式； 注注注注2 2：：：：EulerEuler公式，梯形，右矩形公式为单步法，计算公式，梯形，右矩形公式为单步法，计算公式，梯形，右矩形公式为单步法，计算公式，梯形，右矩形公式为单步法，计算y yn n+1+1只用只用只用只用y yn n，而中点法公式为多步法（还可如上二所述，构造，而中点法公式为多步法（还可如上二所述，构造，而中点法公式为多步法（还可如上二所述，构造，而中点法公式为多步法（还可如上二所述，构造多步法）即必须已知多步法）即必须已知多步法）即必须已知多步法）即必须已知y yn n-1-1，，，，y yn n才才才才 能计算能计算能计算能计算y yn n+1+1，（求，（求，（求，（求y y0 0 , ,y y1 1不能不能不能不能用此公式。
用此公式用此公式用此公式y y0 0 , , y y1 1称为多步法的开始值，称为多步法的开始值，称为多步法的开始值，称为多步法的开始值，y y0 0给定，给定，给定，给定， 而而而而y y1 1必必必必须须须须由其它公式算出，然后才能用中点法）；由其它公式算出，然后才能用中点法）；由其它公式算出，然后才能用中点法）；由其它公式算出，然后才能用中点法）； 注注注注3 3：：：：前面已有前面已有前面已有前面已有EulerEuler法法法法 的局部截断误差：的局部截断误差：的局部截断误差：的局部截断误差： 后退后退后退后退EulerEuler法的局部截断误差：法的局部截断误差：法的局部截断误差：法的局部截断误差： 误差阶误差阶误差阶误差阶：：：：如果局部截断误差如果局部截断误差如果局部截断误差如果局部截断误差 则称方法为则称方法为则称方法为则称方法为P P阶的14 W Y 显然，步长显然，步长显然，步长显然，步长h h越小，阶数越小，阶数越小，阶数越小，阶数P P越高，局部截断误差越小，当越高，局部截断误差越小，当越高，局部截断误差越小，当越高，局部截断误差越小，当然计算精度越高；然计算精度越高；然计算精度越高；然计算精度越高； 注注注注4 4：：：：梯形法是几阶？梯形法精度比梯形法是几阶？梯形法精度比梯形法是几阶？梯形法精度比梯形法是几阶？梯形法精度比EulerEuler法高，阶数肯定法高，阶数肯定法高，阶数肯定法高，阶数肯定比比比比EulerEuler法高，其实我们可以利用数值积分公式的误差估法高，其实我们可以利用数值积分公式的误差估法高，其实我们可以利用数值积分公式的误差估法高，其实我们可以利用数值积分公式的误差估计式，因为我们是用梯形数值积计式，因为我们是用梯形数值积计式，因为我们是用梯形数值积计式，因为我们是用梯形数值积 分公式计算分公式计算分公式计算分公式计算 因此由积分中梯形公式的误差知此因此由积分中梯形公式的误差知此因此由积分中梯形公式的误差知此因此由积分中梯形公式的误差知此时的局部截断误差为：时的局部截断误差为：时的局部截断误差为：时的局部截断误差为：∴∴梯形法为梯形法为2阶方法阶方法 ！！ Euler Euler法，后退法，后退法，后退法，后退EulerEuler法为法为法为法为1 1阶方法阶方法阶方法阶方法，，，，而中点法为而中点法为而中点法为而中点法为2 2 阶阶阶阶，，，，Euler公式注释（续）公式注释（续）15 W Y关于关于Euler法的整体截断误差注释法的整体截断误差注释注注注注5 5：：：：关于关于关于关于EulerEuler法法法法的整体截断误差：的整体截断误差：的整体截断误差：的整体截断误差： EulerEuler方法的局部截断误差公式为：方法的局部截断误差公式为：方法的局部截断误差公式为：方法的局部截断误差公式为： 实际计算时，实际计算时，实际计算时，实际计算时，y yn n是是是是y y( (x xn n) )的近似值，因此，计算过程的近似值，因此，计算过程的近似值，因此，计算过程的近似值，因此，计算过程中除每步所产生的局部截断中除每步所产生的局部截断中除每步所产生的局部截断中除每步所产生的局部截断误差外，还有因前面的计算不准确而引起的误差。
在不考误差外，还有因前面的计算不准确而引起的误差在不考误差外，还有因前面的计算不准确而引起的误差在不考误差外，还有因前面的计算不准确而引起的误差在不考虑舍入误差的情况下，称虑舍入误差的情况下，称虑舍入误差的情况下，称虑舍入误差的情况下，称y y( (x xn n+1+1) )与与与与y yn n+1+1之差为整体截断误差之差为整体截断误差之差为整体截断误差之差为整体截断误差，记为：，记为：，记为：，记为： 下面讨论下面讨论下面讨论下面讨论EulerEuler方法的整体截断误差方法的整体截断误差方法的整体截断误差方法的整体截断误差 为简便起见，假定函数为简便起见，假定函数为简便起见，假定函数为简便起见，假定函数f f ( (x x, , y y) )充分光滑，问题（充分光滑，问题（充分光滑，问题（充分光滑，问题（8-18-1）解）解）解）解y y( (x x) )在在在在[ [a a, , b b] ]上二阶连续可微，于是由式（上二阶连续可微，于是由式（上二阶连续可微，于是由式（上二阶连续可微，于是由式（8-68-6），局部截），局部截），局部截），局部截断断断断误差有界，即存在误差有界，即存在误差有界，即存在误差有界，即存在MM>0>0，，，，使得对任意使得对任意使得对任意使得对任意x x [ [a a, ,b b] ]，都有，都有，都有，都有| |y y′( ′(x x)|≤)|≤MM，从而有：，从而有：，从而有：，从而有： （紧接下屏）（紧接下屏）（紧接下屏）（紧接下屏）16 W Y 式（式（式（式（8-78-7）表明，）表明，）表明，）表明，EulerEuler方法的整体截断误差与方法的整体截断误差与方法的整体截断误差与方法的整体截断误差与h h同阶，当同阶，当同阶，当同阶，当h h0 0时，时，时，时，e en n0 0。
关于关于Euler法的整体截断误差注释（续）法的整体截断误差注释（续） 反复递推得：反复递推得：反复递推得：反复递推得：17 W Y结结 论论 对于实际问题来说，由于对于实际问题来说，由于对于实际问题来说，由于对于实际问题来说，由于L L，，，，M M 难以估计，因此（难以估计，因此（难以估计，因此（难以估计，因此（8-8-6 6））））很难应用，而且上述推导过程中一再放大了误差上限，这很难应用，而且上述推导过程中一再放大了误差上限，这很难应用，而且上述推导过程中一再放大了误差上限，这很难应用，而且上述推导过程中一再放大了误差上限，这样的估计往往也很保守，远远大于实际的误差，但是，从样的估计往往也很保守，远远大于实际的误差，但是，从样的估计往往也很保守，远远大于实际的误差，但是，从样的估计往往也很保守，远远大于实际的误差，但是，从估计式（估计式（估计式（估计式（8-78-7）却可以得到下面很有用处的结论却可以得到下面很有用处的结论却可以得到下面很有用处的结论却可以得到下面很有用处的结论。
1 1））））当当当当h h0 0时，时，时，时，e en n0 0即，即，即，即， 亦即数值解亦即数值解亦即数值解亦即数值解y yn n，，，，一致收敛于初值问题（一致收敛于初值问题（一致收敛于初值问题（一致收敛于初值问题（8-18-1）的真解）的真解）的真解）的真解y y( (x xn n) )，，，， 并且，并且，并且，并且，EulerEuler法的整体截断误差的阶为法的整体截断误差的阶为法的整体截断误差的阶为法的整体截断误差的阶为O O( (h h) )与与与与h h同阶，比局同阶，比局同阶，比局同阶，比局部截断误差低一阶部截断误差低一阶部截断误差低一阶部截断误差低一阶 2 2））））舍入误差舍入误差舍入误差舍入误差 局部截断误差局部截断误差局部截断误差局部截断误差 对实际计算结果有影响，并且随对实际计算结果有影响，并且随对实际计算结果有影响，并且随对实际计算结果有影响，并且随h h减少减少减少减少 而减少或增大而减少或增大而减少或增大而减少或增大18 W Y3））计算结果与解法的阶数计算结果与解法的阶数p，真解的导数，真解的导数y (p+1)有有关，关，p越大，越大，h p+1越小，越小， |y(p+1)(ξ)|的上限越大，的上限越大，M 也越大，因此为保证精度当然应选阶数也越大，因此为保证精度当然应选阶数p较高的较高的方法。
但如果方法但如果M 很大，当很大，当f (x,y)是分段连续的函是分段连续的函数时，则应采用低阶的方法如用数时，则应采用低阶的方法如用Euler法结结 论（续）论（续） 4））计算结果还与开始值的精度有关，为使这种计算结果还与开始值的精度有关，为使这种误差的影响不致于超过局部误差的影响不致于超过局部 截断误差，对多步截断误差，对多步法，应采用跟多步法同阶的方法计算开始值法，应采用跟多步法同阶的方法计算开始值 19 W Y1.2 改进的改进的Euler法法 梯形公式为二阶方法，但却是隐式格式，即若利用梯梯形公式为二阶方法，但却是隐式格式，即若利用梯梯形公式为二阶方法，但却是隐式格式，即若利用梯梯形公式为二阶方法，但却是隐式格式，即若利用梯形公式求形公式求形公式求形公式求y yn n+1+1，就要求解方程（，就要求解方程（，就要求解方程（，就要求解方程（8-58-5）式，计算量较大，通）式，计算量较大，通）式，计算量较大，通）式，计算量较大，通常在实际计算时，将常在实际计算时，将常在实际计算时，将常在实际计算时，将EulerEuler法法法法与与与与梯形公式梯形公式梯形公式梯形公式合起来使用，即合起来使用，即合起来使用，即合起来使用，即先使用先使用先使用先使用EulerEuler公式公式公式公式, ,由（由（由（由（x xn n, ,y yn n）算出）算出）算出）算出y yn n+1+1，记为，记为，记为，记为y yn n+1+1(0)(0)，称为，称为，称为，称为预测值预测值预测值预测值，然后用，然后用，然后用，然后用梯形公式梯形公式梯形公式梯形公式去提高精度，将去提高精度，将去提高精度，将去提高精度，将y yn n+1+1(0)(0) 校正为较校正为较校正为较校正为较准确的值：准确的值：准确的值：准确的值：由于函数由于函数由于函数由于函数f f ( (x x, ,y y) )满足满足满足满足LipschitzLipschitz条件，容易得出条件，容易得出条件，容易得出条件，容易得出: :20 W Y改进的改进的Euler法法(续）续）21 W Y预测预测──校正型公式校正型公式 实际经验表明，式（实际经验表明，式（实际经验表明，式（实际经验表明，式（8-88-8）的迭代效果主要体现在第一次，）的迭代效果主要体现在第一次，）的迭代效果主要体现在第一次，）的迭代效果主要体现在第一次，由此构成如下的预测由此构成如下的预测由此构成如下的预测由此构成如下的预测────校正型公式校正型公式校正型公式校正型公式: :此式称为改进的此式称为改进的此式称为改进的此式称为改进的EulerEuler公式，为上机计算编程方便，常将式公式，为上机计算编程方便，常将式公式，为上机计算编程方便，常将式公式，为上机计算编程方便，常将式（（（（8-98-9）改写为）改写为）改写为）改写为: :下面分析改进的下面分析改进的下面分析改进的下面分析改进的EulerEuler公式的局部截断误差：公式的局部截断误差：公式的局部截断误差：公式的局部截断误差： 22 W Y改进的改进的Euler公式的局部截断误差分析公式的局部截断误差分析假定假定假定假定y yn n = = y y( (x xn n) )，，，，y y( (x xn n+1+1) ) 的的的的TaylorTaylor展式为展式为展式为展式为: : 对于改进的对于改进的对于改进的对于改进的EulerEuler公式，由于公式，由于公式，由于公式，由于 这说明改进的这说明改进的这说明改进的这说明改进的EulerEuler法的局部法的局部法的局部法的局部截断误差为截断误差为截断误差为截断误差为O O( (h h3 3) )，比，比，比，比EulerEuler公式高一阶，是二阶方法。
公式高一阶，是二阶方法公式高一阶，是二阶方法公式高一阶，是二阶方法 23 W Y改进的改进的Euler公式举例公式举例例例1 这些结果在表这些结果在表这些结果在表这些结果在表8-18-1中，可见计算结果的精度，中，可见计算结果的精度，中，可见计算结果的精度，中，可见计算结果的精度，EulerEuler法法法法与与与与后退后退后退后退EulerEuler法法法法差不多，与准确值相比较差不多，与准确值相比较差不多，与准确值相比较差不多，与准确值相比较EulerEuler法法法法偏小，而偏小，而偏小，而偏小，而后后后后退退退退EulerEuler法法法法偏大；偏大；偏大；偏大；中点法中点法中点法中点法与与与与梯形法梯形法梯形法梯形法精度同为精度同为精度同为精度同为2 2阶，但梯形法阶，但梯形法阶，但梯形法阶，但梯形法更好一些，这跟它们局部截断误差的符号，阶数和系更好一些，这跟它们局部截断误差的符号，阶数和系更好一些，这跟它们局部截断误差的符号，阶数和系更好一些，这跟它们局部截断误差的符号，阶数和系数的大小是完全一致的数的大小是完全一致的数的大小是完全一致的。
数的大小是完全一致的 表见下屏：表见下屏：表见下屏：表见下屏：24 W Y表格表格8-1表表表表8-18-1 y y = = y y, , y y(0)=1 (0)=1 的数值解的数值解的数值解的数值解( (h h=0.1)=0.1) x x精精精精 确确确确 解解解解欧拉法欧拉法欧拉法欧拉法后退欧拉后退欧拉后退欧拉后退欧拉中点法中点法中点法中点法梯形法梯形法梯形法梯形法.1 .1.904837.904837.900000.900000.909091.909091.900000.900000.904762.904762.2 .2.808731.808731.810000.810000.826446.826446.820000.820000.818594.818594. .3 3.740818.740818.729000.729000.751315.751315.736000.736000.740633.740633.4 .4.670320.670320.656100.656100.683013.683013.627800.627800.670096.670096.5 .5.060531.060531.590490.590490.620921.620921.601440.601440.606278.606278.6 .6.548812.548812.531441.531441.654474.654474.552512.552512.548537.548537.7 .7.496585.496585. 478298. 478298.5131458.5131458.490938.490938.496295.496295.8 .8.449329.449329.430467.430467.466507.466507.454324.454324.449029.449029.9 .9.406570.406570.387421.387421.424098.424098.400073.400073.406264.4062641 1.367879.367879.348679.348679.385543.385543.374310.374310.367573.36757325 W Y表格表格表格表格8-28-2 而而而而表表表表8-28-2是分别取了不同的是分别取了不同的是分别取了不同的是分别取了不同的h h=0.1 ,=0.1 ,h h=0.01=0.01，，，，h h=0.001,=0.001,h h=0.0001=0.0001，还是利用这些公式，经过若干步的计算（，还是利用这些公式，经过若干步的计算（，还是利用这些公式，经过若干步的计算（，还是利用这些公式，经过若干步的计算（h h越越越越小，计算量越大）算到小，计算量越大）算到小，计算量越大）算到小，计算量越大）算到y y(1)(1)的近似值，的近似值，的近似值，的近似值，可见可见可见可见：：：： 随着随着随着随着h h的减小，的减小，的减小，的减小，y y(1)(1)的近似值的精度在提高，的近似值的精度在提高，的近似值的精度在提高，的近似值的精度在提高，0.010.01比比比比0.0010.001差，即差，即差，即差，即0.0010.001比比比比0.010.01时的时的时的时的y y(1)(1)准确。
准确 Y ′=-y, y(0)=1的解的解y(1)的近似值的近似值(y(1)=0.367879)h欧拉法欧拉法后退欧拉法后退欧拉法中点法中点法梯形法梯形法0.1.348678.385543.374310.3675730.01.366033.369711.367944.3678770.001.367700.368052.367879.3678760.0001.367800.367800.367881.368020（紧接下屏）（紧接下屏）（紧接下屏）（紧接下屏）26 W Y表表8-2计算结果说明（续）计算结果说明（续） 但但h太小，到太小，到h=0.0001时却又变得误差大了，这时却又变得误差大了，这与前面所说与前面所说h越小，越小，p阶越高，应该局部截断误差越阶越高，应该局部截断误差越小，因而计算精度更高矛盾了，为什么会产生这种小，因而计算精度更高矛盾了，为什么会产生这种情况呢？情况呢？ 这是由于这是由于h太小而引起计算量大因而造成了舍入太小而引起计算量大因而造成了舍入误差和截断误差的积累，这种情况由于初值问题不误差和截断误差的积累，这种情况由于初值问题不同可能会影响更大，偏离更严重，同可能会影响更大，偏离更严重，如下面的例如下面的例2 。
这种问题实际上是稳定性问题，我们将会讨论方法这种问题实际上是稳定性问题，我们将会讨论方法的稳定性，由此得出对的稳定性，由此得出对h有一定的要求的稳定性制有一定的要求的稳定性制区域 27 W Y例例2 求解初值问题求解初值问题求解初值问题求解初值问题y y′=-20′=-20y y，，，，y y(0)=1(0)=1，计算，计算，计算，计算y y(1)(1)的近似值的近似值的近似值的近似值 解解解解：：：：类似于例类似于例类似于例类似于例1 1，用欧拉法、后退欧拉法、中点法、梯形，用欧拉法、后退欧拉法、中点法、梯形，用欧拉法、后退欧拉法、中点法、梯形，用欧拉法、后退欧拉法、中点法、梯形法求解，得到如下法求解，得到如下法求解，得到如下法求解，得到如下表表表表8.38.3 表表表表8.38.3 y y = = 20,20,y y(0)= (0)= 1 1的解的解的解的解y(1)y(1)的近似值的近似值的近似值的近似值( (y y(1)=0.206116(1)=0.206116E E 8) 8) h h欧拉法欧拉法欧拉法欧拉法后退欧拉法后退欧拉法后退欧拉法后退欧拉法中点法中点法中点法中点法梯形法梯形法梯形法梯形法.1 .11 1.169351E-4.169351E-4.514229E+6.514229E+60 0.01.01.203704E-9.203704E-9.120746E-7.120746E-7.413244E+7.413244E+7.192743E-8.192743E-8.001.001.168300E-8.168300E-8.251090E-8.251090E-8.484136E+5.484136E+5.205979E-8.205979E-8.0001.0001.202081E-8.202081E-8.210331E-8.210331E-8.475130+3.475130+3.206103E-8.206103E-8 由由由由表表表表8.38.3可见，尽管可见，尽管可见，尽管可见，尽管中点法中点法中点法中点法的阶数与的阶数与的阶数与的阶数与梯形法梯形法梯形法梯形法相同，比相同，比相同，比相同，比欧拉法欧拉法欧拉法欧拉法和和和和后退欧拉法后退欧拉法后退欧拉法后退欧拉法的阶数高，计算结果的精度却很糟的阶数高，计算结果的精度却很糟的阶数高，计算结果的精度却很糟的阶数高，计算结果的精度却很糟糕。
此外，尽管糕此外，尽管糕此外，尽管糕此外，尽管欧拉法欧拉法欧拉法欧拉法与与与与后退欧拉法后退欧拉法后退欧拉法后退欧拉法的阶数相同，但的阶数相同，但的阶数相同，但的阶数相同，但欧欧欧欧拉法拉法拉法拉法计算结果的精度，当计算结果的精度，当计算结果的精度，当计算结果的精度，当h h=0.1=0.1时却比时却比时却比时却比后退欧拉法后退欧拉法后退欧拉法后退欧拉法差 28 W Y§2 龙格龙格-库塔（库塔（Runge-kutta）方法）方法 2.1 2.1 龙格龙格龙格龙格- -库塔方法的基本思想库塔方法的基本思想库塔方法的基本思想库塔方法的基本思想：：：： 因此只要对平均斜率因此只要对平均斜率因此只要对平均斜率因此只要对平均斜率k k* *提供一种算法，由（提供一种算法，由（提供一种算法，由（提供一种算法，由（8-118-11）式）式）式）式便相应地得到一种微分方程的数值计算公式便相应地得到一种微分方程的数值计算公式便相应地得到一种微分方程的数值计算公式便相应地得到一种微分方程的数值计算公式紧接下屏）（紧接下屏）（紧接下屏）（紧接下屏）29 W Y龙格龙格-库塔方法的基本思想库塔方法的基本思想(续）续） 改进欧拉公式比欧拉公式精度高的原因，也就在于确改进欧拉公式比欧拉公式精度高的原因，也就在于确改进欧拉公式比欧拉公式精度高的原因，也就在于确改进欧拉公式比欧拉公式精度高的原因，也就在于确定平均斜率时，多取了一个点的斜率值。
因此它启发我们，定平均斜率时，多取了一个点的斜率值因此它启发我们，定平均斜率时，多取了一个点的斜率值因此它启发我们，定平均斜率时，多取了一个点的斜率值因此它启发我们，如果设法在如果设法在如果设法在如果设法在[ [x xi i, ,x xi i+1+1] ]上多预报几个点的斜率值，然后将它上多预报几个点的斜率值，然后将它上多预报几个点的斜率值，然后将它上多预报几个点的斜率值，然后将它们们们们加权平均作为加权平均作为加权平均作为加权平均作为k*k*的近似值，则有可能构造出更高精度的计的近似值，则有可能构造出更高精度的计的近似值，则有可能构造出更高精度的计的近似值，则有可能构造出更高精度的计算公式，这是算公式，这是算公式，这是算公式，这是龙格龙格龙格龙格- -库塔库塔库塔库塔方法的基本思想方法的基本思想方法的基本思想方法的基本思想 用这个观点来研究欧拉公式与改进欧拉公式，可以发现用这个观点来研究欧拉公式与改进欧拉公式，可以发现用这个观点来研究欧拉公式与改进欧拉公式，可以发现用这个观点来研究欧拉公式与改进欧拉公式，可以发现欧拉公式由于仅取欧拉公式由于仅取欧拉公式由于仅取欧拉公式由于仅取x xn n一个点的斜率值一个点的斜率值一个点的斜率值一个点的斜率值f f ( (x xn n, ,y yn n) )作为平均斜率作为平均斜率作为平均斜率作为平均斜率k k* * 的近似值，因此精度很低。
而改进欧拉公式（的近似值，因此精度很低而改进欧拉公式（的近似值，因此精度很低而改进欧拉公式（的近似值，因此精度很低而改进欧拉公式（8-108-10）却）却）却）却是利用了是利用了是利用了是利用了x xn n与与与与x xn n-1-1两个点的斜率值两个点的斜率值两个点的斜率值两个点的斜率值k k1 1 = = f f ( (x xn n, ,y yn n) )与与与与k k2 2= =f f ( (x xn n+1+1, ,y yn n+ +hkhk1 1) )取算术平均作为平均斜率取算术平均作为平均斜率取算术平均作为平均斜率取算术平均作为平均斜率k k* *的近似值的近似值的近似值的近似值 其中其中其中其中k k2 2是通过已知信息是通过已知信息是通过已知信息是通过已知信息y yn n利利利利用欧拉公式求得的用欧拉公式求得的用欧拉公式求得的用欧拉公式求得的30 W Y2.2 二阶龙格二阶龙格-库塔公式库塔公式 首先推广改进欧拉公式，首先推广改进欧拉公式，首先推广改进欧拉公式，首先推广改进欧拉公式，考察区间考察区间考察区间考察区间[ [x xn n, ,x xn n+1+1] ]内任一点：内任一点：内任一点：内任一点： 我们希望用我们希望用我们希望用我们希望用x xn n和和和和x x n n+1+1两个点的斜率值两个点的斜率值两个点的斜率值两个点的斜率值k k1 1和和和和k k2 2加权平均作为加权平均作为加权平均作为加权平均作为 平均斜率平均斜率平均斜率平均斜率k k* *的近似值的近似值的近似值的近似值 ：：：：其中其中其中其中c c1 1，，，，c c2 2为待定常数，同为待定常数，同为待定常数，同为待定常数，同改进欧拉公式一样，这里仍取：改进欧拉公式一样，这里仍取：改进欧拉公式一样，这里仍取：改进欧拉公式一样，这里仍取：问题在于怎样预测问题在于怎样预测问题在于怎样预测问题在于怎样预测x xn n+ +l l处的斜率值处的斜率值处的斜率值处的斜率值k k2 2。
仿照改进欧拉公式，先用仿照改进欧拉公式，先用仿照改进欧拉公式，先用仿照改进欧拉公式，先用 欧拉公式提供欧拉公式提供欧拉公式提供欧拉公式提供y y ( (x xn n+ +l l) )的预测值的预测值的预测值的预测值 然后再用预测值然后再用预测值然后再用预测值然后再用预测值y yn n+ +l l通过计算通过计算通过计算通过计算f f 产生斜率值产生斜率值产生斜率值产生斜率值k k2 2= =f f ( (x xn n+ +l l , ,y yn n+ +l l) )，，，， 这样设计出的计算公式具有形式这样设计出的计算公式具有形式这样设计出的计算公式具有形式这样设计出的计算公式具有形式 ：：：：31 W Y二阶龙格二阶龙格-库塔公式（续库塔公式（续1）） 公式（公式（公式（公式（8-128-12）中含有三个待定参数）中含有三个待定参数）中含有三个待定参数）中含有三个待定参数c c1 1, ,c c2 2和和和和l l，我们希望适，我们希望适，我们希望适，我们希望适当选取这些参数值，使得公式（当选取这些参数值，使得公式（当选取这些参数值，使得公式（当选取这些参数值，使得公式（8-128-12）具有二阶精度，亦）具有二阶精度，亦）具有二阶精度，亦）具有二阶精度，亦即使：即使：即使：即使：现在仍假定现在仍假定现在仍假定现在仍假定y yn n= =y y( (x xn n) )，即，即，即，即y yn n是准确的，将是准确的，将是准确的，将是准确的，将y y( (x xn n+1+1) )与与与与y yn n+1+1都在都在都在都在x xi i处作泰勒展开：处作泰勒展开：处作泰勒展开：处作泰勒展开：32 W Y二阶龙格二阶龙格-库塔公式（续库塔公式（续2）） 代入（代入（代入（代入（8-128-12）式，得：）式，得：）式，得：）式，得： 比较（比较（比较（比较（8-138-13）与（）与（）与（）与（8-148-14）两式，要使公式具有二阶精度，）两式，要使公式具有二阶精度，）两式，要使公式具有二阶精度，）两式，要使公式具有二阶精度，只有满足下列条件：只有满足下列条件：只有满足下列条件：只有满足下列条件： 这里一共有三个待定参数，但只需满足两个条件，因这里一共有三个待定参数，但只需满足两个条件，因这里一共有三个待定参数，但只需满足两个条件，因这里一共有三个待定参数，但只需满足两个条件，因此有一个自由度，于是满足条件（此有一个自由度，于是满足条件（此有一个自由度，于是满足条件（此有一个自由度，于是满足条件（8-158-15）的参数不止一组，）的参数不止一组，）的参数不止一组，）的参数不止一组，而是一族，相应的公式（而是一族，相应的公式（而是一族，相应的公式（而是一族，相应的公式（8-128-12）也有一族，这些公式统称）也有一族，这些公式统称）也有一族，这些公式统称）也有一族，这些公式统称为为为为二二二二阶龙格阶龙格-库塔公式库塔公式，简称，简称，简称，简称二阶二阶R-K公式公式 。
33 W Y 特别，当特别，当特别，当特别，当l l=1=1即即即即x xn+ln+l=x=xn+n+1 1时，时，时，时， c c1 1= =c c2 2=1/2=1/2，二阶，二阶，二阶，二阶R R- -K K公式就公式就公式就公式就 是是是是改进欧拉公式改进欧拉公式改进欧拉公式改进欧拉公式 如果取如果取如果取如果取l l=1/2=1/2，则，则，则，则c c1 1=0,=0,c c2 2=1=1，，，， 这时二阶这时二阶这时二阶这时二阶R R- -K K公式称为公式称为公式称为公式称为变形变形变形变形的欧拉公式的欧拉公式的欧拉公式的欧拉公式，其形式见左边：，其形式见左边：，其形式见左边：，其形式见左边： 从表面上看，变形的欧拉公式仅含一个斜率值从表面上看，变形的欧拉公式仅含一个斜率值k2，但，但k2是是通过通过k1计算出来的，因此每完成一步，计算出来的，因此每完成一步，仍然需要两次计算函数仍然需要两次计算函数f 的值，工作量和改进欧拉的值，工作量和改进欧拉公式相同。
公式相同 二阶龙格二阶龙格-库塔公式（续库塔公式（续3））34 W Y构造二阶构造二阶R-K公式的主要步骤公式的主要步骤 综上所述，构造二阶综上所述，构造二阶R-K公式主要由以下几步产生：公式主要由以下几步产生： 1)在区间在区间[xn,xn+1]上取二点，预报相应点上取二点，预报相应点 的斜率值；的斜率值；2) 对此两斜率值加权平均作为平均斜率对此两斜率值加权平均作为平均斜率 值的近似值；值的近似值；3) 将将yn+1与与y(xn+1)都在都在xi处作泰勒展开，处作泰勒展开， 为使公式达到二阶精度，比较相应系为使公式达到二阶精度，比较相应系 数，建立有关参数所应满足的方程组；数，建立有关参数所应满足的方程组；4) 解此方程组得一族二阶解此方程组得一族二阶R-K公式35 W Y2.3 高阶高阶R-K公式公式 为了进一步提高精度，在为了进一步提高精度，在为了进一步提高精度，在为了进一步提高精度，在[ [x xn n, ,x xn n+1+1] ]上除上除上除上除x xn n和和和和x xn n+ +l l外再增外再增外再增外再增加加加加一点一点一点一点x xn+mn+m= =x xn n+ +mhmh( (l l m m 1) 1)，并用，并用，并用，并用x xn n, ,x xn n+l+l, ,x xn n+ +mm三处的斜率三处的斜率三处的斜率三处的斜率值值值值k k1 1, ,k k2 2, ,k k3 3加权平均作为加权平均作为加权平均作为加权平均作为k k* *的近似值，这时计算公式为：的近似值，这时计算公式为：的近似值，这时计算公式为：的近似值，这时计算公式为：其中其中其中其中k k1 1, ,k k2 2仍用（仍用（仍用（仍用（8-128-12）式所取的形式。
式所取的形式式所取的形式式所取的形式 为了预测为了预测为了预测为了预测x xn n+ +mm处的斜率值处的斜率值处的斜率值处的斜率值k k3 3，要定出，要定出，要定出，要定出x xn n+ +mm处所对应的处所对应的处所对应的处所对应的y yn n+ +mm，可以看作在区间，可以看作在区间，可以看作在区间，可以看作在区间[ [x xn n, ,x xn n+ +mm] ]上使用二阶上使用二阶上使用二阶上使用二阶R-KR-K公式，从而公式，从而公式，从而公式，从而得到得到得到得到y y( (x xn n+ +mm) )的预测值：的预测值：的预测值：的预测值： 于是，再通过计算函数值于是，再通过计算函数值于是，再通过计算函数值于是，再通过计算函数值f f 得到得到得到得到: :36 W Y高阶高阶R-K公式（续公式（续1））这样设计出的计算公式具有形式这样设计出的计算公式具有形式:运用泰勒展开方法，适当选择参数运用泰勒展开方法，适当选择参数c1,c2,c3,l,m,及及b1,b2使上述公式具有三阶精度采用与使上述公式具有三阶精度采用与(8-12)类似类似的处理方法，得到这些参数需要满足的条件：的处理方法，得到这些参数需要满足的条件： (见下屏见下屏)37 W Y高阶高阶R-K公式（续公式（续2））满足满足满足满足5 5个条件的一族公式（个条件的一族公式（个条件的一族公式（个条件的一族公式（8-188-18）））），统称为三阶，统称为三阶，统称为三阶，统称为三阶R-KR-K公式，其中常公式，其中常公式，其中常公式，其中常见的是库塔公式：见的是库塔公式：见的是库塔公式：见的是库塔公式： c c1 1=1/6 , c=1/6 , c2 2=4/6 , c=4/6 , c3 3=1/6 , =1/6 , l=1/2 , m=1,b l=1/2 , m=1,b1 1=-1 , b=-1 , b2 2=2=238 W Y 若需再将精度提高至四阶，用上述类似处理方法，只若需再将精度提高至四阶，用上述类似处理方法，只若需再将精度提高至四阶，用上述类似处理方法，只若需再将精度提高至四阶，用上述类似处理方法，只是必须在是必须在是必须在是必须在[ [x xi i, ,x xi i+1+1] ]上用四个点处的斜率加权平均作为上用四个点处的斜率加权平均作为上用四个点处的斜率加权平均作为上用四个点处的斜率加权平均作为k k* *的的的的近近近近似值，构成一族四阶龙格似值，构成一族四阶龙格似值，构成一族四阶龙格似值，构成一族四阶龙格- -库塔公式，由于推导复杂，这库塔公式，由于推导复杂，这库塔公式，由于推导复杂，这库塔公式，由于推导复杂，这里从略，只将常用的公式介绍如下：里从略，只将常用的公式介绍如下：里从略，只将常用的公式介绍如下：里从略，只将常用的公式介绍如下： 四阶公式中常用的还四阶公式中常用的还有下面的有下面的Gill公式：公式：高阶高阶R-K公式（续公式（续3））一般称它为标准四阶一般称它为标准四阶R-K公式，其局部截断误公式，其局部截断误差是差是O(h5) 。
39 W YGill 公式公式 在计算机上它具有节省存贮单元和控制舍入在计算机上它具有节省存贮单元和控制舍入误差增长的优点误差增长的优点 40 W Y初值问题两种方法举例初值问题两种方法举例用改进的用改进的用改进的用改进的EulerEuler法和四阶标准法和四阶标准法和四阶标准法和四阶标准R-KR-K公式解初值公式解初值公式解初值公式解初值问题：问题：问题：问题：例例3[ [解解解解] ] 用改进的用改进的用改进的用改进的EulerEuler公式，取公式，取公式，取公式，取h h=0.1=0.1计算公式为：计算公式为：计算公式为：计算公式为： 部分计算结果见表部分计算结果见表部分计算结果见表部分计算结果见表8-48-4：：：： 表表表表8-48-4x xn n改进的欧改进的欧改进的欧改进的欧拉法拉法拉法拉法准准准准 确确确确 解解解解y y( (x xn n) )误差误差误差误差x xn n改进的欧改进的欧改进的欧改进的欧拉法拉法拉法拉法准准准准 确确确确 解解解解y y( (x xn n) )误差误差误差误差0 01. 1.1 10 00.60.61.4859561.4859561.4832401.4832400.0027160.0027160.20.21.1840961.1840961.1832161.1832160.0000880.0000880.80.81.6164761.6164761.6124521.6124520.0040240.0040240.40.41.3433601.3433601.3416411.3416410.0017190.0017191.01.01.7378691.7378691.7320511.7320510.0058180.00581841 W Y对四阶标准对四阶标准R-K法，法，取取h=0.2，计算公，计算公式如下：式如下：例例3（续（续1））42 W Y例例 3 （续（续2））表表表表8-58-5 x xn n四阶四阶四阶四阶R-KR-K法法法法误差误差误差误差x xn n四阶四阶四阶四阶R-KR-K法法法法误差误差误差误差0 01 11 10.60.61.4832811.4832810.0000410.0000410.20.21.1832291.1832290.0000130.0000130.80.81.6125131.6125130.0000610.0000610.40.41.3416671.3416670.000260.000261.01.01.7321401.7321400.0000890.000089 计算结果见计算结果见计算结果见计算结果见计算结果见计算结果见表表表表表表8-58-58-5，与准确解比较，至少有四位有效，与准确解比较，至少有四位有效，与准确解比较，至少有四位有效，与准确解比较，至少有四位有效，与准确解比较，至少有四位有效，与准确解比较，至少有四位有效数字。
而取步长数字而取步长数字而取步长数字而取步长数字而取步长数字而取步长h h h=0.1=0.1=0.1时，用改进的欧拉公式计算，最多时，用改进的欧拉公式计算，最多时，用改进的欧拉公式计算，最多时，用改进的欧拉公式计算，最多时，用改进的欧拉公式计算，最多时，用改进的欧拉公式计算，最多只有三位有效数字只有三位有效数字只有三位有效数字只有三位有效数字只有三位有效数字只有三位有效数字（见表（见表（见表（见表（见表（见表8-48-48-4）））））），虽然改进的欧拉公式每，虽然改进的欧拉公式每，虽然改进的欧拉公式每，虽然改进的欧拉公式每，虽然改进的欧拉公式每，虽然改进的欧拉公式每前进一步只要计算两次前进一步只要计算两次前进一步只要计算两次前进一步只要计算两次前进一步只要计算两次前进一步只要计算两次f f f 值，而四阶值，而四阶值，而四阶值，而四阶值，而四阶值，而四阶R R R- - -K K K法要计算法要计算法要计算法要计算法要计算法要计算4 4 4次次次次次次f f f 值，值，值，值，值，值，但由于后者步长比前者大，所耗费的计算总量还是差不多但由于后者步长比前者大，所耗费的计算总量还是差不多但由于后者步长比前者大，所耗费的计算总量还是差不多但由于后者步长比前者大，所耗费的计算总量还是差不多但由于后者步长比前者大，所耗费的计算总量还是差不多但由于后者步长比前者大，所耗费的计算总量还是差不多的。
这个例子又一次显示了选择算法的重要意义这个例子又一次显示了选择算法的重要意义这个例子又一次显示了选择算法的重要意义这个例子又一次显示了选择算法的重要意义这个例子又一次显示了选择算法的重要意义这个例子又一次显示了选择算法的重要意义 43 W Y表表8-6计算计算f 的次数与精度阶数的关系的次数与精度阶数的关系 每每每每 步步步步 计计计计 算算算算f f 的的的的 次次次次 数数数数2 23 34 45 56 67 78 89 9精度的阶精度的阶精度的阶精度的阶数数数数2 23 34 44 45 56 66 67 7表表表表8-6 8-6 的说明的说明的说明的说明 一般，高精度的方法要求解有较好的光滑性，解的光滑一般，高精度的方法要求解有较好的光滑性，解的光滑一般，高精度的方法要求解有较好的光滑性，解的光滑一般，高精度的方法要求解有较好的光滑性，解的光滑性是由性是由性是由性是由f (x,y)f (x,y)的可微性决定如果解的光滑性差，则用四的可微性决定如果解的光滑性差，则用四的可微性决定如果解的光滑性差，则用四的可微性决定。
如果解的光滑性差，则用四阶阶阶阶R-KR-K法求得的数值解反而不如用改进的欧拉法好因此，法求得的数值解反而不如用改进的欧拉法好因此，法求得的数值解反而不如用改进的欧拉法好因此，法求得的数值解反而不如用改进的欧拉法好因此，一定要针对具体问题的特点来选择求解方法一定要针对具体问题的特点来选择求解方法一定要针对具体问题的特点来选择求解方法一定要针对具体问题的特点来选择求解方法例例 3 （续（续3））44 W Y表表8-6 的说明（续）的说明（续） 龙格龙格-库塔方法的推导是在用泰勒展开方法的库塔方法的推导是在用泰勒展开方法的基础上进行的，因而它要求所求的解具有较好的基础上进行的，因而它要求所求的解具有较好的光滑性质假若解的光滑性差，那么使用四阶光滑性质假若解的光滑性差，那么使用四阶R-K公式求得的数值解，其精度可能反而不如改进公式求得的数值解，其精度可能反而不如改进欧拉公式在实际计算时，应当针对问题的具体欧拉公式在实际计算时，应当针对问题的具体特点，选择合适的算法特点，选择合适的算法 从理论上讲，可以构造任意高阶的龙格从理论上讲，可以构造任意高阶的龙格从理论上讲，可以构造任意高阶的龙格从理论上讲，可以构造任意高阶的龙格- -库塔公式，库塔公式，库塔公式，库塔公式，但只要注意到精度的阶数与计算函数值但只要注意到精度的阶数与计算函数值但只要注意到精度的阶数与计算函数值但只要注意到精度的阶数与计算函数值f(x,y)f(x,y)的次数之的次数之的次数之的次数之间的关系不是等量增加的（见表间的关系不是等量增加的（见表间的关系不是等量增加的（见表间的关系不是等量增加的（见表8-68-6），就可得出再提），就可得出再提），就可得出再提），就可得出再提高公式阶数已没有多大意思了。
因此更说明高公式阶数已没有多大意思了因此更说明高公式阶数已没有多大意思了因此更说明高公式阶数已没有多大意思了因此更说明四阶四阶四阶四阶R-KR-K公公公公式式式式是兼顾了精度及计算量的较理想的计算公式是兼顾了精度及计算量的较理想的计算公式是兼顾了精度及计算量的较理想的计算公式是兼顾了精度及计算量的较理想的计算公式 需要指出的是：需要指出的是：45 W Y2.4 变步长变步长R-K法法 当用数值方法解微分方程时，单从一步来看，步长越当用数值方法解微分方程时，单从一步来看，步长越当用数值方法解微分方程时，单从一步来看，步长越当用数值方法解微分方程时，单从一步来看，步长越小，局部截断误差就越小但从整体来看，步长越小，在小，局部截断误差就越小但从整体来看，步长越小，在小，局部截断误差就越小但从整体来看，步长越小，在小，局部截断误差就越小但从整体来看，步长越小，在求解范围内所要完成的步数就会越多这一方面增加了计求解范围内所要完成的步数就会越多这一方面增加了计求解范围内所要完成的步数就会越多这一方面增加了计求解范围内所要完成的步数就会越多。
这一方面增加了计算量；另一方面也增大了舍入误差的累积，因此有个如何算量；另一方面也增大了舍入误差的累积，因此有个如何算量；另一方面也增大了舍入误差的累积，因此有个如何算量；另一方面也增大了舍入误差的累积，因此有个如何合理选择步长的问题在解决这个问题时，需要考虑两个合理选择步长的问题在解决这个问题时，需要考虑两个合理选择步长的问题在解决这个问题时，需要考虑两个合理选择步长的问题在解决这个问题时，需要考虑两个问题：问题：问题：问题： ①①①① 怎样衡量和检验计算结果的精度？怎样衡量和检验计算结果的精度？怎样衡量和检验计算结果的精度？怎样衡量和检验计算结果的精度？ ②②②②如何依据获得的精度信息及时处理步长如何依据获得的精度信息及时处理步长如何依据获得的精度信息及时处理步长如何依据获得的精度信息及时处理步长 以四阶标准以四阶标准以四阶标准以四阶标准R-KR-K公式为例，从节点公式为例，从节点公式为例，从节点公式为例，从节点x xn n出发，先以步长出发，先以步长出发，先以步长出发，先以步长h h求得求得求得求得x xn n+1+1处解的近似值，记作处解的近似值，记作处解的近似值，记作处解的近似值，记作y yn n+1+1[h][h]，由于公式的局部截断，由于公式的局部截断，由于公式的局部截断，由于公式的局部截断误差误差误差误差为为为为O O( (h h5 5) )，故有：，故有：，故有：，故有： 46 W Y变步长变步长R-K法（续法（续1）） 式中式中式中式中c c与与与与y y(5)(5)( (x x) )在区间在区间在区间在区间[ [x xn n , ,x xn n+1+1] ]内的值有关，假设内的值有关，假设内的值有关，假设内的值有关，假设y y(5)(5)( (x x) )在在在在区间区间区间区间[ [x xn n, ,x xn n+1+1] ]内变化不大，则可将内变化不大，则可将内变化不大，则可将内变化不大，则可将c c 视作常数。
将步长折视作常数将步长折视作常数将步长折视作常数将步长折半，即以半，即以半，即以半，即以h/h/2 2为步长，跨两步到为步长，跨两步到为步长，跨两步到为步长，跨两步到x xn n+1+1，再求得一个近似值，再求得一个近似值，再求得一个近似值，再求得一个近似值每跨一步的局部截断误差为每跨一步的局部截断误差为每跨一步的局部截断误差为每跨一步的局部截断误差为c c( (h/h/2)2)5 5，因此，因此，因此，因此 47 W Y具体做法具体做法 1. 1. 若若若若 < < 且相差不多，则取且相差不多，则取且相差不多，则取且相差不多，则取 ，并以，并以，并以，并以h h为步长为步长为步长为步长计算计算计算计算x xn n+2+2对应的对应的对应的对应的y yn n+2+23. 3. 以以以以h h作基本步长，当作基本步长，当作基本步长，当作基本步长，当 > > 时，反复将步长折半，每折半一时，反复将步长折半，每折半一时，反复将步长折半，每折半一时，反复将步长折半，每折半一 次，由次，由次，由次，由x xn n+1+1的步数增加一倍，直至的步数增加一倍，直至的步数增加一倍，直至的步数增加一倍，直至 2. 2. 如果如果如果如果 <<<< ，又不需要节点，又不需要节点，又不需要节点，又不需要节点x xn n+ +h h处的数值解值，则可将步处的数值解值，则可将步处的数值解值，则可将步处的数值解值，则可将步 长加倍，只要保证所求节点处的长加倍，只要保证所求节点处的长加倍，只要保证所求节点处的长加倍，只要保证所求节点处的y y值符合精度要求即可。
值符合精度要求即可值符合精度要求即可值符合精度要求即可48 W Y变步长四阶标准变步长四阶标准R-K法解初值问题举例法解初值问题举例 例例4解解 取基本步长取基本步长h=0.1，节点，节点xk=kh，， =10 6，按变，按变步长四阶标准步长四阶标准R-K法进行计算，法进行计算， 其结果见其结果见表表8-7，显然，尽管每一步的步长均，显然，尽管每一步的步长均只折半一次，但解的精度可达只折半一次，但解的精度可达10 6，与准确解：，与准确解：比较，说明以上分析是正确的比较，说明以上分析是正确的 表表8-7见下屏见下屏49 W Y步长折半次数步长折半次数步长折半次数步长折半次数x xk ky yk k准确解准确解准确解准确解y y( (x xk k) )1 10.10.10.1049920.1049920.104991710.104991711 10.20.20.2198680.2198680.219868070.219868071 10.30.30.3443410.3443410.344340770.344340771 10.40.40.4779530.4779530.477953480.477953481 10.50.50.6201150.6201150.62011450.62011451 10.60.60.7701350.7701350.770135280.770135281 10.70.70.9272700.9272700.927270240.927270241 10.80.81.0907541.0907541.09075361.09075361 10.90.91.2598311.2598311.25983041.25983041 11.01.01.4337811.4337811.43378081.4337808表表8-7 50 W Y§3 线性多步法线性多步法 前两节介绍的几种差分格式有一共同的前两节介绍的几种差分格式有一共同的特点，在计算特点，在计算yn+1时仅用前一步求得的时仅用前一步求得的yn，而，而对前几步的信息未予启用，因而称之为单步对前几步的信息未予启用，因而称之为单步法。
现讨论多启用一些已知信息的线性多步法现讨论多启用一些已知信息的线性多步法线性多步法公式的一般形式为：法线性多步法公式的一般形式为： 其中其中 i、、 i均为与均为与f,n无关的常数，无关的常数， r、、 r不同时为不同时为0由于求yn+1时要到前时要到前r+1步的信息步的信息yn-r，，…，，yn及对应的及对应的f 值，因称它为值，因称它为r+1步法 -1=0，，式（式（8-21）为显式，）为显式， -1 0时，式（时，式（8-21）为隐式 51 W Y3.1 阿当姆斯阿当姆斯Adams公式公式 当式当式当式当式（（（（8-218-21））））中系数中系数中系数中系数 0 0=1,=1, 1 1=…==…= r r=0=0时，称之为时，称之为时，称之为时，称之为阿当阿当阿当阿当姆斯（姆斯（姆斯（姆斯（AdamsAdams）公式）公式）公式）公式，这类公式很容易借助数值积分导，这类公式很容易借助数值积分导，这类公式很容易借助数值积分导，这类公式很容易借助数值积分导出设y y( (x x) )是初值问题（是初值问题（是初值问题（是初值问题（8-18-1）的准确解，则）的准确解，则）的准确解，则）的准确解，则: : 两边从两边从两边从两边从x xn n到到到到x xn n+1+1积分，得积分，得积分，得积分，得: : 对此式右端应用数值积分公式，则可导出不同的阿当对此式右端应用数值积分公式，则可导出不同的阿当对此式右端应用数值积分公式，则可导出不同的阿当对此式右端应用数值积分公式，则可导出不同的阿当姆斯公式。
例如，记姆斯公式例如，记姆斯公式例如，记姆斯公式例如，记f fj j= =f f ( (x xj j , ,y y( (x xj j)) ))，取数据点，取数据点，取数据点，取数据点( (x xn n-1-1, , f fn n-1-1) )，，，，( (x xn n , ,f fn n) ) 构造构造构造构造f f 的线性插值多项式的线性插值多项式的线性插值多项式的线性插值多项式: : 其中：其中：其中：其中：52 W Y阿当姆斯阿当姆斯Adams公式（续）公式（续）53 W Yr+1步阿当姆斯步阿当姆斯显式显式公式公式 （又称（又称Adams-Bashforth公式）公式） 类似地，取数据点（类似地，取数据点（类似地，取数据点（类似地，取数据点（x xn n- -r r, , f fn n- -r r）、）、）、）、( (x xn n- -r r+1+1, , f fn n- -r r+1+1) )、、、、……、、、、( (x xn n, ,f fn n) )，，，，构造构造构造构造f f ( (x x, , y y( (x x)) ))的的的的r r 次插值多项式，可导出次插值多项式，可导出次插值多项式，可导出次插值多项式，可导出r r+1+1步阿当姆斯步阿当姆斯步阿当姆斯步阿当姆斯显式公式（又称显式公式（又称显式公式（又称显式公式（又称Adams-BashforthAdams-Bashforth公式）及其局部截断误公式）及其局部截断误公式）及其局部截断误公式）及其局部截断误差差差差 : :r rA Ar rB Br r0 0B Br r1 1B Br r2 2B Br r3 3B Br r4 4 r r+1+10 01 11 1 1/21/21 12 23 3 1 1 5/125/122 212122323 16165 5 3/83/83 324245555 59593737 9 9 251/720251/7204 472072019011901 2774277426162616 1274127425125195/28895/288其系数见如下表其系数见如下表其系数见如下表其系数见如下表8-88-8，显然具有，显然具有，显然具有，显然具有r r+1+1阶精度。
阶精度 54 W Yr+1步阿当姆斯步阿当姆斯隐式隐式公式公式 （又称（又称Adams- Moulton公式）公式） 若改取数据点（若改取数据点（xn-r+1, fn-r+1）、）、(xn-r+2,fnr+2)、、…、、(xn+1, fn+1)，构造，构造f (x, y(x))的的r次插值多项式，则次插值多项式，则可可以导出以导出n+1步阿当姆斯隐式公式（又称步阿当姆斯隐式公式（又称Adams-Moulton公式公式）及其局部截断误差：）及其局部截断误差： 55 W Y 其系数见下其系数见下其系数见下其系数见下表表表表8-98-9，显然也具有，显然也具有，显然也具有，显然也具有r r+1+1阶精度 r rA Ar rB Br r-1-1* *B Br r0 0* *B Br r1 1* *B Br r2 2* *B Br r3 3* * r+1r+1* *0 01 11 1 1/21/21 12 21 11 1 1/121/122 212125 58 8 1 1 1/241/243 324249 91919 5 51 1 19/72019/7204 4720720251251646646 264264106106 1919 3/1603/160可见，一阶阿当姆斯显式公式（对应可见，一阶阿当姆斯显式公式（对应r=0）即为欧）即为欧拉公式；一阶与二阶阿当姆斯隐式公式对应拉公式；一阶与二阶阿当姆斯隐式公式对应r=0与与r=1，分别为后退欧拉公式与梯形公式。
分别为后退欧拉公式与梯形公式 r+1步阿当姆斯步阿当姆斯隐式隐式公式公式 （又称（又称Adams- Moulton公式）（续）公式）（续） 56 W Y四阶阿当姆斯显式公式四阶阿当姆斯显式公式 高阶阿当姆斯公式中最常用的是四阶公式，对应高阶阿当姆斯公式中最常用的是四阶公式，对应高阶阿当姆斯公式中最常用的是四阶公式，对应高阶阿当姆斯公式中最常用的是四阶公式，对应r r=3=3若选选选选x xn-3 n-3 , x, xn-2n-2，，，，x xn-1 n-1 , x, xn n作为插值节点，则四阶阿当姆斯显式作为插值节点，则四阶阿当姆斯显式作为插值节点，则四阶阿当姆斯显式作为插值节点，则四阶阿当姆斯显式公式（又称为外推公式）及其局部截断误差为公式（又称为外推公式）及其局部截断误差为公式（又称为外推公式）及其局部截断误差为公式（又称为外推公式）及其局部截断误差为 ：：：： 此公式的优点是精度高，其局部截断误差与四阶此公式的优点是精度高，其局部截断误差与四阶此公式的优点是精度高，其局部截断误差与四阶此公式的优点是精度高，其局部截断误差与四阶R-KR-K法法法法同阶，但每前进一步只要计算一次函数值同阶，但每前进一步只要计算一次函数值同阶，但每前进一步只要计算一次函数值同阶，但每前进一步只要计算一次函数值f f ( (x xn n, ,y yn n) )，其余，其余，其余，其余的值的值的值的值f f 在前面的计算中已求出，计算量远小于在前面的计算中已求出，计算量远小于在前面的计算中已求出，计算量远小于在前面的计算中已求出，计算量远小于四阶四阶四阶四阶R-KR-K法法法法。
不足之处是不能从不足之处是不能从不足之处是不能从不足之处是不能从n=1n=1开始使用，必须用其它方法（例如开始使用，必须用其它方法（例如开始使用，必须用其它方法（例如开始使用，必须用其它方法（例如四阶四阶四阶四阶R-KR-K法）求得出发值法）求得出发值法）求得出发值法）求得出发值y y1 1、、、、y y2 2、、、、y y3 3此外，若要中途改变此外，若要中途改变此外，若要中途改变此外，若要中途改变 步长也较麻烦步长也较麻烦步长也较麻烦步长也较麻烦 57 W Y四阶阿当姆隐式公式内插公式及其局部截断误差为四阶阿当姆隐式公式内插公式及其局部截断误差为四阶阿当姆隐式公式内插公式及其局部截断误差为四阶阿当姆隐式公式内插公式及其局部截断误差为（取（取（取（取x xn-2n-2,x ,xn-1n-1,x ,xn n,x ,xn+1n+1为插值节点）：为插值节点）：为插值节点）：为插值节点）： 它也具有公式（它也具有公式（它也具有公式（它也具有公式（8-268-26）的优缺点，只是对）的优缺点，只是对）的优缺点，只是对）的优缺点，只是对y yn+1n+1必须提供必须提供必须提供必须提供一个预测值，然后通过迭代才能求得一个预测值，然后通过迭代才能求得一个预测值，然后通过迭代才能求得一个预测值，然后通过迭代才能求得y yn+1n+1。
将公式（将公式（将公式（将公式（8-8-2626）与公式（）与公式（）与公式（）与公式（8-278-27）联合使用，可组成如下预测一校正系）联合使用，可组成如下预测一校正系）联合使用，可组成如下预测一校正系）联合使用，可组成如下预测一校正系统：统：统：统： 四阶阿当姆斯显式公式（续）四阶阿当姆斯显式公式（续）58 W Y阿当姆斯公式举例阿当姆斯公式举例例例例例5 5 分别用四阶阿当姆斯显式公式（分别用四阶阿当姆斯显式公式（分别用四阶阿当姆斯显式公式（分别用四阶阿当姆斯显式公式（8-268-26）与阿当姆斯预）与阿当姆斯预）与阿当姆斯预）与阿当姆斯预测一校正公式测一校正公式测一校正公式测一校正公式 （（（（8-288-28）解例）解例）解例）解例3 3中的初值问题中的初值问题中的初值问题中的初值问题 解解解解 取步长取步长取步长取步长h h=0.1=0.1，由于公式（，由于公式（，由于公式（，由于公式（8-98-9）、（）、（）、（）、（8-288-28）均为四步法，）均为四步法，）均为四步法，）均为四步法，因此必须用其它方法求出发值。
现采用四阶标准龙格因此必须用其它方法求出发值现采用四阶标准龙格因此必须用其它方法求出发值现采用四阶标准龙格因此必须用其它方法求出发值现采用四阶标准龙格- -库库库库塔法求塔法求塔法求塔法求y y1 1 ,y ,y2 2 , y, y3 ‘3 ‘然后再分别用指定的方法进行计算，计然后再分别用指定的方法进行计算，计然后再分别用指定的方法进行计算，计然后再分别用指定的方法进行计算，计算结果见表算结果见表算结果见表算结果见表8-108-10 表表表表8-108-10x xn ny yn n四阶标准四阶标准四阶标准四阶标准R-KR-K法法法法显式公式显式公式显式公式显式公式（（（（8-258-25））））预测预测预测预测- -校正公式（校正公式（校正公式（校正公式（8-288-28））））准确解准确解准确解准确解0.10.11.0954461.095446 1.0954451151.0954451150.20.21.1832171.183217 1.1832159561.1832159560.30.31.2649121.264912 1.2649110461.2649110460.40.4 1.3415521.3415521.3416411.3416411.3416407861.3416407860.50.5 1.4140461.4140461.4142131.4142131.4142135621.4142135620.60.6 1.4830191.4830191.4832391.4832391.4832396871.4832396870.70.7 1.5489191.5489191.5491921.5491921.5481833381.5481833380.80.8 1.6121161.6121161.6124501.6124501.6124515491.6124515490.90.9 1.6729171.6729171.6733781.6733781.633200531.633200531.01.0 1.7315691.7315691.7320481.7320481.7320508071.73205080759 W Y阿当姆斯预测阿当姆斯预测—校正公式校正公式 为了进一步提高为了进一步提高为了进一步提高为了进一步提高上述上述上述上述预测一校正法的计算精度，可以预测一校正法的计算精度，可以预测一校正法的计算精度，可以预测一校正法的计算精度，可以利用预测公式与校正公式同阶的特点，导出局部截断误差利用预测公式与校正公式同阶的特点，导出局部截断误差利用预测公式与校正公式同阶的特点，导出局部截断误差利用预测公式与校正公式同阶的特点，导出局部截断误差估计公式，对预测值与校正值分别用各自的局部截断误差估计公式，对预测值与校正值分别用各自的局部截断误差估计公式，对预测值与校正值分别用各自的局部截断误差估计公式，对预测值与校正值分别用各自的局部截断误差估计值进行修正，以此来提高解的精度。
由式估计值进行修正，以此来提高解的精度由式估计值进行修正，以此来提高解的精度由式估计值进行修正，以此来提高解的精度由式（（（（8-268-26））））、、、、（（（（8-278-27））））可知，公式可知，公式可知，公式可知，公式（（（（8-288-28））））中预测公式与校正公式的局中预测公式与校正公式的局中预测公式与校正公式的局中预测公式与校正公式的局部截断误差分别为部截断误差分别为部截断误差分别为部截断误差分别为 ：：：：若记预测值为若记预测值为若记预测值为若记预测值为p pn n+1+1，校正值为，校正值为，校正值为，校正值为c cn n+1+1，则：，则：，则：，则：紧接下屏紧接下屏紧接下屏紧接下屏 60 W Y于是得于是得pn+1与与cn+1的事后误差估计：的事后误差估计： 阿当姆斯预测阿当姆斯预测—校正公式（续）校正公式（续）分别代替分别代替pn+1与与cn+1能能提高解的精度为简提高解的精度为简单计，在未求出单计，在未求出cn+1以前，用以前，用pn cn代替代替pn+1 cn+1 61 W Y阿当姆斯带误差修正的预测阿当姆斯带误差修正的预测—校正公式校正公式据此，得到下述带误差修正的预测一校正公式据此，得到下述带误差修正的预测一校正公式据此，得到下述带误差修正的预测一校正公式据此，得到下述带误差修正的预测一校正公式 ：：：：式中初始值式中初始值式中初始值式中初始值p p3 3与与与与c c3 3可取为可取为可取为可取为0 0。
62 W Y3.2 哈明哈明(Hamming)公式公式 一般的线性多步公式（一般的线性多步公式（一般的线性多步公式（一般的线性多步公式（8-218-21）还可用待系数法并借助泰）还可用待系数法并借助泰）还可用待系数法并借助泰）还可用待系数法并借助泰勒展开导出例如考虑确定如下形式的四阶公式勒展开导出例如考虑确定如下形式的四阶公式勒展开导出例如考虑确定如下形式的四阶公式勒展开导出例如考虑确定如下形式的四阶公式 ：：：：其中其中其中其中y y n n表示表示表示表示f fn n只需求系数只需求系数只需求系数只需求系数 j j、、、、 j j，使其局部截断误差，使其局部截断误差，使其局部截断误差，使其局部截断误差T T= =O O( (h h5 5) ) 即可局部截断误差为局部截断误差为局部截断误差为局部截断误差为 ：：：：将其在将其在将其在将其在x xn n点进点进点进点进行泰勒展开，得：行泰勒展开，得：行泰勒展开，得：行泰勒展开，得： （紧接下屏）（紧接下屏）（紧接下屏）（紧接下屏）63 W Y哈明哈明(Hamming)公式（续公式（续1））令令令令h hk k( (k k=0,1,…4) =0,1,…4) 的系数为的系数为的系数为的系数为0, 0,得关于得关于得关于得关于 j j、、、、 j j的方程组的方程组的方程组的方程组 : :方程组有六个未知数，而只有五个方程，因此有无穷多方程组有六个未知数，而只有五个方程，因此有无穷多方程组有六个未知数，而只有五个方程，因此有无穷多方程组有六个未知数，而只有五个方程，因此有无穷多解，取解，取解，取解，取 1 1为自由参量，则有为自由参量，则有为自由参量，则有为自由参量，则有: : 64 W Y哈明哈明(Hamming)公式（续公式（续2））当当当当 1 1=0=0时，得公式：时，得公式：时，得公式：时，得公式： 这是一个四阶隐式公式。
这是一个四阶隐式公式这是一个四阶隐式公式这是一个四阶隐式公式 完全类似，可导出如下形式的四阶显式公式：完全类似，可导出如下形式的四阶显式公式：完全类似，可导出如下形式的四阶显式公式：完全类似，可导出如下形式的四阶显式公式：其中最常其中最常其中最常其中最常用的是：用的是：用的是：用的是：65 W YHamming公式公式（预测一校正公式）（预测一校正公式）（预测一校正公式）（预测一校正公式） 用显式公式（用显式公式（用显式公式（用显式公式（8-348-34）与隐式公式（）与隐式公式（）与隐式公式（）与隐式公式（8-338-33）联合可组成）联合可组成）联合可组成）联合可组成预测一校正公式预测一校正公式预测一校正公式预测一校正公式 此式称为此式称为此式称为此式称为哈明（哈明（哈明（哈明（HammingHamming）公式）公式）公式）公式，大量数值实验表，大量数值实验表，大量数值实验表，大量数值实验表明，此公式的数值稳定性在同类公式中几乎是最好的明，此公式的数值稳定性在同类公式中几乎是最好的明，此公式的数值稳定性在同类公式中几乎是最好的。
明，此公式的数值稳定性在同类公式中几乎是最好的 为了进一步改善哈明公式计算结果的精度而又不致为了进一步改善哈明公式计算结果的精度而又不致为了进一步改善哈明公式计算结果的精度而又不致为了进一步改善哈明公式计算结果的精度而又不致增加过多的计算量，可仿照带误差修正的阿当姆斯预测增加过多的计算量，可仿照带误差修正的阿当姆斯预测增加过多的计算量，可仿照带误差修正的阿当姆斯预测增加过多的计算量，可仿照带误差修正的阿当姆斯预测一校正公式，建立一个带误差修正的哈明公式一校正公式，建立一个带误差修正的哈明公式一校正公式，建立一个带误差修正的哈明公式一校正公式，建立一个带误差修正的哈明公式 ：：：：（紧接下屏）（紧接下屏）（紧接下屏）（紧接下屏）66 W Y 前面介绍的四套预测一校正公式（前面介绍的四套预测一校正公式（前面介绍的四套预测一校正公式（前面介绍的四套预测一校正公式（8-288-28）、（）、（）、（）、（8-298-29）、）、）、）、（（（（8-328-32）、（）、（）、（）、（8-338-33）都是实际计算中常用的方法，它们的）都是实际计算中常用的方法，它们的）都是实际计算中常用的方法，它们的）都是实际计算中常用的方法，它们的 共同特点是计算量省而精度高，但都必须用别的单步法共同特点是计算量省而精度高，但都必须用别的单步法共同特点是计算量省而精度高，但都必须用别的单步法共同特点是计算量省而精度高，但都必须用别的单步法 为其提供出发值，由于解题过程是递推进行的，所以出为其提供出发值，由于解题过程是递推进行的，所以出为其提供出发值，由于解题过程是递推进行的，所以出为其提供出发值，由于解题过程是递推进行的，所以出 发值的精度对以后的计算影响很大，故一般采用发值的精度对以后的计算影响很大，故一般采用发值的精度对以后的计算影响很大，故一般采用发值的精度对以后的计算影响很大，故一般采用四阶龙四阶龙四阶龙四阶龙 格一库塔法求其出发值。
格一库塔法求其出发值格一库塔法求其出发值格一库塔法求其出发值 Hamming公式公式（预测一校正公式）（续）（预测一校正公式）（续）（预测一校正公式）（续）（预测一校正公式）（续） 67 W Y§4 一阶方程组与高阶方程初值问题一阶方程组与高阶方程初值问题 4.1 4.1 一阶方程组一阶方程组一阶方程组一阶方程组 下面以两个方程的情形为例，介绍一阶方程组的数值解下面以两个方程的情形为例，介绍一阶方程组的数值解下面以两个方程的情形为例，介绍一阶方程组的数值解下面以两个方程的情形为例，介绍一阶方程组的数值解法设初值问题：法设初值问题：法设初值问题：法设初值问题： 若采用向量记号，记：若采用向量记号，记：若采用向量记号，记：若采用向量记号，记：则初值问题（则初值问题（则初值问题（则初值问题（8-348-34）可表示为：）可表示为：）可表示为：）可表示为： 它与初值问题（它与初值问题（它与初值问题（它与初值问题（8-18-1）形式上完全一致，前面介绍的解初）形式上完全一致，前面介绍的解初）形式上完全一致，前面介绍的解初）形式上完全一致，前面介绍的解初值问题（值问题（值问题（值问题（8-18-1）的各种数值方法均可用于解一阶方程组初）的各种数值方法均可用于解一阶方程组初）的各种数值方法均可用于解一阶方程组初）的各种数值方法均可用于解一阶方程组初值问题（值问题（值问题（值问题（8-358-35），只要将公式中的），只要将公式中的），只要将公式中的），只要将公式中的y y与与与与f f换成向量换成向量换成向量换成向量y y与与与与f f即可。
即可68 W Y一阶方程组（续一阶方程组（续1））例如：解初值问题例如：解初值问题例如：解初值问题例如：解初值问题（（（（8-358-35）的欧拉公式为：）的欧拉公式为：）的欧拉公式为：）的欧拉公式为： 即：即：即：即：写成分量形式即为：写成分量形式即为：写成分量形式即为：写成分量形式即为：又如解初值问题又如解初值问题又如解初值问题又如解初值问题（（（（8-358-35）的四阶）的四阶）的四阶）的四阶标准龙格标准龙格标准龙格标准龙格- -库塔库塔库塔库塔公式为：公式为：公式为：公式为：69 W Y一阶方程组（续一阶方程组（续2）） 写成分量形式为：写成分量形式为：写成分量形式为：写成分量形式为： 70 W Y4.2 高阶方程初值问题高阶方程初值问题 高阶方程初值问题常转化为一阶方程组求解例如：高阶方程初值问题常转化为一阶方程组求解例如：高阶方程初值问题常转化为一阶方程组求解例如：高阶方程初值问题常转化为一阶方程组求解例如： 引进新变量引进新变量引进新变量引进新变量z=z=y y ，则此问题等价于，则此问题等价于，则此问题等价于，则此问题等价于 对其使用四阶标龙格对其使用四阶标龙格对其使用四阶标龙格对其使用四阶标龙格- -库塔法，计算公式为库塔法，计算公式为库塔法，计算公式为库塔法，计算公式为 71 W Y高阶方程初值问题举例高阶方程初值问题举例例例例例5 5 用四阶龙格用四阶龙格用四阶龙格用四阶龙格- -库塔方法在库塔方法在库塔方法在库塔方法在[0[0，，，，1]1]上取步长上取步长上取步长上取步长h h=0.2=0.2，求解，求解，求解，求解二阶方程初值问题：二阶方程初值问题：二阶方程初值问题：二阶方程初值问题： 解解解解 先将二阶方程化为一先将二阶方程化为一先将二阶方程化为一先将二阶方程化为一 阶方程组。
令阶方程组令阶方程组令阶方程组令z=yz=y ，，，， 则得方程组：则得方程组：则得方程组：则得方程组：使用四阶龙格使用四阶龙格使用四阶龙格使用四阶龙格- -库塔公库塔公库塔公库塔公式，其相应的形式为：式，其相应的形式为：式，其相应的形式为：式，其相应的形式为：其中其中其中其中 72 W Y例例5（续）（续）表表表表8-118-11 y y -3-3y y +2+2y y=0=0解的近拟值解的近拟值解的近拟值解的近拟值 x xi iy yi iz zi ik k1 1l l1 1k k2 2l l2 2k k3 3l l3 3k k4 4l l4 40 01.00001.00001.00001.00001.00001.00001.00001.00001.10001.10001.10001.10001.11001.11001.11001.11001.22201.22201.22201.22200.20.21.22141.22141.22141.22141.22141.22141.22141.22141.34351.34351.34351.34351.35581.35581.35581.35581.49261.49261.49261.49260.40.41.19181.19181.19181.19181.19181.19181.19181.19181.64101.64101.64101.64101.65591.65591.65591.65591.82301.82301.82301.82300.60.61.82211.82211.82211.82211.82211.82211.82211.82212.00432.00432.00432.00432.02252.02252.02252.02252.22662.22662.22662.22660.80.82.22552.22552.22552.22552.22552.22552.22552.22552.44812.44812.44812.44812.47032.47032.47032.47032.71962.71962.71962.71961.01.02.71822.7182 计算结果列于表计算结果列于表计算结果列于表计算结果列于表8-118-11中：中：中：中： 与准确解与准确解与准确解与准确解y y= =e e x x比较，可知所得数值解除比较，可知所得数值解除比较，可知所得数值解除比较，可知所得数值解除x x=1=1外，外，外，外，均有五位有效数。
均有五位有效数均有五位有效数均有五位有效数 73 W Y§5 收敛性与稳定性收敛性与稳定性 通过前面的讨论可以看到，微分方通过前面的讨论可以看到，微分方程数值解法的程数值解法的基本思想基本思想是：是： 通过某种离散化手段，将微分方程通过某种离散化手段，将微分方程转化为差分格式求解从理论上说，这转化为差分格式求解从理论上说，这里需要解决两个问题，一是当步长里需要解决两个问题，一是当步长h0时，差分解，这就是差分格式的收敛性时，差分解，这就是差分格式的收敛性问题；二是利用差分格式求解时，初始问题；二是利用差分格式求解时，初始误差及计算过程中的舍入误差能否得到误差及计算过程中的舍入误差能否得到控制，这就是差分格式的稳定性问题控制，这就是差分格式的稳定性问题本节对这两个问题作一简单介绍本节对这两个问题作一简单介绍 74 W Y5.1 收敛性收敛性 用某种差分格式解微分方程时，若对求解区间用某种差分格式解微分方程时，若对求解区间用某种差分格式解微分方程时，若对求解区间用某种差分格式解微分方程时，若对求解区间[a,b][a,b]中的中的中的中的任一点任一点任一点任一点x x，，，，当当当当h h0 0时差分解时差分解时差分解时差分解y yn n一致收敛一致收敛一致收敛一致收敛于微分方程之解于微分方程之解于微分方程之解于微分方程之解y(x)y(x)，，，，对对对对x x一致成立，则称该差分格式是收敛的。
一致成立，则称该差分格式是收敛的一致成立，则称该差分格式是收敛的一致成立，则称该差分格式是收敛的 显然，只有收敛的差分格式才有实用价值可以证明，显然，只有收敛的差分格式才有实用价值可以证明，显然，只有收敛的差分格式才有实用价值可以证明，显然，只有收敛的差分格式才有实用价值可以证明，前几节介绍的单步法与多步法对满足解的存在唯一性定理前几节介绍的单步法与多步法对满足解的存在唯一性定理前几节介绍的单步法与多步法对满足解的存在唯一性定理前几节介绍的单步法与多步法对满足解的存在唯一性定理条件的微分方程都是收敛的条件的微分方程都是收敛的条件的微分方程都是收敛的条件的微分方程都是收敛的 对于收敛的差分格式，从理论上说只要对于收敛的差分格式，从理论上说只要对于收敛的差分格式，从理论上说只要对于收敛的差分格式，从理论上说只要h h充分小，差充分小，差充分小，差充分小，差分解分解分解分解y yn n的精度就可以任意高的精度就可以任意高的精度就可以任意高的精度就可以任意高75 W Y5.2 稳定性稳定性 实际计算中舍入误差几乎不可避免，这类误差在递实际计算中舍入误差几乎不可避免，这类误差在递实际计算中舍入误差几乎不可避免，这类误差在递实际计算中舍入误差几乎不可避免，这类误差在递推计算过程中会不会恶性增长，以致淹没了差分方程的推计算过程中会不会恶性增长，以致淹没了差分方程的推计算过程中会不会恶性增长，以致淹没了差分方程的推计算过程中会不会恶性增长，以致淹没了差分方程的真解，这是微分方程数值解法中的另一个重问题。
由于真解，这是微分方程数值解法中的另一个重问题由于真解，这是微分方程数值解法中的另一个重问题由于真解，这是微分方程数值解法中的另一个重问题由于研究的对象不同，关于稳定性的定义很多，这里只介绍研究的对象不同，关于稳定性的定义很多，这里只介绍研究的对象不同，关于稳定性的定义很多，这里只介绍研究的对象不同，关于稳定性的定义很多，这里只介绍其中一种其中一种其中一种其中一种────绝对稳定性绝对稳定性绝对稳定性绝对稳定性 由于实际问题中微分方程（由于实际问题中微分方程（由于实际问题中微分方程（由于实际问题中微分方程（8-18-1）的右端）的右端）的右端）的右端f(x,y)f(x,y)复杂多样，复杂多样，复杂多样，复杂多样，这给研究稳定性带来了困难，因此通常取以下模型方程这给研究稳定性带来了困难，因此通常取以下模型方程这给研究稳定性带来了困难，因此通常取以下模型方程这给研究稳定性带来了困难，因此通常取以下模型方程为基础进行研究为基础进行研究为基础进行研究为基础进行研究 其中为常数，其中为常数，其中为常数，其中为常数， λ λ <0 <0表示微分方程本身是稳定的表示微分方程本身是稳定的。
表示微分方程本身是稳定的表示微分方程本身是稳定的 若取步长若取步长若取步长若取步长h h，用某种差分格式解模型方程时，任何一，用某种差分格式解模型方程时，任何一，用某种差分格式解模型方程时，任何一，用某种差分格式解模型方程时，任何一步产生的舍入误差在以后的计算中引起的误差逐渐减弱，步产生的舍入误差在以后的计算中引起的误差逐渐减弱，步产生的舍入误差在以后的计算中引起的误差逐渐减弱，步产生的舍入误差在以后的计算中引起的误差逐渐减弱，则称该差分格式关于则称该差分格式关于则称该差分格式关于则称该差分格式关于z=z=h h是绝对稳定的，使差分格式稳是绝对稳定的，使差分格式稳是绝对稳定的，使差分格式稳是绝对稳定的，使差分格式稳定的所有定的所有定的所有定的所有z z值的集合称为绝对值的集合称为绝对值的集合称为绝对值的集合称为绝对稳定区间稳定区间稳定区间稳定区间 76 W Y稳定性（续稳定性（续1））例如，将欧拉公式用于模型方程例如，将欧拉公式用于模型方程例如，将欧拉公式用于模型方程例如，将欧拉公式用于模型方程y y = = y y，得，得，得，得 不妨设在不妨设在不妨设在不妨设在y yn n上有误差上有误差上有误差上有误差 n n，则由此引起，则由此引起，则由此引起，则由此引起y yn+1n+1的误差的误差的误差的误差 n+1n+1满足满足满足满足: : 显然，要保证误差逐渐减弱，必须使显然，要保证误差逐渐减弱，必须使显然，要保证误差逐渐减弱，必须使显然，要保证误差逐渐减弱，必须使z=z=h h满足满足满足满足: : 77 W Y表表8-12列出了几个常用公式的绝对稳定区间。
列出了几个常用公式的绝对稳定区间从表中可以看到，隐式公式比同阶显式公式的稳从表中可以看到，隐式公式比同阶显式公式的稳定区间大，这就是为什么最好立足于隐式公式求定区间大，这就是为什么最好立足于隐式公式求解的原因解的原因 龙格龙格龙格龙格- -库塔公式库塔公式库塔公式库塔公式（显式）（显式）（显式）（显式）阿当姆斯公式阿当姆斯公式阿当姆斯公式阿当姆斯公式显式显式显式显式隐式隐式隐式隐式2 2阶阶阶阶3 3阶阶阶阶4 4阶阶阶阶（－（－（－（－2, 02, 0））））（－（－（－（－2.51, 02.51, 0））））（－（－（－（－2.78, 02.78, 0））））（－（－（－（－1, 01, 0））））（－（－（－（－6/11, 06/11, 0））））（－（－（－（－3/10, 03/10, 0））））（－（－（－（－ , 0, 0））））（－（－（－（－6, 06, 0））））（－（－（－（－3, 03, 0））））表表表表8-128-12稳定性（续稳定性（续1））78 W Y稳定性（续稳定性（续3））研究一般方程的稳定性时，研究一般方程的稳定性时，研究一般方程的稳定性时，研究一般方程的稳定性时， 相当于模型方程中的相当于模型方程中的相当于模型方程中的相当于模型方程中的λ λ ，，，，例如对欧拉公式例如对欧拉公式例如对欧拉公式例如对欧拉公式: : 设设设设y yn n有误差有误差有误差有误差 n n，则由此引起，则由此引起，则由此引起，则由此引起y yn n+1+1的误差的误差的误差的误差 n n+1+1满足满足满足满足: :所以有：所以有：所以有：所以有：（紧接下屏）（紧接下屏）（紧接下屏）（紧接下屏）79 W Y综上所述，综上所述，综上所述，综上所述，当用数值方法解微分方程时，步长当用数值方法解微分方程时，步长当用数值方法解微分方程时，步长当用数值方法解微分方程时，步长h h的选择由的选择由的选择由的选择由两个因素确定。
一是精度要求，由收敛性可知，两个因素确定一是精度要求，由收敛性可知，两个因素确定一是精度要求，由收敛性可知，两个因素确定一是精度要求，由收敛性可知，h h越小精越小精越小精越小精度越高 二是稳定性要求，二是稳定性要求，二是稳定性要求，二是稳定性要求， 应落在所选公式的应落在所选公式的应落在所选公式的应落在所选公式的绝对稳定域中虽然绝对稳定域中虽然绝对稳定域中虽然绝对稳定域中虽然h h十分小时两个条件都满足，但十分小时两个条件都满足，但十分小时两个条件都满足，但十分小时两个条件都满足，但h h太小太小太小太小计算步数势必大增，每步的微小误差累积起来也很可观计算步数势必大增，每步的微小误差累积起来也很可观计算步数势必大增，每步的微小误差累积起来也很可观计算步数势必大增，每步的微小误差累积起来也很可观因此，实际计算中在满足精度要求与稳定性的条件下，步因此，实际计算中在满足精度要求与稳定性的条件下，步因此，实际计算中在满足精度要求与稳定性的条件下，步因此，实际计算中在满足精度要求与稳定性的条件下，步长长长长h h应尽可能取得大些。
应尽可能取得大些应尽可能取得大些应尽可能取得大些 式中式中式中式中 n n在在在在y yn n+ + n n与与与与y yn n之间，与式（之间，与式（之间，与式（之间，与式（8-378-37）相比，可知）相比，可知）相比，可知）相比，可知 起着起着起着起着 的作用，因此为了保证数值方法的绝对稳定，的作用，因此为了保证数值方法的绝对稳定，的作用，因此为了保证数值方法的绝对稳定，的作用，因此为了保证数值方法的绝对稳定， 应落在它的绝对稳定域中应落在它的绝对稳定域中应落在它的绝对稳定域中应落在它的绝对稳定域中 稳定性（续稳定性（续4））80 W Y第八章第八章结结 束束81 。