数值分析09-常微方程数值解法.ppt

W Y 第八章 常微分方程常微分方程 数值解法数值解法 8-1阜师院数科院第八章 常微分方程数值解法 W Y 第八章目录 1 1 欧拉（欧拉（EulerEuler）方法方法 1.1 1.1 EulerEuler法及其简单改进法及其简单改进 1.2 1.2 改进的改进的EulerEuler法法 2 2 龙格库塔（龙格库塔（Runge-kuttaRunge-kutta）方法方法 2.1 2.1 龙格龙格- -库塔方法的基本思想库塔方法的基本思想 2.2 2.2 二阶龙格二阶龙格- -库塔公式库塔公式 2.3 2.3 高阶高阶R-KR-K公式公式 2.4 2.4 变步长变步长R-KR-K法法 3 3 线性多步法线性多步法 4 4 一阶方程组与高阶方程初值问题一阶方程组与高阶方程初值问题 5 5 收敛性与稳定性收敛性与稳定性 2阜师院数科院第八章 常微分方程数值解法 W Y 第八章 序 许多科学技术问题，例如天文学中的星体运动，空间许多科学技术问题，例如天文学中的星体运动，空间 技术中的物体飞行，自动控制中的系统分析，力学中的振技术中的物体飞行，自动控制中的系统分析，力学中的振 动，工程问题中的电路分析等，都可归结为常微分方程的动，工程问题中的电路分析等，都可归结为常微分方程的 初值问题。初值问题。 所谓初值问题，是函数及其必要的导数在积分的起始所谓初值问题，是函数及其必要的导数在积分的起始 点为已知的一类问题，一般形式为：点为已知的一类问题，一般形式为： 我们先介绍我们先介绍 简单的一阶问题：简单的一阶问题： 3阜师院数科院第八章 常微分方程数值解法 W Y 第八章 序 由由常微分方程常微分方程的理论可知：上述问题的解唯一存在。的理论可知：上述问题的解唯一存在。 常微分方程常微分方程求解求什么？应求一满足初值问题（求解求什么？应求一满足初值问题（8181） 的解函数的解函数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(x i i ) ) （离散点）离散点）: : 上上 y(xy(x i i ) ) 的的近似值近似值y y i i （ i=i=1,2,n1,2,n），），并且还可由这些（并且还可由这些（ x x i i , ,y y i i ) ) （i=i=1,2,n1,2,n）构造插值函数作为近似函数。上述离散点相构造插值函数作为近似函数。上述离散点相 邻两点间的距离邻两点间的距离h h i i =x=xi i-1 -1 - -x x i i 称为步长，若称为步长，若h h i i 都相等为一定数都相等为一定数 h, h, 则称为定步长，否则为变步长。则称为定步长，否则为变步长。 由于在实际问题和科学研究中遇到的微分方程往往很由于在实际问题和科学研究中遇到的微分方程往往很 复杂，复杂，绝大多数很难，甚至不可能求出解析绝大多数很难，甚至不可能求出解析函数函数y y( (x x) )，因因 此只能考虑求其数值解。此只能考虑求其数值解。 本章重点讨论如下本章重点讨论如下 一阶微分方程：一阶微分方程： 在此基础上介绍一阶微分方程组与在此基础上介绍一阶微分方程组与 高阶微分方程的数值解法。高阶微分方程的数值解法。 5阜师院数科院第八章 常微分方程数值解法 W Y 1 欧拉（Euler）法 以以EulerEuler法及其改进方法为例法及其改进方法为例, ,说明说明 常微分方程常微分方程初值问题数值解法的一般概初值问题数值解法的一般概 念，念，EulerEuler法很简单，准确度也不高，法很简单，准确度也不高， 介绍此方法的目的，是由于对它的分析介绍此方法的目的，是由于对它的分析 讨论能够比较清楚地显示出方法的一些讨论能够比较清楚地显示出方法的一些 特点，而这些特点及基本方法反映了其特点，而这些特点及基本方法反映了其 它方法的特点。它方法的特点。 EulerEuler法法用于求用于求 解一阶微分方解一阶微分方 程初值问题：程初值问题： 6阜师院数科院第八章 常微分方程数值解法 W Y 1.1 Euler法及其简单改进 EulerEuler公公 式为：式为： 由由x x 0 0 出发出发x x 1 1, ,x x2 2 ,x x N N ，而利用此式可算出对应的而利用此式可算出对应的 y y1 1, ,y y2 2 ,y y N N ，式（式（8-28-2）称为）称为差分方程差分方程（序列（序列 y y n 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 P 0 0 ( (x x 0 0, ,y y0 0 ) ) 出发，以出发，以P P 0 0 为切点，作切线，切线斜率为曲线为切点，作切线，切线斜率为曲线y y( (x x) )的切线斜率的切线斜率 y y = =f f ( (x x 0 0, ,y y0 0 ) )，因此可因此可 得切线：（点斜式）得切线：（点斜式） P1 P2 y(x) P0 x2x1x0 7阜师院数科院第八章 常微分方程数值解法 W Y Euler公式的推导(续1） 几何意义几何意义：用折线近似解曲线：用折线近似解曲线y y = = y y( (x x) )，折线不会偏离太远折线不会偏离太远 ，因为每项以，因为每项以f f ( (x x, , y y) )（斜率）修正。斜率）修正。 切线与切线与x x = = x x 1 1 交于交于P P 1 1 ( (x x 1 1, ,y y1 1 ) )，在在 x x 0 0, ,x x1 1 上以切线上以切线 近似曲线，近似曲线， 8阜师院数科院第八章 常微分方程数值解法 W Y Euler公式的推导（续2） 二二、利用、利用TaylorTaylor级数：将级数：将y y( (x x) )在在x x n n 处展开：处展开： 9阜师院数科院第八章 常微分方程数值解法 W Y Euler公式的推导（续3） 公式（公式（8-38-3）称为后退）称为后退EulerEuler公式公式 所谓局部载断误差是指以所谓局部载断误差是指以y y n n 代替代替y y( (x x n n ) )而得到而得到y y( (x xn n+1 +1) )的近 的近 似值似值y yn n+1 +1的误差（只估计这一步的误差）。 的误差（只估计这一步的误差）。 三、利用数值微分公式：利用两点公式三、利用数值微分公式：利用两点公式 并且并且EulerEuler公式的公式的 局部截断误差为局部截断误差为: : 后退后退EulerEuler公式的公式的 局部截断误差为局部截断误差为: : 10阜师院数科院第八章 常微分方程数值解法 W Y Euler公式的推导（续4） 11阜师院数科院第八章 常微分方程数值解法 W Y Euler公式的推导（续5） 四四、利用数值积分公式：在利用数值积分公式：在 x x n n, ,x x n n+1+1 上对 上对y y ( (x x)=)=f f ( (x x, ,y y( (x x) ) 积分积分 对右端积分项采用不同的数值积分公式，便可得到各种对右端积分项采用不同的数值积分公式，便可得到各种 不同的求解不同的求解dEdE初值问题的计算公式。初值问题的计算公式。 如，以矩形面积代替曲边梯形面积如，以矩形面积代替曲边梯形面积 1 1）以左矩形面积代替曲边梯形面积如图以左矩形面积代替曲边梯形面积如图8-28-2，亦即以，亦即以 yf(x, y) xn x xn+1 图8-2 12阜师院数科院第八章 常微分方程数值解法 W Y y f(x, y) xn x xn+1图8-3 y f(x, y) xn x xn+1 图8-4 3 3）以梯形公式（面积）代替曲边形如图）以梯形公式（面积）代替曲边形如图8-48-4则有则有 式（式（8-58-5）称为求）称为求dEdE初值问题的梯形公式，梯形公式看作初值问题的梯形公式，梯形公式看作 是以（是以（x x n n, ,y yn n ）( (x xn n+1 +1 , ,y y n n+1+1) )构造的插值多项式代替被积函数 构造的插值多项式代替被积函数 得到的，而得到的，而EulerEuler公式公式则是以左端点函数值近似被积函数则是以左端点函数值近似被积函数 而得到，还可以用多个点做插值多项式近似被积函数构造而得到，还可以用多个点做插值多项式近似被积函数构造 另一些精度更高的解微分方程的数值公式，梯形公式比另一些精度更高的解微分方程的数值公式，梯形公式比 EulerEuler公式公式更准确一些，误差更小。更准确一些，误差更小。 EulerEuler公式的推导（续公式的推导（续6 6） 2）以右矩形面积代替曲边梯形：如图8-3亦即以 13阜师院数科院第八章 常微分方程数值解法 W Y Euler公式注释 注注1 1 ：EulerEuler公式为显式，右矩形，梯形公式为隐式；公式为显式，右矩形，梯形公式为隐式； 注注2 2：EulerEuler公式，梯形，右矩形公式为单步法，计算公式，梯形，右矩形公式为单步法，计算y yn n+1 +1 只用只用y y n n ，而中点法公式为多步法（还可如上二所述，构造而中点法公式为多步法（还可如上二所述，构造 多步法）即必须已知多步法）即必须已知y yn n-1 -1， ，y y n n 才才 能计算能计算y yn n+1 +1，（求 ，（求y y 0 0 , ,y y1 1 不能不能 用此公式。用此公式。y y 0 0 , , y y1 1 称为多步法的开始值，称为多步法的开始值，y y 0 0 给定，给定， 而而y y 1 1 必须必须 由其它公式算出，然后才能用中点法）；由其它公式算出，然后才能用中点法）； 注注3 3：前面已有前面已有EulerEuler法法 的局部截断误差：的局部截断误差： 后退后退EulerEuler法的局部截断误差：法的局部截断误差： 误差阶误差阶：如果局部截断误差如果局部截断误差 则称方法为则称方法为P P阶的。阶的。 14阜师院数科院第八章 常微分方程数值解法 W Y 显然，步长显然，步长h h越小，阶数越小，阶数P P越高，局部截断误差越小，当越高，局部截断误差越小，当 然计算精度越高；然计算精度越高； 注注4 4：梯形法是几阶？梯形法精度比梯形法是几阶？梯形法精度比EulerEuler法高，阶数肯定法高，阶数肯定 比比EulerEuler法高，其实我们可以利用数值积分公式的误差估法高，其实我们可以利用数值积分公式的误差估 计式，因为我们是用梯形数值积计式，因为我们是用梯形数值积 分公式计算分公式计算 因此由积分中梯形公式的误差知此因此由积分中梯形公式的误差知此 时的局部截断误差为：时的局部截断误差为： 梯形法为梯形法为2 2阶方法阶方法 ！ EulerEuler法，后退法，后退EulerEuler法为法为1 1阶方法阶方法，而中点法为而中点法为2 2 阶阶， Euler公式注释（续） 15阜师院数科院第八章 常微分方程数值解法 W Y 关于Euler法的整体截断误差注释 注注5 5：关于关于EulerEuler法法 的整体截断误差：的整体截断误差： EulerEuler方法的局部截断误差公式为：方法的局部截断误差公式为： 实际计算时，实际计算时，y y n n 是是y y( (x x n 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），局部截断），局部截断 误差有界，即存在误差有界，即存在MM00， 使得对任意使得对任意x x a a, ,b b ，都有都有 | |y y( (x x)|)|MM，从而有：从而有： （紧接下屏）紧接下屏） 16阜师院数科院第八章 常微分方程数值解法 W Y 式（式（8-78-7）表明，）表明，EulerEuler方法的整体截断误差与方法的整体截断误差与h h同阶，当同阶，当 h h0 0时，时，e e n n 0 0。 关于Euler法的整体截断误差注释（续 ） 反复递推得：反复递推得： 17阜师院数科院第八章 常微分方程数值解法 W Y 结 论 对于实际问题来说，由于对于实际问题来说，由于L L，M M 难以估计，因此（难以估计，因此（8-68-6 ） 很难应用，而且上述推导过程中一再放大了误差上限，这很难应用，而且上述推导过程中一再放大了误差上限，这 样的估计往往也很保守，远远大于实际的误差，但是，从样的估计往往也很保守，远远大于实际的误差，但是，从 估计式（估计式（8-78-7）却可以得到下面很有用处的结论。）却可以得到下面很有用处的结论。 1 1）当当h h0 0时，时，e e n n 0 0即，即， 亦即数值解亦即数值解y y n n ，一致收敛于初值问题（ 一致收敛于初值问题（8-18-1）的真解）的真解y y( (x x n n ) ) ， 并且，并且，EulerEuler法的整体截断误差的阶为法的整体截断误差的阶为O O( (h h) )与与h h同阶，比局同阶，比局 部截断误差低一阶。部截断误差低一阶。 2 2）舍入误差舍入误差 局部截断误差局部截断误差 对实际计算结果有影响，并且随对实际计算结果有影响，并且随h h减少减少 而减少或增大。而减少或增大。 18阜师院数科院第八章 常微分方程数值解法 W Y 3 3）计算结果与解法的阶数计算结果与解法的阶数p p，真解的导数真解的导数y y ( (p p+1) +1)有 有 关，关，p p越大，越大，h h p+ p+1 1越小， 越小， | |y y( (p p+1) +1)( ( )| )|的上限越大， 的上限越大，M M 也越大，因此为保证精度当然应选阶数也越大，因此为保证精度当然应选阶数p p较高的方较高的方 法。但如果法。但如果M M 很大，当很大，当f f ( (x x, ,y y) )是分段连续的函数时是分段连续的函数时 ，则应采用低阶的方法如用，则应采用低阶的方法如用EulerEuler法。法。 结 论（续） 4 4）计算结果还与开始值的精度有关，为使这种计算结果还与开始值的精度有关，为使这种 误差的影响不致于超过局部误差的影响不致于超过局部 截断误差，对多步截断误差，对多步 法，应采用跟多步法同阶的方法计算开始值。法，应采用跟多步法同阶的方法计算开始值。 19阜师院数科院第八章 常微分方程数值解法 W Y 1.2 改进的Euler法 梯形公式为二阶方法，但却是隐式格式，即若利用梯梯形公式为二阶方法，但却是隐式格式，即若利用梯 形公式求形公式求y yn n+1 +1， ，就要求解方程（就要求解方程（8-58-5）式，计算量较大，通）式，计算量较大，通 常在实际计算时，将常在实际计算时，将EulerEuler法法与与梯形公式梯形公式合起来使用，即合起来使用，即 先使用先使用EulerEuler公式公式, ,由（由（x x n 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 y n n = = y y( (x x n n ) ) ，y y( (x xn n+1 +1) ) 的 的 TaylorTaylor展式为展式为: : 对于改进的对于改进的EulerEuler公式，由于公式，由于 这说明改进的这说明改进的EulerEuler法的局部法的局部 截断误差为截断误差为O O( (h h 3 3 ) )，比，比EulerEuler公式高一阶，是二阶方法。公式高一阶，是二阶方法。 23阜师院数科院第八章 常微分方程数值解法 W Y 改进的Euler公式举例 例例1 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.406264 1 1 .367879.367879.348679.348679.385543.385543.374310.374310.367573.367573 25阜师院数科院第八章 常微分方程数值解法 W Y 表格8-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.367573 0.01.366033.369711.367944.367877 0.001.367700.368052.367879.367876 0.0001.367800.367800.367881.368020 （紧接下屏）紧接下屏） 26阜师院数科院第八章 常微分方程数值解法 W Y 表8-2计算结果说明（续） 但但h h太小，到太小，到h h=0.0001=0.0001时却又变得误差大了，这时却又变得误差大了，这 与前面所说与前面所说h h越小，越小，p p阶越高，应该局部截断误差阶越高，应该局部截断误差 越小，因而计算精度更高矛盾了，为什么会产生这越小，因而计算精度更高矛盾了，为什么会产生这 种情况呢？种情况呢？ 这是由于这是由于h h太小而引起计算量大因而造成了舍入太小而引起计算量大因而造成了舍入 误差和截断误差的积累，这种情况由于初值问题不误差和截断误差的积累，这种情况由于初值问题不 同可能会影响更大，偏离更严重，同可能会影响更大，偏离更严重，如下面的例如下面的例2 2 。 这种问题实际上是稳定性问题，我们将会讨论方法这种问题实际上是稳定性问题，我们将会讨论方法 的稳定性，由此得出对的稳定性，由此得出对h h有一定的要求的稳定性制有一定的要求的稳定性制 区域。区域。 27阜师院数科院第八章 常微分方程数值解法 W Y 例例2 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 .1 1 1 .169351.169351E-4E-4.514229.514229E+6E+6 0 0 .01.01.203704.203704E-9E-9.120746.120746E-7E-7.413244.413244E+7E+7.192743.192743E-8E-8 .001.001.168300.168300E-8E-8.251090.251090E-8E-8.484136.484136E+5E+5.205979.205979E-8E-8 .0001.0001.202081.202081E-8E-8.210331.210331E-8E-8.475130+3.475130+3.206103.206103E-8E-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 x i i , ,x x i i+1+1 上多预报几个点的斜率值，然后将它 上多预报几个点的斜率值，然后将它 们们 加权平均作为加权平均作为k*k*的近似值，则有可能构造出更高精度的计的近似值，则有可能构造出更高精度的计 算公式，这是算公式，这是龙格龙格- -库塔库塔方法的基本思想。方法的基本思想。 用这个观点来研究欧拉公式与改进欧拉公式，可以发现用这个观点来研究欧拉公式与改进欧拉公式，可以发现 欧拉公式由于仅取欧拉公式由于仅取x x n n 一个点的斜率值一个点的斜率值f f ( (x x n n, ,y yn n ) )作为平均斜率作为平均斜率 k k* * 的近似值，因此精度很低。而改进欧拉公式（的近似值，因此精度很低。而改进欧拉公式（8-108-10）却）却 是利用了是利用了x x n n 与与x xn n-1 -1两个点的斜率值 两个点的斜率值k k 1 1 = = f f ( (x x n n, ,y yn n ) )与与 k k2 2 = =f f ( (x xn n+1 +1 , ,y yn n + +hkhk 1 1 ) )取算术平均作为平均斜率取算术平均作为平均斜率k k* *的近似值。的近似值。 其中其中k k 2 2 是通过已知信息是通过已知信息y y n n 利利 用欧拉公式求得的。用欧拉公式求得的。 30阜师院数科院第八章 常微分方程数值解法 W Y 2.2 二阶龙格-库塔公式 首先推广改进欧拉公式，首先推广改进欧拉公式， 考察区间考察区间 x x n n, ,x x n n+1+1 内任一点： 内任一点： 我们希望用我们希望用x x n n 和和x x n n+1 +1两个点的斜率值 两个点的斜率值k k 1 1 和和k k 2 2 加权平均作为加权平均作为 平均斜率平均斜率k k* *的近似值的近似值 ： 其中其中c c 1 1 ，c c 2 2 为待定常数，同为待定常数，同 改进欧拉公式一样，这里仍取：改进欧拉公式一样，这里仍取： 问题在于怎样预测问题在于怎样预测x xn n+ +l l处的斜率值处的斜率值k k 2 2 。 仿照改进欧拉公式，先用仿照改进欧拉公式，先用 欧拉公式提供欧拉公式提供y y ( (x xn n+ +l l) )的预测值的预测值 然后再用预测值然后再用预测值y yn n+ +l l通过计算通过计算f f 产生斜率值产生斜率值k k 2 2 = =f f ( (x xn n+ +l l , ,y y n n+ +l l) )， ， 这样设计出的计算公式具有形式这样设计出的计算公式具有形式 ： 31阜师院数科院第八章 常微分方程数值解法 W Y 二阶龙格-库塔公式（续1） 公式（公式（8-128-12）中含有三个待定参数）中含有三个待定参数c c 1 1, ,c c2 2 和和l l，我们希望适我们希望适 当选取这些参数值，使得公式（当选取这些参数值，使得公式（8-128-12）具有二阶精度，亦）具有二阶精度，亦 即使：即使： 现在仍假定现在仍假定y y n n = =y y( (x x n n ) )，即，即y y n n 是准确的，将是准确的，将y y( (x xn n+1 +1) )与 与y yn n+1 +1 都在都在x x i i 处作泰勒展开：处作泰勒展开： 32阜师院数科院第八章 常微分方程数值解法 W Y 二阶龙格-库塔公式（续2） 代入（代入（8-128-12）式，得：）式，得： 比较（比较（8-138-13）与（）与（8-148-14）两式，要使公式具有二阶精度，）两式，要使公式具有二阶精度， 只有满足下列条件：只有满足下列条件： 这里一共有三个待定参数，但只需满足两个条件，因这里一共有三个待定参数，但只需满足两个条件，因 此有一个自由度，于是满足条件（此有一个自由度，于是满足条件（8-158-15）的参数不止一组）的参数不止一组 ，而是一族，相应的公式（，而是一族，相应的公式（8-128-12）也有一族，这些公式统）也有一族，这些公式统 称为称为二二阶龙格阶龙格- -库塔公式库塔公式，简称，简称二阶二阶R R- -K K公式公式 。 33阜师院数科院第八章 常微分方程数值解法 W Y 特别，当特别，当l l=1=1即即x xn n+l +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 c 1 1 =0,=0,c c 2 2 =1=1， 这时二阶这时二阶R R- -K K公式称为公式称为变形变形 的欧拉公式的欧拉公式，其形式见左边，其形式见左边 ： 从表面上看，变形的欧拉公式仅含一个斜率值从表面上看，变形的欧拉公式仅含一个斜率值 k k2 2 ，但，但k k 2 2 是是通过通过k k 1 1 计算出来的，因此每完成一步，计算出来的，因此每完成一步， 仍然需要两次计算函数仍然需要两次计算函数f f 的值，工作量和改进欧拉的值，工作量和改进欧拉 公式相同。公式相同。 二阶龙格-库塔公式（续3） 34阜师院数科院第八章 常微分方程数值解法 W Y 构造二阶R-K公式的主要步骤 综上所述，构造二阶综上所述，构造二阶R-KR-K公式主要由以下几步产生公式主要由以下几步产生 ： 1)1) 在区间在区间 x x n n, ,x x n n+1+1 上取二点，预报相应点 上取二点，预报相应点 的斜率值；的斜率值； 2) 2) 对此两斜率值加权平均作为平均斜率对此两斜率值加权平均作为平均斜率 值的近似值；值的近似值； 3) 3) 将将y yn n+1 +1与 与y y( (x xn n+1 +1) )都在 都在x x i i 处作泰勒展开，处作泰勒展开， 为使公式达到二阶精度，比较相应系为使公式达到二阶精度，比较相应系 数，建立有关参数所应满足的方程组；数，建立有关参数所应满足的方程组； 4) 4) 解此方程组得一族二阶解此方程组得一族二阶R-KR-K公式。公式。 35阜师院数科院第八章 常微分方程数值解法 W Y 2.3 高阶R-K公式 为了进一步提高精度，在为了进一步提高精度，在 x x n n, ,x x n n+1+1 上除 上除x x n n 和和x xn n+ +l l外再增加外再增加 一点一点x xn n+m +m = =x x n n + +mhmh( (l l m m 1) 1)，并用并用x x n n, ,x x n n+l+l , ,x x n n+ +mm三处的斜率值 三处的斜率值 k k1 1, ,k k2 2, ,k k3 3 加权平均作为加权平均作为k k* *的近似值，这时计算公式为：的近似值，这时计算公式为： 其中其中k k 1 1, ,k k2 2 仍用（仍用（8-128-12）式所取的形式。）式所取的形式。 为了预测为了预测x xn n+ +m m处的斜率值 处的斜率值k k 3 3 ，要定出要定出x xn n+ +m m处所对应的 处所对应的 y y n n+ +mm， ，可以看作在区间可以看作在区间 x x n n, ,x x n n+ +mm 上使用二阶 上使用二阶R-KR-K公式，从而公式，从而 得到得到y y( (x xn n+ +m m) )的预测值： 的预测值： 于是，再通过计算函数值于是，再通过计算函数值f f 得到得到: : 36阜师院数科院第八章 常微分方程数值解法 W Y 高阶R-K公式（续1） 这样设计出的计算公式具有形式这样设计出的计算公式具有形式: : 运用泰勒展开方法，适当选择参数运用泰勒展开方法，适当选择参数c c 1 1, ,c c2 2, ,c c3 3 , ,l l, ,mm, ,及及 b b1 1, ,b b2 2 使上述公式具有三阶精度采用与使上述公式具有三阶精度采用与(8-12)(8-12)类似类似 的处理方法，得到这些参数需要满足的条件：的处理方法，得到这些参数需要满足的条件： ( (见下屏见下屏) ) 37阜师院数科院第八章 常微分方程数值解法 W Y 高阶R-K公式（续2） 满足满足5 5个条件的一族公式（个条件的一族公式（8-188-18） ，统称为三阶，统称为三阶R-KR-K公式，其中常公式，其中常 见的是库塔公式：见的是库塔公式： c c 1 1 =1/6 , c=1/6 , c 2 2 =4/6 , c=4/6 , c 3 3 =1/6 , =1/6 , l=1/2 , m=1,b l=1/2 , m=1,b 1 1 =-1 , b=-1 , b 2 2 =2=2 38阜师院数科院第八章 常微分方程数值解法 W Y 若需再将精度提高至四阶，用上述类似处理方法，只若需再将精度提高至四阶，用上述类似处理方法，只 是必须在是必须在 x x i i , ,x x i i+1+1 上用四个点处的斜率加权平均作为 上用四个点处的斜率加权平均作为k k* *的近的近 似值，构成一族四阶龙格似值，构成一族四阶龙格- -库塔公式，由于推导复杂，这库塔公式，由于推导复杂，这 里从略，只将常用的公式介绍如下：里从略，只将常用的公式介绍如下： 四阶公式中常用的还四阶公式中常用的还 有下面的有下面的GillGill公式：公式： 高阶R-K公式（续3） 一般称它为标准四阶一般称它为标准四阶 R-KR-K公式，其局部截断误公式，其局部截断误 差是差是O O( (h h 5 5 ) ) 。 39阜师院数科院第八章 常微分方程数值解法 W Y Gill 公式 在计算机上它具有节省存贮单元和控制舍入在计算机上它具有节省存贮单元和控制舍入 误差增长的优点。误差增长的优点。 40阜师院数科院第八章 常微分方程数值解法 W Y 初值问题两种方法举例 用改进的用改进的EulerEuler 法法和四阶标准和四阶标准 R-KR-K公式解初值公式解初值 问题：问题： 例例3 3 解解 用改进的用改进的EulerEuler公式，取公式，取h h=0.1=0.1计算公式为：计算公式为： 部分计算结果见表部分计算结果见表8-48-4： 表表8-48-4 x xn n 改进的欧改进的欧 拉法拉法 准准 确确 解解 y y( (x x n n ) ) 误差误差 x xn n 改进的欧改进的欧 拉法拉法 准准 确确 解解 y y( (x x n n ) ) 误差误差 0 0 1. 1. 1 1 0 0 0.60.61.4859561.485956 1.4832401.4832400.0027160.002716 0.20.21.1840961.184096 1.1832161.1832160.0000880.000088 0.80.81.6164761.616476 1.6124521.6124520.0040240.004024 0.40.41.3433601.343360 1.3416411.3416410.0017190.001719 1.01.01.7378691.737869 1.7320511.7320510.0058180.005818 41阜师院数科院第八章 常微分方程数值解法 W Y 对四阶标准对四阶标准R-KR-K法法 ， 取取h h=0.2=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 0 1 1 1 1 0.60.61.4832811.4832810.0000410.000041 0.20.2 1.1832291.1832290.0000130.0000130.80.81.6125131.6125130.0000610.000061 0.40.4 1.3416671.3416670.000260.000261.01.01.7321401.7321400.0000890.000089 计算结果见计算结果见表表8-58-5，与准确解比较，至少有四位有效，与准确解比较，至少有四位有效 数字。而取步长数字。而取步长h h=0.1=0.1时，用改进的欧拉公式计算，最多时，用改进的欧拉公式计算，最多 只有三位有效数字只有三位有效数字（见表（见表8-48-4），虽然改进的欧拉公式每，虽然改进的欧拉公式每 前进一步只要计算两次前进一步只要计算两次f f 值，而四阶值，而四阶R R- -K K法要计算法要计算4 4次次f f 值值 ， 但由于后者步长比前者大，所耗费的计算总量还是差不多但由于后者步长比前者大，所耗费的计算总量还是差不多 的。这个例子又一次显示了选择算法的重要意义。的。这个例子又一次显示了选择算法的重要意义。 43阜师院数科院第八章 常微分方程数值解法 W Y 表表8-68-6计算计算f f 的次数与精度阶数的关系的次数与精度阶数的关系 每每 步步 计计 算算 f f 的的 次次 数数 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 精度的阶精度的阶 数数 2 2 3 3 4 4 4 4 5 5 6 6 6 6 7 7 表表8-6 8-6 的说明的说明 一般，高精度的方法要求解有较好的光滑性，解的光滑一般，高精度的方法要求解有较好的光滑性，解的光滑 性是由性是由f (x,y)f (x,y)的可微性决定。如果解的光滑性差，则用四的可微性决定。如果解的光滑性差，则用四 阶阶R-KR-K法求得的数值解反而不如用改进的欧拉法好。因此法求得的数值解反而不如用改进的欧拉法好。因此 ，一定要针对具体问题的特点来选择求解方法。，一定要针对具体问题的特点来选择求解方法。 例 3 （续3） 44阜师院数科院第八章 常微分方程数值解法 W Y 表8-6 的说明（续） 龙格龙格- -库塔方法的推导是在用泰勒展开方法的库塔方法的推导是在用泰勒展开方法的 基础上进行的，因而它要求所求的解具有较好的基础上进行的，因而它要求所求的解具有较好的 光滑性质。假若解的光滑性差，那么使用四阶光滑性质。假若解的光滑性差，那么使用四阶R-R- K K公