第九章.常微分方程初值问题数值解法_文档资料库
当前位置： >>
第九章.常微分方程初值问题数值解法
第九章 常微分方程初值问题数值解法9.1 引言科学技术中常常需要求解常微分方程的定解问题．这类问题最简单的形式， 是本章将要着重考察的一阶方程的初值问题? y ' = f ( x, y ) ? ? ? y ( x0 ) = y0 ?(1.1) (1.2 )我们只要，
只要函数 f ( x, y ) 适当光滑―譬如关于 y 满足利普希茨 ( Lipschitz ) 条件f ( x, y ) ? f x, y ≤ L y ? y .理论上就可以保证初值问题 (1.1) ， (1.2 ) 的解 y = y ( x ) 存在并且唯一.( )(1.3)虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些 特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法. 所谓数值解法 数值解法，就是寻求解 y ( x ) 在一系列离散节点 数值解法 x1 & x2 & L & xn & xn +1 & L 上的近似值 y1 , y2 ,L yn , yn +1 ,L .相邻两个节点的间距 hn = xn +1 ? xn 成为步长 步长.今后如不特别说明，总是假定 hi = h ( i = 1, 2,L) 为定数，这时节点为 步长 xn = x0 + nh, n = 0,1, 2, L. 初值问题 (1.1) ， (1.2 ) 的数值解法有个基本特点，它们都采取“步进式”，即 求解过程顺着节点排列的次序一步一步地向前推进．描述这类算法，只要给出用 已知信息 yn , yn ?1 , yn ? 2 ,L 计算 yn +1 的递推公式. 首先，要对方程 (1.1) 离散化，建立求数值解的递推公式．一类是计算 yn +1 时 只 用 到 前 一 点 的 值 yn ， 称 为 单 步 法 ． 另 一 类 是 用 到 yn +1 前 面 k 点 的 值yn , yn ?1 ,L yn ? k +1 , 称为 k 步法． 其次， 要研究公式的局部截断误差和阶， 数值解 yn 与 精确解 y ( xn ) 的误差估计及收敛性，还有递推公式的计算稳定性等问题．9.2 9.2.1简单的数值方法与基本概念欧拉法与后退欧拉法我们知道，在 xy 平面上，微分方程 (1.1) 的解 y = y ( x ) 称为它的积分曲线．积 分曲线上一点 ( x, y ) 的切线斜率等于函数 f ( x, y ) 的值． 如果按函数 f ( x, y ) 在 xy 平 面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点方 向相一致． 基于上述几何解释，我们从初始点 p0 ( x0 , y0 ) 出发，先依方向场在该点的方向 推进到 x = x1 上一点 p1 ，然后再从 p1 依方向场的方向推进到 x = x2 上一点 p2 ，循 此前进做出一条折线 p0 p1 p2 L （图 9 ? 1 ）． 一般地，设已做出该折线的顶点 pn ，过 pn ( xn , yn ) 依方向场的方向再推进到pn +1 ( xn +1 , yn +1 ) ，显然两个定点 pn ， pn +1 的坐标有关系yn +1 ? yn = f ( xn , yn ) , xn +1 _ xn即yn +1 = yn + hf ( xn , yn )( 2.1)这就是著名的欧拉 ( Euler ) 公式．若初值 y0 已知，则依公式 ( 2.1) 可逐步算出y1 = y0 + hf ( x0 , y0 ) , y2 = y1 + hf ( x1 , y1 ) , L例1求解初值问题2x ? ' ?y = y ? y ( 0 & x & 1) ， ? ? y ( 0) = 1 ?( 2.2 )解 为便于进行比较，本章将用多种数值方法求解上述初值问题．这里先用 欧拉方法，欧拉公式的具体形式为? 2x ? yn +1 = yn + h ? yn ? n ? . yn ? ?取步长 h = 0.1, 计算结果见表 9 ? 1 表 9 ? 1 计算结果对比 xn0.1 0.2 0.3 0.4 0.5yn1.8 1.2 1.4351y ( xn )1.2 1.6 1.4142xn0.6 0.7 0.8 0.9 1.0yn1.3 1.8 1.7848y ( xn )1.2 1.3 1.7321初值问题 ( 2.2 ) 有解 y = 1 + 2 x ，按这个解析式子算出的准确值 y ( xn ) 同近似 值 yn 一起列在表 9 ? 1 中，两者相比较可以看出欧拉方法的精度很差． 还可以通过几何直观来考察欧拉方法的精度．假设 yn = y ( xn ) ，即顶点 pn 落 在积分曲线 y = y ( x ) 上，那么，按欧拉方法做出的切线 pn pn +1 便是 y = y ( x ) 过点pn 的切线．从图形上看，这样定出顶点 pn +1 显著地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的． 为了分析计算公式的精度，通常可采用泰勒展开将 y ( xn +1 ) 在 xn 处展开，则有y ( xn +1 ) = y ( xn + h )= y ( xn ) + y ' ( xn ) h +h 2 '' y (ξ n ) , 2ξ n ∈ ( xn , xn+1 ) .在 yn = y ( xn ) 的前提下， f ( xn , yn ) = f ( xn , y ( xn ) ) = y ' ( xn ) . 于是可得欧拉法 ( 2.1) 的 公式误差y ( xn +1 ) ? yn +1 = 称为此方法的局部截断误差．h 2 '' h 2 '' y (ξn ) ≈ y ( xn ) ， 2 2( 2.3)如果对方程 (1.1) 从 xn 到 xn +1 积分，得 y ( xn +1 ) = y ( xn ) + ∫ xn +1 f ( t , y ( t ) ) dt ． xn( 2.4 )右端积分用左矩形公式 hf ( xn , y ( xn ) ) 近似， 再以 yn 代替 y ( xn ) ，yn +1 代替 y ( xn +1 ) 也 得到 ( 2.1) ，局部截断误差也是 ( 2.3) ． 如果在 ( 2.4 ) 中右端积分用右矩形公式 hf ( xn +1 , y ( xn +1 ) ) 近似，则得另一个式yn +1 = yn + hf ( xn +1 , y ( xn +1 ) )称为后退得欧拉法． 后退得欧拉法． 后退得欧拉法( 2.5)后退的欧拉公式与欧拉公式有着本质的区别，后者是关于 yn +1 的一个直接的 计算公式，这类公式称作是显式的 显式的；然而公式 ( 2.5 ) 的右端含有未知的 yn +1 ，它实 显式的 际上是关于 yn +1 的一个函数方程，这类公式称作是隐式的 隐式的． 隐式的 显式与隐式两类方法各有特点．考虑到数值稳定性等其他因素，人们有时需 要选用隐式方法，但使用显式算法远比隐式方便． 隐式方程 ( 2.5 ) 通常用迭代法求解，而迭代过程的实质是逐步显示化． 设用欧拉公式( yn0 )1 = yn + hf ( xn , yn ) +( 给出迭代初值 yn0 )1 ，用它代入 ( 2.5 ) 式的右端，使之转化为显示，直接计算得 +( ( yn1+)1 = yn + hf xn +1 , yn0)1 ， +()然后再用 y1 +1 代入 ( 2.5 ) 式，又有 n( ( yn +)1 = yn + hf xn +1 , yn +)1 ．2 1()如此反复进行，得( ( yn +1 ) = yn + hf xn +1 , yn +)1k +1 k()k( k = 0,1,L) ．( 2.6 )由于 f ( x, y ) 对 y 满足利普希茨 ( Lipschitz ) 条件 (1.3) ．由 ( 2.6 ) 减 ( 2.5 ) 得( ( yn +1 ) ? yn +1 = h f xn +1 , yn +)1 ? f ( xn +1 , yn +1 )k +1()( ≤ hL ynk )1 ? yn +1 . +由此可知，只要 hL & 1 迭代法 ( 2.6 ) 就收敛到解 yn +1 ．关于后退欧拉法的公式 误差，从积分公式看到它与欧拉法是相似的．9.2.2梯形方法为得到比欧拉法精度高的计算公式， 在等式 ( 2.4 ) 右端积分中若用梯形求积公 式近似，若用 yn 代替 y ( xn ) ， yn +1 代替 y ( xn +1 ) ，则得h yn +1 = yn + ? f ( xn , yn ) + f ( xn +1 , yn +1 ) ? ? 2?( 2.7 )称为梯形方法． 梯形方法 梯形方法 梯形方法是隐式单步法，可用迭代法求解．同后退的欧拉法一样，仍用欧拉 方法提供迭代初值，则梯形法的迭代公式为( ? yn0 )1 = yn + hf ( xn , yn ) ; + ? h ? ( k +1) (k ) ? yn +1 = yn + ? f ( xn , yn ) + f xn +1 , yn +1 ? ; ? ? 2 ? ?( k = 0,1, 2, L) . ?()( 2.8 )为了分析迭代过程的收敛性，将 ( 2.7 ) 式与 ( 2.8 ) 式相减，得( yn +1 ? yn +1 ) =k +1h? (k f x ,y ? f xn +1 , yn +)1 ? ， ? ( n +1 n +1 ) ? 2k +1()于是有( yn +1 ? yn +1 ) ≤hL (k yn +1 ? yn +)1 , 2式中 L 为 f ( x, y ) 关于 y 的利普希茨常数．如果选取 h 充分小，使得hL ≤ 1， 2( 则当 k → ∞ 时有 yn +)1 → yn +1 ，这说明迭代过程 ( 2.8 ) 是收敛的．k9.2.3单步法的局部截断误差与阶初值问题 (1.1) , (1.2 ) 的单步法可用一般形式表示为yn +1 = yn + h? ( xn , yn , yn +1 , h ) ,( 2.9 )其中多元函数 ? 与 f ( x, y ) 有关，当 ? 含有 yn +1 时，方法是隐式的，若不含 yn +1 则 为显式方法，所以显式单步法可表示为yn +1 = yn + h? ( xn , yn , h ) ,( 2.10 )? ( x, y, h ) 称为增量函数，例如对欧拉法 ( 2.1) 有 ? ( x , y , h ) = f ( x, y ) ．它的局部截断误差已由 ( 2.3) 给出，对一般显式单步法则可如下定义． 定义 1 设 y ( x ) 是初值问题 (1.1) , (1.2 ) 的准确解，称Tn +1 = y ( xn +1 ) ? y ( xn ) ? h? ( xn , y ( xn ) , h )为显式单步法 ( 2.10 ) 的局部截断误差． 局部截断误差． 局部截断误差( 2.11)Tn +1 之所以称为局部的，是假设在 xn 前各步没有误差．当 yn = y ( xn ) 时，计算 一步，则有 y ( xn +1 ) ? yn +1 = y ( xn +1 ) ? ? yn + h? ( xn , yn , h ) ? ? ?= y ( xn +1 ) ? y ( xn ) ? h? ( xn , y ( xn ) , h ) = Tn +1 ．所以，局部截断误差可理解为用方法 ( 2.10 ) 计算一步的误差，也即公式 ( 2.10 ) 中 用准确解 y ( x ) 代替数值解产生的公式误差．根据定义，显然欧拉法的局部截断误 差Tn +1 = y ( xn +1 ) ? y ( xn ) ? h? ( xn , y ( xn ) ) = y ( xn +1 ) ? y ( xn ) ? hy ' ( xn )= 即为 ( 2.3) 的结果．这里 般情形的定义如下， 定义 2h 2 '' y ( xn ) + Ο ( h3 ) , 2h 2 '' y ( xn ) 称为局部截断误差主项．显然 Tn +1 = Ο ( h 2 ) ．一 2设 y ( x ) 是初值问题 (1.1) , (1.2 ) 的准确解，若存在最大整数 p 使显式单步法 ( 2.10 ) 的局部截断误差满足Tn +1 = y ( x + h ) ? y ( x ) ? h? ( x, y, h ) = Ο ( h p +1 ) ，( 2.12 )则称方法 ( 2.10 ) 具有 p 阶精度 阶精度． 若将 ( 2.12 ) 展开式写成Tn +1 = ψ ( xn , y ( xn ) ) h p +1 + Ο ( h p + 2 ) ，则ψ ( xn , y ( xn ) ) h p +1 称为局部截断误差主项 局部截断误差主项． 局部截断误差主项 以上定义对隐式单步法 ( 2.9 ) 也是适用的．例如，对后退欧拉法 ( 2.5 ) 其局部 截断误差为Tn +1 = y ( xn +1 ) ? y ( xn ) ? hf ( xn +1 , y ( xn +1 ) )= hy ' ( xn ) +h 2 '' y ( xn ) + Ο ( h3 ) 2? h ? y ' ( xn ) + hy '' ( xn ) + Ο ( h 2 ) ? ? ? =?h 2 '' y ( xn ) + Ο ( h 3 ) . 2这里 p = 1 ，是 1 阶方法，局部截断误差主项为 ? 同样对梯形法 ( 2.7 ) 有h 2 '' y ( xn ) ． 2h Tn +1 = y ( xn +1 ) ? y ( xn ) ? ? y ' ( xn ) + y ' ( xn +1 ) ? ? 2? = hy ' ( xn ) + h 2 '' h3 y ( xn ) + y ''' ( xn ) 2 3!? h? h2 ? ? y ' ( xn ) + y ' ( xn ) + hy '' ( xn ) + y ''' ( xn ) ? + Ο ( h4 ) 2? 2 ?=? h3 ''' y ( xn ) + Ο ( h 4 ) . 12 h3 ''' y ( xn ) ． 12所以梯形方法 ( 2.7 ) 是二阶的，其局部误差主项为 ?9.2.4改进的欧拉公式我们看到， 梯形方法虽然提高了精度， 但其算法复杂， 在应用迭代公式 ( 2.9 ) 进行实际计算时，每迭代一次，都要重新计算函数 f ( x, y ) 的值，而迭代又要反复 进行若干次，计算量很大，而且往往难以预测．为了控制计算量，通常只迭代一 两次就转入下一步的计算，这就简化了算法． 具体地说，我们先用欧拉公式求得一个初步的近似值 y n +1 ，称之为预测值， 预测值， 预测值 预测值 y n +1 的精度可能很差， 再用梯形公式 ( 2.7 ) 将它校正一次， 即按 ( 2.8 ) 式迭代 一次得 yn +1 ，这个结果称校正值 校正值，而这样建立的预测-校正系统通常称为改进的欧 校正值 改进的欧 拉公式： 拉公式： 预测 校正y n +1 = yn + hf ( xn , yn ) ,yn +1 = yn +h? f ( xn , yn ) + f xn +1 , y n +1 ? . ? 2?()( 2.13)或表为下列平均化形式? ? y p = yn + hf ( xn , yn ) , ? ? ? yc = yn + hf ( xn +1 , y p ) , ? ? y = 1 ( y + y ). c ? n +1 2 p ?例 2 用改进的欧拉方法求解初值问题 ( 2.2 ) ． 解 改进的欧拉公式为? ? 2 xn ? ? y p = yn + hf ? yn ? ?, yn ? ? ? ? ? 2 xn +1 ? ? ?, ? yc = yn + hf ? y p ? ? ? yp ? ? ? ? 1 ? yn +1 = ( y p + yc ) . 2 ? ? 仍取 h = 0.1 ，计算结果见表 9 ? 2 ．同例 1 中欧拉法的计算结果比较，改进欧拉法 明显改善了精度． 表9 ? 2 计算结果对比xnyny ( xn )1.2 1.6 1.xnyny ( xn )1.2 1.3 1.73210.1 0.2 0.3 0.4 0.51.1 1.4 1.41640.6 0.7 0.8 0.9 1.01.5 1.2 1.7379龙格-库塔方法 龙格 库塔方法9.3.1显式龙格-库塔法的一般形式 显式龙格 库塔法的一般形式上节给出了显式单步法的表达式 ( 2.10 ) ，其局部截断误差为 ( 2.12 ) ，对欧拉 法 Tn +1 = Ο ( h 2 ) ，即方法为 p = 1 阶，若用改进欧拉法 ( 2.13) ，它可表为h yn +1 = yn + ? f ( xn , yn ) + f ( xn + h, yn + hf ( xn , yn ) ) ? ． ? 2?( 3.1) ( 3.2 )此时增量函数? ( xn , yn , h ) = ? f ( xn , yn ) + f ( xn + h, yn + hf ( xn , yn ) ) ? ． ? ?1 2它比欧拉法的 ? ( xn , yn , h ) =f ( xn , yn ) ，增加了计算一个右函数 f 的值，可望p = 2 ．若要使得到的公式阶数 p 更大，? 就必须包含更多的 f 值．实际上从方程(1.1) 等价的积分形式 ( 2.4 ) ，即y ( xn +1 ) ? y ( xn ) = ∫xn +1xnf ( x, y ( x ) )dx ，( 3.3)若要使公式阶数提高，就必须使右端积分得数值求积公式精度提高，它必然要增 加求积节点，为此可将 ( 3.3) 的右端用积分公式表示为∫xn+1xnf ( x, y ( x ) )dx ≈ ∑ ci f ( xn + λi h, y ( xn + λi h ) ) .i =1r一般来说，点数 r 越多，精度越高，上式右端相当于增量函数 ? ( x, y, h ) ，为得到 便于计算的显示方法，可类似于改进欧拉法 ( 3.1) , ( 3.2 ) ，将公式表示为yn +1 = yn + h? ( xn , yn , h ) ,其中( 3.4 )? ( xn , yn , h ) = ∑ ci K i ，i =1r( 3.5)K i = f ( xn , yn ) ,i ?1 ? ? K i = f ? xn + λi h, yn + ∑ ?ij K j ? j =1 ? ?i = 2,L , r ,这里 ci ， i ， ij 均为常数． 3.4 ) 和 ( 3.5 ) 称为 r 级显式龙格-库塔 ( Runge ? Kutta ) 法， λ ? ( 简称 R ? K 方法． 当 r = 1 ， ? ( xn , yn , h ) = f ( xn , yn ) 时，就是欧拉法，此时方法的阶为 p = 1 ．当r = 2 时，要改进欧拉法 ( 3.1) ， ( 3.2 ) 就是其中的一种，下面将证明阶 p = 2 ．要使公式 ( 3.4 ) ， 3.5 ) 具有更高的阶 p ， 就要增加点数 r ． 下面我们只就 r = 2 推导 R ? K ( 方法． 并给出 r = 3, 4 时的常用公式， 其推导方法与 r = 2 时类似， 只是计算较复杂．9.3.2二阶显式 R ? K 方法对 r = 2 的 R ? K 方法，由 ( 3.4 ) ， ( 3.5 ) 可得到如下的计算公式? yn +1 = yn + h ( c1 K1 + c2 K 2 ) , ? ? K1 = f ( xn , yn ) , ? ? K 2 = f ( xn + λ2 h, yn + ? 21hK1 ) .( 3.6 )这里 c1 ，c2 ，λ2 ， ?21 均为待定常数，我们希望适当选取这些系数，使公式阶数 p 尽量高．根据局部截断误差定义， ( 3.6 ) 的局部截断误差为Tn +1 = y ( xn +1 ) ? y ( xn )= ? h ?c1 f ( xn , yn ) + c2 f ( xn + λ2 h, yn + ? 21hf n ) ? , ? ?( 3.7 )这里 yn = y ( xn ) , f n = f ( xn , yn ) . 为得到 Tn +1 的阶 p ，要将上式各项在 ( xn , yn ) 处做泰 勒展开，由于 f ( x, y ) 是二元函数，估要用到二元泰勒展开，各项展开式为' y ( xn +1 ) = yn + hyn +h 2 '' h3 ''' yn + yn + Ο ( h 4 ) , 2 3!其中' ? yn = f ( xn , yn ) = f n , ? ? '' d f ( xn , y ( xn ) ) = f x' ( xn , yn ) + f y' ( xn , yn ) ? f n , ? yn = dx ? ''' '' ' '' ? yn = f xx ( xn , yn ) + 2 f n f xy ( xn , yn ) + f n2 f yy ( xn , yn ) + f y' ( xn , yn ) ? f x' ( xn , yn ) + f n f y' ( xn , yn ) ? ; ? ? ?( 3.8)f ( xn + λ2 h, yn + ? 21hf n )= f n + f x' ( xn , yn ) λ2 h + f y' ( xn , yn ) ?21hf n + Ο ( h 2 ) .将以上结果代入 ( 3.7 ) 则有Tn +1 = hf n + h2 ' ? f x ( xn , yn ) + f y' ( xn , yn ) f n ? ? 2 ?? h ?c1 f n + c2 ( f n + λ2 f x' ( xn , yn ) h ) + ? 21 f y' ( xn , yn ) f n h ? + O ( h 3 ) ? ??1 ? = (1 ? c1 ? c2 ) f n h + ? ? c2 λ2 ? f x' ( xn , yn ) h 2 ?2 ? ?1 ? + ? ? c2 ?21 ? f y' ( xn , yn ) f n h 2 + O ( h3 ) ． ?2 ?要使公式 ( 3.6 ) 具有 p = 2 阶，必须使1 1 1 ? c1 ? c2 = 0, ? c2 λ2 = 0, ? c2 ?21 = 0. 2 2( 3.9 )即 1 1 c2 λ2 = , c2 ?21 = , c1 + c2 = 1. 2 2( 3.9 ) 的解是不唯一的．可令 c2 = a ≠ 0 ，则得c1 = 1 ? a, λ2 = ?21 = 1 . 2a这样得到的公式称为二阶 R-K 方法， 如取 a = 1/ 2 ， c1 = c2 = 1 / 2 ， 2 = ?21 = 1 ． 则 λ 这 就是改进欧拉法 ( 3.1) ． 若取 a = 1 ，则 c2 = 1, c1 = 0, λ2 = ?21 = 1/ 2 ．得计算公式? ? y = y + hK , n 2 ? n +1 ? ? K1 = f ( xn , yn ) , ? ? K 2 = f ? xn + h , yn + h K1 ? . ? ? ? 2 2 ? ? ?称为中点公式 中点公式，相当于数值积分公式的中矩形公式． ( 3.10 ) 也可表示为 中点公式h h ? ? yn +1 = yn + hf ? xn + , yn + f ( xn , yn ) ? . 2 2 ? ?( 3.10 )对 r = 2 的 R-K 公式 ( 3.6 ) 能否使局部误差提高到 O ( h 4 ) ？为此需把 K 2 多展开''' 一项，从 ( 3.8 ) 的 yn 看到展开式中 f y' f x' + ( f y' ) f 的项是不能通过选择参数消掉的， 2实际上要使 h3 的项为零，需增加 3 个方程，要确定 4 个参数 c1 , c2 , λ2 及 ?21 ，这是 不可能的．故 r = 2 的显示方法的阶只能是 p = 2 ，而不能得到三阶公式．9.3.3三阶与四阶显式 R ? K 方法要得到三阶显式 R ? K 方法，必须 r = 3 ．此时 ( 3.4 ) , ( 3.5) 的公式表示为? yn +1 = yn + h ( c1 K1 + c2 K 2 + c3 K 3 ) , ? ? K1 = f ( xn , yn ) , ? ? K 2 = f ( xn + λ2 h, yn + ?21hK1 ) , ? K = f ( x + λ h, ? hK + ? hK ) . n 3 31 1 32 2 ? 3( 3.11)其中 c1 , c2 , c3 及 λ2 , ?21 , λ3 , ?31 , ?32 均为待定参数，公式 ( 3.11) 的局部截断误差为Tn +1 = y ( xn +1 ) ? y ( xn ) ? h [ c1 K1 + c2 K 2 + c3 K 3 ] .只要将 K 2 , K 3 按二元函数泰勒展开，使 Tn +1 = O ( h 4 ) ，可得待定参数满足方程? c1 + c2 + c3 = 1, ? λ =? , 2 21 ? ? λ3 = ?31 + ?32 , ? ?c λ + c λ = 1 , ? 2 2 3 3 2 ? ?c2 λ22 + c3λ32 = 1 , ? 3 ? 1 ? c3λ2 ?32 = . 6 ?( 3.12 )这是 8 个未知数 6 个方程的方程组，解也不是唯一的．可以得到很多公式．满足 条件的公式统称为三阶 R ? K 公式．下面只给出其中一个常见的公式．h ? ? yn +1 = yn + 6 ( K1 + 4 K 2 + K 3 ) , ? ? K1 = f ( xn , yn ) , ? ? ? K 2 = f ? xn + h , yn + h K1 ? , ? ? 2 2 ? ? ? ? K 3 = f ( xn + h, yn ? hK1 + 2hK 2 ) .此公式称为库塔 库塔三阶方法． 库塔 继续上述过程，经过较复杂的数学演算，可以导出各种四阶龙格-库塔公式， 下列经典公式是其中常用的一个：? h ? yn +1 = yn + ( K1 + 2 K 2 + 2 K 3 + K 4 ) , 6 ? ? K1 = f ( xn , yn ) , ? h h ? ? ? . ? K 2 = f ? xn + , yn + K1 ? , 2 2 ? ? ? ? h h ? ? ? K 3 = f ? xn + , yn + K 2 ? , 2 2 ? ? ? ? K = f ( x + h, y + hK ) . n n 3 ? 4( 3.13)四阶龙格 - 库塔方法的每一步需要四次函数值 f ，可以证明其截断误差为O ( h5 ) ．不过证明极其繁琐，这里从略．例3 题 ( 2.2 ) ．设取步长 h = 0.2 ，从 x = 0 直到 x = 1 用四阶龙格-库塔方法求解初值问解 这里，经典的四阶龙格-库塔公式 ( 3,13) 具有形式h ? ? yn +1 = yn + 6 ( K1 + 2 K 2 + 2 K 3 + K 4 ) , ? 2x ? K1 = yn ? n , ? yn ? 2 xn + h h ? , K 2 = yn + K1 ? ? h 2 yn + K1 ? 2 ? 2 xn + h h ? , K 3 = yn + K 2 ? ? h 2 yn + K 2 ? 2 ? ? 2 ( xn + h ) . K 4 = yn + hK 3 ? ? yn + hK 3 ?右面列出计算结果 yn ，表中 y ( xn ) 仍表示准确解． 比较例 3 和例 2 的计算结果，显然以龙格-库 塔方法的精度为高．要注意，虽然四阶龙 格-库塔方法的计算量（每一步要 4 次计算 函数 f ）比改进的欧拉方法（它是一种二 阶龙格-库塔方法，每一步只要 2 次计算函 数 f ）打一倍，但由于这里放大了步长， 表和表所耗费的计算量几乎相同．这个例 子又一次显示了选择算法的重要意义． 然后值得指出的是，龙格-库塔方法的推导基于泰勒展开方法，因而它要求所 求的解具有较好的光滑性质．反之，如果解得光滑性差，那么，使用四阶龙格库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法．实际计算时，我 们应当针对问题的具体特点选择合适的算法．9.3.4表9 ?3计算结果 计算结果xn0.2 0.4 0.6 0.8 1.0yn1.7 1.5 1.7321y ( xn )1.7 1.5 1.7321变步长的龙格-库塔方法 变步长的龙格 库塔方法单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加了．步数的增加不但引起计算量的增大，而且可 能导致舍入误差的严重积累．因此同积分的数值计算一样，微分发出的数值解法 也有个选择步长的问题． 在选择步长时，需要考虑两个问题： 1o 怎样衡量和检验计算结果的精度？ 2o 如何依据所获得的精度处理步长？ 我们考察经典的四阶龙格-库塔公式 ( 3,13) ．从节点 xn 出发，先以 h 为步长求( 出一个近似值，记为 yn +)1 ，由于公式的局部截断误差为 O ( h5 ) ，故有h( y ( xn +1 ) ? yn +)1 ≈ ch 5 ，h( 3.14 )?h?? ? h ?2 然后将步长折半，即取 为步长从 xn 跨两步到 xn +1 ，再求得一个近似值 yn +1? ，每 2?h? 跨一步的截断误差是 c ? ? ，因此有 ?2?? ? ?h? ?2 y ( xn +1 ) ? yn +1? ≈ 2c ? ? ， ?2? ?h? 55( 3.15)1 ，即有 16比较 ( 3.14 ) 式和 ( 3.15 ) 式我们看到，步长折半后，误差大约减少到?2 y ( xn +1 ) ? yn +1? ?h? ? ?y ( xn +1 ) ? yn +1(h)≈1 . 16由此易得下列事后估计式y ( xn +1 ) ? y?h? ? ? ?2? n +1 ?h? ? ? 1 ? ?2? (h ? ≈ ? yn +1 ? yn +)1 ? . 15 ? ? ? ?这样，我们可以通过检查步长，折半前后两次计算结果的偏差( ?2 ? = yn +1? ? yn +)1h?h? ? ?来判定所选的步长是否合适，具体地说，将区分以下两种情况处理： ⒈ 对于给定的精度 ε ， 如果 ? & ε ， 我们反复将步长折半进行计算， 直至 ? & ε为止，这样取最终得到的 y?h? ? ? ?2? n +1作为结果；⒉ 如果 ? & ε ，我们将反复将步长加倍，直到 ? & ε 位止，这时再将步长折 半一次，就得到所要的结果． 这种通过加倍或折半步长的方法称为变步长方法 表面上看， 变步长方法． 为了选择步长， 变步长方法 每一步的计算量增加了，但总体考虑往往是合算的．9.4 单步法的收敛性与稳定性9.4.1 收敛性与相容性数值解法的基本思想是，通过某种离散化手段将微分方程 (1.1) 转化为差分方 程，如单步法 ( 2.10 ) ，即yn +1 = yn + h? ( xn , yn , h ) .( 4.1)xn ? x0 →0 n它 在 xn 处 的 解 为 yn ， 而 初 值 问 题 (1.1) , (1.2 ) 在 xn 处 的 精 确 解 为 y ( xn ) ， 记en = y ( xn ) ? yn 称为整体截断误差． 收敛性就是讨论当 x = xn 固定且 h =时 en → 0 的问题．定义 3 若一种数值方法（如单步法 ( 4.1) ）对于固定的 xn = x0 + nh ，当 h → 0 时有 yn → y ( xn ) ，其中 y ( x ) 是 (1.1) , (1.2 ) 的准确解，则称该方法是收敛 收敛的． 收敛 显然数值方法收敛是指 en = y ( xn ) ? yn → 0 ，对单步法 ( 4.1) 有下述收敛定理： 定理 1 茨条件 假设单步法 ( 4.1) 具有 p 阶精度，且增量函数 ? ( x, y, h ) 关于 y 满足利普希? ( x, y, h ) ? ? ? x, y, h ? ≤ L? y ? y ， ? ?? ? 又设是准确的，即 y0 = y ( x0 ) ，则其整体截断误差 整体截断误差y ( xn ) ? yn = O ( h p ) .???( 4.2 )( 4.3)证明设以 y n +1 表示取 yn = y ( xn ) 公式 ( 4.1) 求得的结果，即y n +1 = y ( xn ) + h? ( xn , y ( xn ) , h ) ，??( 4.4 )则 y ( xn +1 ) ? y n+1 为局部截断误差，由于所给方法具有 p 阶精度，按定义 2 ，存在定 数 C ，使y ( xn +1 ) ? y n +1 ≤ Ch p +1 ．?又由式 ( 4.4 ) 与 ( 4.1) ，得y n +1 ? yn +1 ≤ y ( xn ) ? yn+ h ? ( xn , y ( xn ) , h ) ? ? ( xn , yn , h ) .?利用假设条件 ( 4.2 ) ，有y n +1 ? yn +1 ≤ (1 + hL? ) y ( xn ) ? yn ,?从而有y ( xn +1 ) ? yn +1 ≤ y n +1 ? yn +1 + y ( xn +1 ) ? y n +1≤ (1 + hL? ) y ( xn ) ? yn + Ch p +1.? ?即对整体截断误差 en = y ( xn ) ? yn 成立下列递推关系式en +1 ≤ (1 + hL? ) en + Ch p +1 ，( 4.5)据此不等式反复递推，可得en ≤ (1 + hL? ) e0 +n n Ch p ? ? ?(1 + hL? ) ? 1? . ? L? ?( 4.6 )再注意到当 xn ? x0 = nh ≤ T 时(1 + hL ) ≤ ( e )n?hL?n≤eTL?，最终得下列估计式 en ≤ e0 eTL?+Ch p TL? e ?1 . L?()( 4.7 )由此可以断定，如果初值是准确的，即 e0 = 0 ，则 ( 4.3) 式成立，定理证毕． 依据这一定理，判断单步法 ( 4.1) 的收敛性，归结为验证增量函数能否满足利 普希茨条件． 对于欧拉方法， 由于其增量函数 ? 就是 f ( x, y ) ， 故当 f ( x, y ) 关于 y 满足利普 希茨条件时它是收敛的． 再考察改进的欧拉方法，其增量函数已由 ( 3.2 ) 式给出，这时有? ( x, y , h ) ? ? ? x, y , h ? ≤ [ f ( x, y ) ? f ( x, y ) ? ?? ? 1 2+ f (( x + h, y + hf ( x, y )) ? f ( x + h, y + hf ( x, y )) ] ．? ???假设 f ( x, y ) 关于 y 满足利普希茨条件，记利普希茨常数 L ，则由上式推得? ( x, y , h ) ? ? ? x , y , h ? ≤ L ? 1 + L ? y ? y . ? ? 2? ? 设限定 h ≤ h0 （ h0 为定数）上式表明 ? 关于 y 的利普希茨常数? h ? L? = L ?1 + 0 L ? , 2 ? ??? ?h ? ??因此改进的欧拉方法也是收敛的． 类似地，不难验证其他龙格-库塔方法的收敛性． 定理 1 表明时单步法收敛，并且当是初值问题的解，具有阶精度时，则有展开 式Tn +1 = y ( x + h) ? y ( x ) ? h? ( x, y ( x), h) = y ' ( x)h + y '' ( x ) 2 h +L 2' ? h[? ( x, y ( x ), 0) + ? x ( x, y ( x ), 0)h +L]= h[ y ' ( x ) ? ? ( x, y ( x ), h)] + O ( h 2 ).所以 p ≥ 1 的充要条件是 y ' ( x) ? ? ( x, y ( x), 0) = 0 ，而 y ' ( x) = f ( x, y ( x)) ，于是可给出如下定义： 定义 4 若单步法 ( 4.1) 的增量函数 ? 满足? ( x, y , 0) = f ( x, y ) ，则称单步法 ( 4.1) 与初值问题 (1.1) , (1.2 ) 相容． 以上讨论表明 p 阶方法 ( 4.1) 当 p ≥ 1 时与 (1.1) , (1.2 ) 相容，反之相容方法至少 是 1 阶的． 于是由定理 1 可知方法 ( 4.1) 收敛的充分必要条件是此方法是相容的．9.4.2绝对稳定性与绝对稳定前面关于收敛性的讨论有个前提， 必须假定数值方法本身的计算是准确的． 实 际情形并不是这样，差分方程的求解还会有计算误差，譬如由于数字舍入而引起 的小扰动．这类小扰动在传播过程中会不会恶性增长，以至于“淹没”了差分方 程的“真解”呢？这就是差分方法的稳定性问题．在实际计算时，我们希望某一 步产生的扰动值，在后面的计算中能够被控制，甚至是逐步衰减的． 定义 5 若一种数值方法在节点值 yn 上大小为 δ 的扰动，于以后各节点值ym ( m & n ) 上产生的误差均不超过 δ ，则称该方法是稳定 稳定的． 稳定下面先以欧拉法为例考察计算稳定性． 例4 考察初值问题? y ' = ?100 y, ? ? ? y ( 0 ) = 1. ?其准确解 y ( x ) = e?100x 是一个按指数曲线衰减得很快的函数，用欧拉法解方程y ' = ?100 y 得yn+1 = (1 ? 100h ) yn ．若取 h = 0.025 ，则欧拉公式的具体形式为yn +1 = ?1.5 yn ,计算结果列于表 9 ? 4 的第 2 列．我们看到，欧拉方法的解 yn 在准确值 y ( xn ) 的上下波动，计算过程明显地不稳定．但若取 h = 0.005, yn +1 = 0.5 yn 则计算过程稳定． 在考察后退的欧拉方法，取 h = 0.025 时计算公式为yn +1 = 1 yn ． 3.5计算结果列于表 9 ? 4 的第 3 列，这时计算过程是稳定的． 表9? 4 节 点 计算结果对比 后退欧拉方法欧拉方法0.025 0.050 0.075 0.100?1.5 2.25 ?3.375 5.06250.6 0.7例题表明稳定性不但与方法有关，也与步长 h 的大小有关，当然也与方程中 的 f ( x, y ) 有关．为了只考察数值方法本身．通常只检验将数值方法用于解模型方 程的稳定性，模型方程为y' = λ y ，( 4.8)其中 λ 为复数，这个方程分析较简单．对一般方程可以通过局部线性化化方程为? ?这种形式，例如在 ( x, y ) 的邻域，可展开为y ' = f ( x, y )= f ( x, y ) + f x' ( x, y )( x ? x) + f y' ( x, y )( y ? y ) +L,略去高阶项，再做变换即可得到 u ' = λ u 的形式．对于 m 个方程的方程组，可线性? ? ? ? ? ? ? ?? ?f 化为 y ' = Ay ，这里 A 为 m × m 的雅可比矩阵 ? i ? ?yi? ? ．若 A 有 m 个特征值 λ1 ,L , λm ， ?其中 λi 可能是复数．所以，为了使模型方程结果能推广到方程组，方程 ( 4.8) 中 λ 为复数．为保证微分方程本身的稳定性，还应假定 Re ( λ ) & 0 ． 下面先研究欧拉方法的稳定性．模型方程 y ' = λ y 的欧拉公式为yn+1 = (1 + hλ ) yn .( 4.9 )设在节点值 yn 上有一扰动值 ε n ，它的传播使节点值 yn +1 产生大小为 ε n +1 的扰动值，* * 假设用 yn = yn + ε n 按欧拉公式得出 yn +1 = yn +1 + ε n +1 的计算过程不再有新的误差，则扰动值满足ε n+1 = (1 + hλ ) ε n .可见扰动值满足原来的差分方程 ( 4.9 ) ．这样，如果差分方程的解是不增长的，即 有yn +1 ≤ yn ，则它就是稳定的．这一论断对于下面将要研究的其他方法同样适用． 显然，为要保证差分方程 ( 4.9 ) 的解是不增长的，只要选取 h 充分小，使1 + hλ ≤ 1.( 4.10 )在 ? = hλ 的复平面上，这是以 ( ?1, 0 ) 为圆心， 1 为半径的单位圆域．称为欧 拉法的绝对稳定域，一般情形可由下面定义． 定义 6 单步法 ( 4.1) 用于解模型方程 ( 4.8) ，若得到的解 yn+1 = E ( hλ ) yn ，满足 E ( hλ ) & 1 ， 则称方法 ( 4.1) 是绝对稳定 绝对稳定的． ? = hλ 的平面上， E ( hλ ) & 1 的 在 使 绝对稳定 变量围成的区域，称为绝对稳定域 绝对稳定域，它与实轴的交称为绝对稳定区间． 绝对稳定区间． 绝对稳定域 绝对稳定区间 对欧拉法 E ( hλ ) = 1 + hλ ，其绝对稳定域已由 ( 4.10 ) 给出，绝对稳定区间为?2 & hλ & 0 ，在例 5 中 λ = ?100, ?2 & ?100h & 0 ，即 0 & h & 2 /100 = 0.02 为绝对稳定区间，例 4 中取 h = 0.025 故它是不稳定的，当取 h = 0.005 时它是稳定的． 对二阶 R ? K 方法，解模型方程 ( 4.8) 可得到2 ? ( hλ ) ? y , yn +1 = ?1 + hλ + ? n 2 ? ? ? ?故( hλ ) E ( hλ ) = 1 + hλ +22.绝对稳定域由( hλ ) 1 + hλ +22&1 得 到 ， 于 是 可 得 到 绝 对 稳 定 区 间 为?2 & hλ & 0 ，即 0 & h & ?2 / λ ．类似可得到三阶及四阶的 R ? K 方法的 E ( hλ ) 分别为( hλ ) E ( hλ ) = 1 + hλ +2!2( hλ ) +3!33，4( hλ ) E ( hλ ) = 1 + hλ +2!2( hλ ) +3!( hλ ) +4!.由 E ( hλ ) & 1 可得到相应的绝对稳定域．当 λ 为实数时则得绝对稳定区间．它们 分别为 三阶 R ? K 显示方法： ?2.51 & hλ & 0, 即 0 & h & ?2.51/ λ. 四阶 R ? K 显示方法： ?2.78 & hλ & 0, 即 0 & h & ?2.78 / λ ． 从以上讨论可知显示的 R ? K 方法的绝对稳定域均为有限域，都对步长 h 有限 制．如果 h 不在所给的绝对稳定区间内，方法就不稳定． 例5 y ' = ?20 y( 0 ≤ x ≤ 1) ， y ( 0 ) = 1 分 别 取 h = 0.1 及 h = 0.2 用 经 典 的R ? K 四阶方法 ( 3.13) 计算．解本例 λ = ?20 ， λ h 分别为 ?2 及 ?4 ，前者在绝对稳定区间内，后者则不在，用四阶 R ? K 方法计算其误差见下表：xn0.20.93 × 10?10.40.12 × 10?10.60.14 × 10?20.80.15 × 10?31.00.17 × 10 ?4h = 0.1 h = 0.24.9825.0125.0625.03125.0从以上结果看到，如果步长 h 不满足绝对稳定条件，误差增长很快． 对隐式单步法，可以同样讨论方法的稳定性，例如对后退欧拉法，用它解模 型方程可得故 由 E ( hλ ) =1 &1 1 ? hλ1 yn , 1 ? hλ 1 ． E ( hλ ) = 1 ? hλ yn +1 =可得绝对稳定域为 1 ? hλ & 1 ， 它是以 (1, 0 ) 为圆心，故绝对稳定区间为 ?∞ & λ h & 0 ． λ & 0 时， 0 & h & ∞ ， 当 则 1 为半径的单位圆外部， 即对任何步长均为稳定的． 对隐式梯形法，它用于解模型方程 ( 4.8) 得 hλ 2 y, yn +1 = hλ n 1? 2 1 + hλ / 2 故 E ( hλ ) = . 1 ? hλ / 2 hλ 1+ 2 & 1 ，故绝对稳定域为 ? = hλ 的左半平面，绝 对 Re ( λ ) & 0 有 E ( hλ ) = hλ 1? 21+对稳定区间为 ?∞ & λ h & 0 ，即 0 & h & ∞ 时梯形法均是稳定的．9.5线性多步法在逐步推进的求解过程中，计算 yn +1 之前事实上已经求出了一系列的近似值 y0 , y1 ,L , yn ，如果充分利用前面多步的信息来预测，则可以期望会获得较高的精 度．这就是构造所谓线性多步法的基本思想． 构造多步法的主要途径是基于数值积分方法和基于泰勒展开方法，前者可直 接由方程 (1.1) 两端积分后利用插值求积公式得到．本节主要介绍基于泰勒展开的 构造方法．9.5.1线性多步法的一般公式如果计算 yn + k 时，除用 yn + k ?1 的值，还用到 yn+i ( i = 0,1,L , k ? 2 ) 的值，则称此 方法为线性多步法 线性多步法．一般的线性多步法公式可表示为 线性多步法yn + k = ∑ α i yn + i + h ∑ β i f n + i ,i=0 i=0k ?1k( 5.1)其中 yn +i 为 y ( xn+i ) 的近似，f n +i = f ( xn +i , yn +i ) , xn +i = x0 + ih, α i , βi 为常数， 0 及 β 0 不 α 全为零，则称 ( 5.1) 为线性 k 步法，计算时需先给出前面 k 个近似值 y0 , y1 ,L , yk ?1 ， 再由 ( 5.1) 逐次求出 yk , yk +1 ,L ．如果 β k = 0 ，称 ( 5.1) 为显示 k 步法，这时 yn + k 可直 接由 ( 5.1) 算出;如果 β k ≠ 0 ，则 ( 5.1) 称为隐式 k 步法，求解时与梯形法 ( 2.7 ) 相同， 要用迭代方法可算出 yn + k ．( 5.1) 中系数 αi 及 βi 可根据方法的局部截断误差及阶确 定，其定义为： 定义 7 设 y ( x ) 是初值问题 (1.1) , (1.2 ) 的准确解，线性多步法 ( 5.1) 在 xn + k 上的局部截断误差为Tn + k = L ? y ( xn ) ; h ? ? ?= y ( xn + k ) ? ∑ α i y ( xn +i ) ? h∑ β i y ' ( xn +i ) .i =0 i =0k ?1k( 5.2 )若 Tn + k = O ( h p +1 ) ， 则称方法 ( 5.1) 是 p 阶的，p ≥ 1 则称方法 ( 5.1) 与方程 (1.1) 是 相容的． 由定义 7 ，对 Tn + k 在 xn 处做泰勒展开，由于 y ( xn + ih ) = y ( xn ) + ihy ( xn )'( ih ) +2!2y '' ( xn )( ih ) +3!' '3y ''' ( xn ) + L ，y ( xn + ih ) = y ( xn ) + ihy ( xn )''( ih ) +2!2y ''' ( xn ) + L代入 ( 5.2 ) 得Tn + k = c0 y ( xn ) + c1hy ' ( xn ) + c2 h2 y '' ( xn )+ L + c p h p y ( p ) ( xn ) + L ，( 5.3)其中c0 = 1 ? (α 0 + L + α k ?1 ) ,c1 = k ? ?α1 + 2α 2 + L + ( k ? 1) α k ?1 ? ? ( β 0 + β1 + L + β k ) , ? ?cp =1? q q k ? α1 + 2q α 2 + L + ( k ? 1) α k ?1 ? ? ? q!()?1 ? β + 2q ?1 β 2 + L + k q ?1β k ? ? ( q ? 1)! ? 1q = 2,3,L.( 5.4 )若在公式 ( 5.1) 中选择系数 αi 及 βi ，使它满足c0 = c1 = L = c p = 0, c p +1 ≠ 0由定义可知此时所构造的多步法是 p 阶的，且Tn + k = c p +1h p +1 y ( p +1) ( xn ) + O ( h p + 2 ) .( 5.5)称右端第一项为局部截断误差主项 c p +1 称为误差常数． 局部截断误差主项， 局部截断误差主项 根据相容性定义， p ≥ 1 ，即 c0 = c1 = 0 ，由 yn +1 = yn + hf n , ( 5.4 ) 得 ?α 0 + α1 + L + α k ?1 = 1, ? k ?1 k ? ∑ iαi + ∑ βi = k. ? i =0 ? i =1 故方法与微分方程相容的充分必要条件是成立． 显然，当 k = 1 时，若 β1 = 0 ，则由 ( 5.6 ) 可求得( 5.6 )α 0 = 1, β0 = 1.此时公式 ( 5.1) 为 yn +1 = yn + hf n , 即为欧拉法．从 ( 5.4 ) 可求得 c2 = 1/ 2 ≠ 0 ，故方法为 1 阶精度，且局部截断误差为Tn +1 = 1 2 '' h y ( xn ) + O ( h 3 ) ， 2这和第 2 节给出的定义及结果是一致的．对 k = 1 ， 若 β1 ≠ 0 ， 此 时 方 法 为 隐 式 公 式 ， 为 了 确 定 系 数 α 0 , β 0 , β1 ， 可 由 c0 = c1 = c2 = 0 解得 α 0 = 1, β 0 = β1 = 1/ 2 ．于是得到公式yn +1 = yn + h ( f n + f n+1 ) ， 2即为梯形法．由 ( 5.4 ) 可求得 c3 = ?1/12 ，故 p = 2 ，所以梯形法是二阶方法，其局 部截断误差主项是 ?h3 y ''' ( xn ) /12 ，这与第 2 节中的讨论也是一致的． 对 k ≥ 2 的多步法公式都可利用 ( 5.4 ) 确定系数 α i , β i ，并由 ( 5.5 ) 给出局部截断 误差，下面只就若干常用的多步法导出具体公式．9.5.2阿当姆斯显示与隐式公式考虑形如yn + k = yn + k ?1 + h∑ β i f n + ii =0 k( 5.7 )的 k 步法，称为阿当姆斯 Adams）方法 β k = 0 为显示方法，β k ≠ 0 为隐式方法， 阿当姆斯（ 方法． 阿当姆斯 方法 通常称为阿当姆斯显示与隐式公式， 也称 Adams-Bashforth 公式与 Adams-Monlton 公式．这类公式可直接由方程 (1.1) 两端积分（从 xn + k ?1 到 xn + k 积分）求得．下面可 利用 ( 5.4 ) 由 c1 = L = c p = 0 推出，对比 ( 5.7 ) 与 ( 5.1) 可知此时系数α 0 = α1 = L = α k ? 2 = 0 ，显然 c0 = 0 成立，下面只需确定系数 β 0 , β1 ,L , β k ，故可令c1 = L = ck +1 = 0 ，则可求得 β 0 , β1 ,L , β k （若 β k = 0 ，则令 c0 = L = ck = 0 来求得β 0 , β1 ,L , β k ?1 ） ．下面以 k = 3 为例，由 c1 = c2 = c3 = c4 = 0 ，根据 ( 5.4 ) 得? β 0 + β1 + β 2 + β3 = 1, ? 2( β + 2β + 3β ) = 5, ? 1 2 3 ? ? 3( β1 + 4β 2 + 9β3 ) = 19, ?4( β1 + 8β 2 + 27 β3 ) = 65. ?若 β 3 = 0 ，则由前三个方程解得β0 =5 , 12β1 = ?16 , 12β2 =5 12得到 k = 3 的阿当姆斯显示公式是yn + 3 = yn + 2 +h (23 f n + 2 ? 16 f n +1 + 5 f n ). 12( 5.8)由 ( 5.4 ) 求得 c4 = 3 / 8 ，所以 ( 5.8 ) 是三阶方法，局部截断误差是3 Tn +3 = h 4 y ( 4) ( xn ) + O( h5 ). 8若 β 3 ≠ 0 ，则可解得β0 =1 , 24β1 = ?5 , 24β2 =19 , 24β3 = ．3 8于是得 k = 3 的阿当姆斯隐式公式为 yn + 3 = yn + 2 + h (9 f n +3 + 19 f n + 2 ? 5 f n +1 + f n ), 2419 5 ( 5) h y ( xn ) + O( h 6 ). 720( 5.9 ) ( 5.10 )它是四阶方法，局部截断误差是 Tn +3 = ?用类似的方法可求得阿当姆斯显示方法和隐式方法的公式， 9 ? 5 及表 9 ? 6 表 分别列出了 k = 1, 2, 3, 4 时的阿当姆斯显示公式与阿当姆斯隐式公式，其中 k 为步 数， p 为方法的阶， c p +1 为误差常数． 表9?5 阿当姆斯显示公式kp公式c p +11 21 2yn +1 = yn + hf n h yn + 2 = yn +1 + (3 f n +1 ? f n ) 2 h yn +3 = yn + 2 + (23 f n + 2 ? 16 f n +1 + 5 f n ) 12 h yn +3 = yn + 2 + (9 f n +3 + 19 f n + 2 ? 5 f n +1 + f n ) 2434341 2 5 12 3 8 251 720表9? 6阿当姆斯隐式公式kp公式c p +11 21 23434h yn +1 = yn + ( f n +1 + f n ) 2 h yn +3 = yn + 2 + (5 f n + 2 + 8 f n +1 ? f n ) 12 h yn +3 = yn + 2 + (9 f n +3 + 19 f n + 2 ? 5 f n +1 + f n ) 24 h yn + 4 = yn + 3 + (251 f n + 4 + 646 f n +3 ? 264 f n + 2 7201 2 5 12 3 8 251 720+106 f n +1 ? 19 f n )例6用四阶阿当姆斯显示和隐式方法解初值问题 y ' = ? y + x + 1,y (0) = 1.取步长 h = 0.1 ． 解 本题 f n = ? yn + xn + 1, xn = nh = 0.1n. 从四阶阿当姆斯显示公式得到 yn + 4 = yn + 3 + = h (55 f n +3 ? 59 f n + 2 + 37 f n +1 ? 9 f n ) 241 (18.5 yn +3 + 5.9 yn + 2 ? 3.7 yn +1 + 0.9 yn 24+0.24n + 3.24).对于四阶阿当姆斯隐式公式得到 yn + 3 = yn + 2 + = h (9 f n +3 + 19 f n + 2 ? 5 f n +1 + f n ) 241 ( ?0.9 yn +3 + 22.1 yn + 2 + 0.5 yn +1 ? 0.1 yn + 0.24n + 3). 24由此可以直接解出 yn +3 而不用迭代，得到 yn + 3 =1 (22.1 yn + 2 + 0.5 yn +1 ? 0.1 yn + 0.24n + 3). 24.9计算结果见表 9 ? 7 ，其中显示方法中的 y0 , y1 , y2 , y3 及隐式方法中的 y0 , y1 , y2 均用 准确解 y ( x) = e ? x + x 计算得到，对一般方程，可用四阶 R ? K 方法计算初始近似．表 9 ? 7 计算结果 精确解 y ( xn )xn= e ? xn + xn阿当姆斯显示方法y n阿当姆斯隐式方法y n y ( xn ) ? yny ( xn ) ? yn上例子看到同阶的阿当姆斯方法， 隐式方法要比显示方法误差小， 这可从以上 例子看到同阶的阿当姆斯方法，隐式方法要比显示方法误差小，这可以从两种方 法的局部截断误差主项 c p +1h p +1 y ( p +1) ( xn ) 的系数大小得到解释，这里 c p +1 分别为251/ 720 及 ?19 / 720 ． 9.5.3米尔尼方法与辛普森方法考虑与 ( 5.7 ) 不同的另一个 k = 4 的显示公式yn +1 = yn + h ( β3 f n +3 + β 2 f n + 2 + β1 f n +1 + β 0 f n ) ,其中 β 0 , β1 , β 2 , β 3 为待定常数，可根据使公式的阶尽可能高这一条件来确定其数 值．由 ( 5.4 ) 可知 c0 = 0 ，再令 c1 = c2 = c3 = c4 = 0 得到? β 0 + β1 + β 2 + β3 = 4, ? 2 β + 2β + 3β = 16, ? ( 1 2 3) ? ? 3 ( β1 + 4β 2 + 9β3 ) = 64, ?4 ( β1 + 8β 2 + 27 β3 ) = 256. ?解此方程组得β 3 = , β 2 = ? , β1 = , β 0 = 0于是得到四步显示公式yn + 4 = yn + 4h ( 2 f n +3 ? f n + 2 + 2 f n +1 ) , 38 34 38 3( 5.11)称为米尔尼 Milne）方法 米尔尼（ 方法.由于 c5 = 14 / 45 ，故方法为 4 阶，其局部截断误差为 米尔尼 方法Tn + 4 = 14 5 ( 5) h y ( xn ) + O ( h 6 ) . 45( 5.12 )米尔尼方法也可以通过方程 (1.1) 两端积分y ( xn + 4 ) ? y ( xn ) = ∫xn+ 4 xnf ( x, y ( x ) )dx得到．若将方程 (1.1) 从 xn 到 xn + 2 积分，可得y ( xn + 2 ) ? y ( xn ) = ∫xn +2 xnf ( x, y ( x ) )dx.右端积分通过辛普森求积公式就有yn + 2 = yn + h ( f n + 4 f n+1 + f n + 2 ) . 3( 5.13)称为辛普森方法 辛普森方法．它是隐式二步四阶方法，其局部截断误差为 辛普森方法 Tn + 2 = ? 汉明方法 辛普森公式是二步方法中阶数最高的， 但它的稳定性较差， 为了改善稳定性， 我们考察另一类三步法公式 h5 ( 5) y ( xn ) + O ( h 6 ) . 90( 5.14 )yn +3 = α 0 yn + α1 yn+1 + α 2 yn + 2 + h ( β1 f n +1 + β 2 f n + 2 + β3 f n +3 ) ,其中系数 α 0 , α1 , α 2 及 β1 , β 2 , β3 为常数，如果希望导出的公式是四阶的，则系数中 至少有一个自由参数．若取 α1 = 1 ，则可得到辛普森公式．若取 α1 = 0 ，仍利用泰 勒展开，由 ( 5.4 ) ，令 c1 = c2 = c3 = c4 = 0 ，则可得到α 0 + α 2 = 1, ? ? 2α 2 + β1 + β 2 + β3 = 3, ? ? ? 4α 2 + 2 ( β1 + 2β 2 + 3β3 ) = 9, ? 8α + 3 ( β + 4β + 9β ) = 27, 2 1 2 3 ? ?16α 2 + 4 ( β1 + 8β 2 + 27 β3 ) = 81. ?解此方程组得1 8 9 8 3 8 6 8 3 8α0 = ? ,于是有yn + 3 =α2 = ,β1 = ? ,β2 = ,β3 = .1 3h ( 9 yn+ 2 ? yn ) + ( f n +3 + 2 f n +2 ? f n +1 ) , 8 8( 5.15)称为汉明(Hamming)方法 汉明 方法．由于 c3 = ?1 / 40 ，故方法是四阶的，而局部截断误差为 Tn +3 = ?9.5.5h 5 ( 5) y ( xn ) + O ( h 6 ) . 40( 5.16 )预测-校正方法 预测 校正方法对于隐式的线性多步法，计算时要进行迭代，计算量较大．为了避免进行迭( 代， 通常采用显示公式给出 yn + k 的一个初始近似， 记为 yn0 )k ， 称为预测 Predictor） （ ， 预测 +接着计算 f n + k 的值（evaluation） ，再用隐式公式计算 yn + k ，称为校正(corrector)．例 校正 如在 ( 2.13) 中用欧拉法做预测，再用梯形法校正，得到改进欧拉法，它就是一个 二阶预测-校正方法．一般情况下，预测公式与校正公式都取同阶的显示方法与隐 式方法相匹配．例如用四阶的阿当姆斯显示方法做预测，再用四阶阿当姆隐式公 式做校正，得到以下格式： 预测 P： ynp+ 4 = yn +3 +h ( 55 f n +3 ? 59 f n + 2 + 37 f n +1 ? 9 f n ) , 24求值 E： f.np+ 4 = f ( xn + 4 , ynp+ 4 ) , 校正 C： yn + 4 = yn +3 +h ( 9 fnp+4 + 19 f n+3 ? 5 fn+2 + fn+1 ) , 24求值 E： f n + 4 = f ( xn + 4 , yn + 4 ) . 此公式称为阿当姆斯四阶预测 校正格式（PECE） 阿当姆斯四阶预测-校正格式 ． 阿当姆斯四阶预测 校正格式 依据四阶阿当姆斯公式的截断误差，对于 PECE 的预测步 P 有y ( xn + 4 ) ? ynp+ 4 ≈251 5 ( 5) h y ( xn ) , 720对校正步 C 有y ( xn + 4 ) ? ynp+ 4 ≈ ?19 5 ( 5) h y ( xn ) . 720两式相减得h5 y (5)( xn ) ≈ ?720 p ( yn + 4 ? yn + 4 ) , 270于是有下列事件后误差估计y ( xn + 4 ) ? ynp+ 4 ≈ ? y ( xn + 4 ) ? ynp+ 4251 p ( yn + 4 ? yn + 4 ) , 720 19 ≈ ( ynp+4 ? yn+4 ) . 270容易看出 ynpm4 = ynp+ 4 + +_y n+4251 ( yn+4 ? ynp+4 ) ,? ? ? 720 ? 19 p = yn + 4 ? ( yn + 4 ? y n + 4 ) . ? ? 270 ?( 5.17 )比 ynp+ 4 , yn + 4 更好．但在 ynpm4 的表达式中 yn + 4 是未知的，因此计算时用上一步代替， + 从而构造一种修正预测 校正格式（PMECME） 修正预测-校正格式 ： 修正预测 校正格式（P： ynp+ 4 = yn +3 + M： ynpm4 + h ( 55 f n +3 ? 59 f n + 2 + 37 f n +1 ? 9 f n ) , 24 251 c = ynp+ 4 + ( yn+3 ? ynp+3 ) , 720E： f npm = f ( xn + 4 , ynpm4 ) , +4 +c C： y n + 4 = y n + 3 +M： y n + 4h ( 9 fnpm4 + 19 fn+3 ? 5 f n+2 + f n+1 ) , + 24 19 c c = yn + 4 ? ( yn+4 ? ynp+4 ) , 270E： f n + 4 = f ( xn + 4 , yn + 4 ) .c 注意：在 PMECME 格式中已将 ( 5.17 ) 的 yn + 4 及 y n + 4 分别改为 yn + 4 及 yn + 4 ． _利 用 米 尔 尼 公 式 ( 5.11) 和 汉 明 公 式 ( 5.15 ) 相 匹 配 ， 并 利 用 截 断 误 差 修正米尔尼-汉明预测 校正格式 修正米尔尼 汉明预测-校正格 ( 5.12 ) , ( 5.16 ) 改进计算结果，可类似地建立四阶修正米尔尼 汉明预测 校正格式 ： （PMECME）4 P： ynp+ 4 = yn + h ( 2 f n +3 ? f n + 2 + 2 f n +1 ) , 3 112 c M： ynpm4 = ynp+ 4 + ( yn+3 ? ynp+3 ) , + 121 E： f npm = f ( xn + 4 , ynpm4 ) , +4 +c C： y n + 4 =M： y n + 41 3 ( 9 yn +3 ? yn+1 ) + h ( f npm + 2 f n +3 ? f n+ 2 ) , +4 8 8 9 c = yn + 4 ? ( ync+4 ? ynp+4 ) , 121E： f n + 4 = f ( xn + 4 , yn + 4 ) .例7将例 6 的初值问题用修正的米尔尼-汉明预测-校正公式计算 y5 及 y6 ，初 值 y0 , y1 , y2 , y3 仍 用 已 算 出 的 精 确 解 ， 即 y0 = 1, y1 = 1., y2 = 1., y3 = 1. ，给出计算结果及误差． 解 根据修正的米尔尼-汉明预测-校正公式可得y5p = 1., y5pm = y5p + 112 c × ( y4 ? y4p ) = 1.. 121c （注： y4p = 1., y4 = 1. ）f ( x5 , y5pm ) = ? y5pm + x5 + 1 = 0., 1 0.3 c y5 = × ( 9 y4 ? y2 ) + × ? f ( x5 , y5pm ) + 2 f 4 ? f3 ? ? 8 8 ?= 1.,c y5 = y5 ?9 c × ( y5 ? y5p ) = 1.1f 5 = ? y5 + x5 + 1 = 0., y6p = y2 + y6pm0.4 × ( 2 f5 ? f 4 + 2 f 3 ) = 1. 112 c = y6p + × ( y5 ? y5p ) = 1.1f 6pm = ? y6pm + x6 + 1 = 0.,1 0.3 c y6 = × ( 9 y5 ? y3 ) + × ( f 6pm + 2 f 5 ? f 4 ) = 1. 8c y6 = y6 ?9 c × ( y6 ? y6p ) = 1.1y ( x5 ) ? y5 = 4.97 × 10 ?8 , y ( x6 ) ? y6 = 2.61× 10?8.误差从结果看，此方法的误差比四阶阿当姆斯隐式法和四阶汉明方法小，这与理论分 析一致．9.5.6构造多步伐公式的注记和例前面已指出构造多步法公式有基于数值积分和泰勒展开两种途径，只对能将 微分方程 (1.1) 转化为等价的积分方程的情形方可利用数值积分方法建立多步法 公式，它是有局限性的．即前种途径只对部分方法适用．而用泰勒展开则可构造 任意多步法公式，其做法是根据多步法公式的形式，直接在 xn 处做泰勒展开即 可．不必套用系数公式 ( 5.4 ) 确定多步法 ( 5.1) 的系数 α i 及 βi （ i = 0,1, L , k ） ，因为 多步法公式不一定如 ( 5.1) 的形式．另外，套用公式容易记错．具体做法见下面例 子． 例 8 解 初 值 问 题 y ' = f ( x, y ) , y ( x0 ) = y0 ． 用 显 示 二 步 法yn +1 = α 0 yn + α1 yn ?1 + h ( β 0 f n + β1 f n ?1 ) ，其中 f n = f ( xn , yn ) , f n ?1 = f ( xn ?1 , yn ?1 ) ．试确定参数 α 0 , α1 , β 0 , β1 使方法阶数尽可能高，并求局部截断误差． 用泰勒展开确定参数满足的方程， 由于 解 本题仍根据局部截断误差定义，Tn +1 = y ( xn + h ) ? α 0 y ( xn ) ? α1 y ( xn ? h ) ? h ? β 0 y ' ( xn ) + β1 y ' ( xn ? h ) ? ? ?= y ( xn ) + hy ' ( xn ) + +h 2 '' h3 y ( xn ) + y ''' ( xn ) 2 3!h 4 ( 4) h2 y ( xn ) + O ( h5 ) ? α 0 y ( xn ) ?α1[ y ( xn ) ? hy ' ( xn ) + y '' ( xn ) 4! 2h3 ''' h 4 ( 4) ? y ( xn ) + y ( xn ) + O ( h5 )] ? β 0 hy ' ( xn ) ? β1h[ y ' ( xn ) ? hy '' ( xn ) 3! 4!+h 2 ''' h3 4 y ( xn ) ? y ( ) ( xn ) + O ( h 4 )] 2 3!= (1 ? α 0 ? α1 ) y ( xn ) + (1 + α1 ? β 0 ? β1 ) hy ' ( xn )1 ? ?1 1 ? ?1 1 + ? ? α1 + β1 ? h 2 y '' ( xn ) + ? + α1 ? β1 ? h3 2 ? ?2 2 ? ?6 6 1 1 ? ? 1 . y ''' ( xn ) + ? ? α1 + β1 ? h 4 y ( 4 ) ( xn ) + O ( h5 ) 6 ? ? 24 24为求参数 α 0 , α1 , β 0 , β1 使方法阶数尽量高，可令1 ? α 0 ? α1 = 0 ，1 1 ? α1 + β1 = 0 ， 2 21 + α1 ? β 0 ? β1 = 0 ，1 1 1 + α1 ? β = 0 ． 6 6 2即得方程组? α 0 + α1 = 1, ? ??α1 + β 0 + β1 = 1, ? ? α1 ? 2 β1 = 1, ? ?α1 + 3β1 = 1. ?解得 α 0 = ?4, α1 = 5, β 0 = 4, β1 = 2 ，此时公式为三阶，而且1 4 Tn +1 = h 4 y ( ) ( xn ) + O ( h5 ) 6即为所求局部解得误差．而所得二步法为yn +1 = ?4 yn + 5 yn ?1 + 2h ( 2 f n + f n ?1 ) .例9 证明存在 α 的一个值，使线性多步法yn +1 + α ( yn ? yn ?1 ) ? yn ? 2 =1 ( 3 + α ) h ( f n + f n?1 ) 2是四阶的． 证明 由于Tn +1 = y ( xn + h ) + α ? y ( xn ) ? y ( xn ? h ) ? ? y ( xn ? 2h ) ? ? ? 1 ( 3 + α ) h ? y ' ( xn ) + y ' ( xn ? h ) ? ? ? 2只要证明局部截断误差 Tn +1 = O ( h5 ) ，则方法为四阶．仍用泰勒展开，= y ( xn ) + hy ' ( xn ) +h 2 '' h3 h 4 ( 4) y ( xn ) + y ''' ( xn ) + y ( xn ) 2 3! 4!? ? h2 h3 h4 +O ( h5 ) ? α ?( ? h ) y ' ( xn ) + y '' ( xn ) ? y ''' ( xn ) + y ( 4 ) ( xn ) + O ( h5 ) ? 2 3! 4! ? ?2 3 4 ? ( 2h ) y '' x ? ( 2h ) y ''' x + ( 2h ) y ( 4) x + O h5 ? ' ? ? y ( xn ) ? 2hy ( xn ) + ( n) ( n) ( n ) ( )? 2 3! 4! ? ? ? ?? ' ? h h ''' h3 ( 4 ) ' '' ? ( 3 + α ) ? y ( xn ) + y ( xn ) ? hy ( xn ) + y ( xn ) ? y ( xn ) + O ( h 4 ) ? 2 2 3! ? ?1 ?1 1 ? = ?1 + α + 2 ? ( 3 + α ) ? hy ' ( xn ) + ? ? α ? 2 + ( 3 + α ) ? h 2 y '' ( xn ) ? ? 2 ?2 2 ? 4 1 1 2 1 ?1 1 ? ?1 ? + ? + α + ? ( 3 + α ) ? h3 y ''' ( xn ) + ? ? α ? + ( 3 + α ) ? 3 4 3 12 ?6 6 ? ? 24 24 ?×h 4 y ( 4 ) ( xn ) + O ( h5 )1 ?4 1 ? 4 = ? ? α ? h3 y ''' ( xn ) + ( ?9 + α ) h 4 y ( ) ( xn ) + O ( h5 ) . 24 ? 3 12 ?当 α = 9 时， Tn +1 = O ( h5 ) ，故方法是四阶的．9.6 9.6.1方程组和高阶方程一阶方程组前面我们研究了单个方程 y ' = f 的数值解法，只要把 y 和 f 理解为向量，那 么，所提供的各种计算公式即可应用到一阶方程组的情形． 考察一阶方程组y 'i = f i ( x, y1 , y2 ,L, yN )的初值问题，初始条件给为( i = 1, 2,L, N )yi ( x0 ) = yi0若采用向量的记号，记( i = 1, 2,L, N )Ty = ( y1 , y2 , L , y N ) ,0 0 y0 = ( y10 , y2 , L , y N ) , Tf = ( f1 , f 2 , L , f N ) .T则上述方程组的初值问题可表示为 y ' = f ( x, y ) , ? ? ? y ( x0 ) = y0 . ? ? 求解这一初值问题的四阶龙格-库塔公式为yn +1 = yn + h ( k1 + 2k2 + 2k3 + k4 ) , 6( 6.1)式中k1 = f ( xn , yn ) ,h h ? ? k2 = f ? xn + , yn + k1 ? , 2 2 ? ? h h ? ? k3 = f ? xn + , yn + k2 ? , 2 2 ? ?k4 = f ( xn + h, yn + hk3 ) .或表示为yi ,n +1 = yin + h ( Ki1 + 2 Ki 2 + 2 Ki 3 + Ki 4 ) 6( i = 1, 2,L, N )其中K i1 = fi ( xn , y1n , y2 n ,L, y Nn ) ,h h h h ? ? K i 2 = f i ? xn + , y1n + K11 , y2 n + K 21 , L , y Nn + K N 1 ? , 2 2 2 2 ? ? h h h h ? ? K i 3 = fi ? xn + , y1n + K12 , y2 n + K 22 , L , y Nn + K N 2 ? , 2 2 2 2 ? ? h ? ? K i 4 = f i ? xn + , y1n + hK13 , y2 n + hK 23 , L , yNn + hK N 3 ? . 2 ? ?这里 yin 是第 i 个因变量 yi ( x ) 在节点 xn = x0 + nh 的近似值． 为了帮助理解这一公式的计算过程，我们再考察两个方程的特殊情形：? y ' = f ( x, y , z ) , ? ' ? z = g ( x, y , z ) , . ? ? y ( x0 ) = y0 , ?z x = z ? ( 0) 0这时四阶龙格-库塔公式具有形式 h ( K1 + 2 K 2 + 2 K3 + K 4 ) , ? ? ? 6 ? h zn +1 = zn + ( L1 + 2 L2 + 2 L3 + L4 ) . ? ? 6 ? yn +1 = yn + 其中( 6.2 )K1 = f ( xn , yn , zn ) ,? ? h h h ? ? ? K 2 = f ? xn + , yn + K1 , zn + L1 ? , 2 2 2 ? ? ? ? h h h ?? ? K 3 = f ? xn + , yn + K 2 , zn + L2 ? , ? 2 2 2 ? ? ? K 4 = f ( xn + h, yn + hK 3 , zn + hL3 ) , ? ? L1 = g ( xn , yn , zn ) , ? h h h ? ? ? L2 = g ? xn + , yn + K1 , zn + L1 ? , ? 2 2 2 ? ? ? ? h h h ? ? L3 = g ? xn + , yn + K 2 , zn + L2 ? , ? 2 2 2 ? ? ? ? L4 = g ( xn + h, yn + hK 3 , zn + hL3 ) . ?( 6.3)这 是 一 步 法 ， 利 用 节 点 xn 上 的 值 yn , zn ， 由 ( 6.3) 式 顺 序 计 算 K1 , L1 , K 2 , L2 , K 3 , L3 , K 4 , L4 ，然后代入 ( 6.2 ) 式即可求得节点 xn +1 上的 yn +1 , zn +1 ．9.6.2化高阶方程为一阶方程组关于高阶微分方程（或方程组）的初值问题，原则上总可以归结为一阶方程 组来求解．例如，考察下列阶微分方程y ( m ) = f ( x, y, y ' , L , y m ?1 ) ,( 6.4 ) ( 6.5)初始条件为y ( x0 ) = y0 , y ' ( x0 ) = y '0 , L , y (m ?1)( x0 ) = y0( m?1) .只要引进新的变量y1 = y, y2 = y ' ,L, ym = y ( m ?1) ,即可将 m 阶方程 ( 6.4 ) 化为如下的一阶方程组：? ? y 2 = y3 , ? ? LLL ? ? ' y m ?1 = ym , ? y 'm = f ( x, y1 , y2 , L , ym ) .? ? y '1 = y2 ,'( 6.6 )初始条件 ( 6.5) 则相应地化为 y1 ( x0 ) = y0 , y2 ( x0 ) = y '0 , L , ym ( x0 ) = y0 ( m ?1) . 不难证明初值问题 ( 6.4 ) , ( 6.5) 和 ( 6.6 ) , ( 6.7 ) 是彼此等价的． 特别的，对于下列二阶方程的初值问题：? y '' = f ( x, y, y ' ) , ? ? ? y ( x0 ) = y0 , ? ' ' ? y ( x0 ) = y 0 . ?( 6.7 )引进新的变量 z = y ' ，即可化为下列方程组的初值问题：? y ' = z, ? ' ? z = f ( x, y, z ) , ? ' ? z ( x0 ) = y 0 .针对这个问题应用四阶龙格―库塔公式 ( 6.2 ) ，有h ? ? yn +1 = yn + 6 ( K1 + 2 K 2 + 2 K 3 + K 4 ) , ? ? ? z = z + h ( L + 2L + 2L + L ) . n 1 2 3 4 ? n +1 6 ?由 ( 6.3) 式可得K1 = zn , K 2 = zn +L1 = f ( xn , yn , zn ) ;h h h h ? ? L1 , L2 = f ? xn + , yn + K1 , zn + L1 ? ; 2 2 2 2 ? ?K 3 = zn +h h h h ? ? L2 , L3 = f ? xn + , yn + K 2 , zn + L2 ? ; 2 2 2 2 ? ?K 4 = zn + hL3 , L4 = f ( xn + h, yn + hK 3 , zn + hL3 ) .如果消去 K1 , K 2 , K 3 , K 4 ，则上述格式可表示为 ? h2 ? yn +1 = yn + hzn + ( L1 + L2 + L3 ) , ? 6 ? ? z = z + h ( L + 2L + 2L + L ) . n 1 2 3 4 ? n +1 6 ? 这里L1 = f ( xn , yn , zn ) ,h h h ? ? L2 = f ? xn + , yn + zn , zn + L1 ? , 2 2 2 ? ?? h h h2 h ? L3 = f ? xn + , yn + zn + L1 , zn + L2 ? , 2 2 4 2 ? ? ? ? h2 L4 = f ? xn + h, yn + hzn + L2 , zn + hL3 ? . 2 ? ?9.6.3刚性方程组在求解方程组 ( 6.1) 时，经常出现解的分量数量级差别很大的情形，这给数值 求解带来很大的困难，这种问题称为刚性 stiff）问题 刚性（ 问题，在化学反应、电子网络 问题 刚性 和自动控制等领域中都是常见的，先考察以下例子． 给定系统u ' = ?10000.25u + 999.75v + 0.5, ? ? v ' = 999.75u ? 10000.25v + 0.5, ? ? u ( 0 ) = 1, ? ? v ( 0 ) = ?1. ?它可用解析方法求出准确解，方程右端系数矩阵 ? ?.75 ? A=? ? ? 999.75 ?1000.25 ?( 6.8)的特征值为 λ1 = ?0.5, λ2 = ?2000 ，方程的准确解为 ?u ( t ) = ?e ?0.5t + e?2000 t + 1, ? ? ?0.5 t ? e?2000 t + 1. ? v ( t ) = ?e ? 当 t → ∞ 时，u ( t ) → 1, v ( t ) → 1 称为稳态解，u , v 中均含有快变分量 e?2000t 及慢变分 量 ?e ?0.5t ．对应于 λ2 的快速衰减的分量在 t = 0.005 秒时已衰减到 e?10 ≈ 0 ，称τ2 = ?1λ2=1 = 0.0005 为时间常数 时间常数．当 t = 10τ 2 时快变分量即可被忽略，而对应 时间常数 2000 1 = 1 = 2 ，它要计算到 t = 10τ 1 = 20 时， 0.5于 λ1 的慢变分量，它的时间常数 τ 1 = ?λ1才能衰减到 e?10 ≈ 0 ， 也就是说解 u , v 必须计算到 t = 20 才能达到稳态解． 它表明方 程 ( 6.8) 的解分量变化速度相差很大，是一个刚性方程组．如果用四阶龙格―库塔 法求解．步长选取要满足 h & ?2.78 / λ ，即 h & ?2.78 / λ2 = 0.00139 ，才能使计算稳 定．而要计算到稳态解至少需要算到 t = 20 ，则需计算 14388 步．这种用小步长计 算长区间的现象是刚性方程数值求解出现的困难， 它是系统本身病态性质引起的． 对一般的线性系统dy = Ay ( t ) + g ( t ) , dtT T( 6.9 )其 中 y = ( y1 , L , yN ) ∈ R N , g = ( g1 , L , g N ) ∈ R N , A ∈ R N × N ． 若 A 的 特 征 值λ j = α j + i β j j = 1,L, N , i = ?1 相 应 的 特 征 向 量 为 ? j ( j = 1,L, N ) ， 则 方 程 组()( 6.9 ) 的通解为y ( t ) = ∑ c j e j ? j +ψ ( t ) ,λtj =1 N( 6.10 )其中 c j 为任意常数，可由初始条件 y ( a ) = y 0 确定，ψ ( t ) 为特解． 假定 λ j 的实部 α j = Re ( λ j ) & 0 ，则当 t → ∞ 时， y ( t ) → ψ ( t ) ,ψ ( t ) 为稳态解． 定义 8 且 若线性系统 ( 6.9 ) 中 A 的特征值 λ j 满足条件 Re ( λ j ) & 0 ( j = 1, L , N ) ，s = max Re ( λ j ) / min Re ( λ j ) && 1,1≤ j ≤ N 1≤ j ≤ N则称系统 ( 6.9 ) 为刚性方程 刚性方程，称 s 为刚性比 刚性比． 刚性方程 刚性比 刚性比 s && 1 时， A 为病态矩阵，故刚性方程也称病态方程．通常 s ≥ 10 就认 为是刚性的． s 越大病态越严重．方程 ( 6.8) 的刚性比 s = 4000 ，故它是刚性的． 对一般非线性方程组 ( 6.1) ，可类似定义 8 ，将 f 在点 ( t , y ( t ) ) 处展开，记 J (t ) = ?f ∈ R N × N ，假定 J ( t ) 的特征值为 λ j ( t ) , j = 1,L , N ，于是由定义 8 可知，当 ?yλ j ( t ) 满足条件 Re ( λ j ) & 0 ( j = 1,L , N ) ，且s ( t ) = max Re ( λ j ) / min Re ( λ j ) && 1,1≤ j ≤ N 1≤ j ≤ N则称系统 ( 6.1) 是刚性的， s ( t ) 称为方程 ( 6.1) 的局部刚性比． 求刚性方程数值解时，若用步长受限制的方法就将出现小步长计算大区间的 问题，因此最好使用对步长 h 不加限制的方法，如前面已介绍的欧拉后退法及梯 形法，即 A-稳定的方法，所谓稳定就是指数值方法的绝对稳定域包含了 ? = hλ 平 面的左半平面．这种方法当然对步长 h 没有限制，但 A-稳定方法要求太苛刻，已 证明所有显式方法都不是 A-稳定的， 而隐式的 A-稳定多步法阶数最高为 2 ， 且以 梯形法误差常数为最小．这就表明本章所介绍的方法中能用于解刚性方程的方法 很少．通常求解刚性方程的高阶线性多步法事吉尔 Gear）方法 吉尔（ 方法，还有隐式龙格吉尔 方法 库塔法（见文献 [ 21] ） ，这些方法都有现成的数学软件可供使用．本书不再介绍．评注本章研究求解常微分方程初值问题的数值方法． 构造数值方法主要有两条途 径：基于数值积分的构造方法和基于泰勒展开的构造方法．后一种方法更灵活， 也更具有一般性．泰勒展开方法还有一个优点，它在构造差分公式的同时可以得 到关于截断误差的估计． 基于泰勒展开构造出的四阶龙格-库塔方法（见 3.3 节）则是计算机上的常用 算法．四阶龙格-库塔方法的优点是精度高，程序简单，计算过程稳定，并且易于调节步长． 四阶龙格-库塔方法也有不足之处：它要求函数 f ( x, y ) 具有较高的光滑 性．如果 f ( x, y ) 的光滑性差，那么，它的精度可能还不如欧拉公式（ 2.1 节）或 改进的欧拉公式（ 2.4 节） ． 四阶龙格-库塔方法的另一个缺点是计算量比较大，需要耗费较多的机器时 间（每一步需四次计算函数 f ( x, y ) 的值） ．相比之下，汉明方法（ 5.4 节）可以节 ．但汉明方法是一种四步法， 省计算量（每一步只需两次计算函数 f ( x, y ) 的值） 它不是自开始的，需要借助于四阶龙格-库塔方法提供开始值． 对数值方法的分析还涉及到局部截断误差，整体误差，收敛性，相容性及稳 定性等概念，特别是有关绝对稳定性的讨论，它涉及计算时步长的选取，如步长 选的不合适，舍入误差恶性增长，结果完全错误．本章只对单步法的收敛性和稳 定性做了讨论．对多步法的有关理论没有涉及，读者可参阅文献 [ 2] 及 [ 21] ． 刚性方程组是具有重要应用价值的问题，它的理论和解法很多，并且还在进 一步发展，详细介绍可参阅文献 [ 21] 和 [ 22] ．习1 ． 用欧拉法解初值问题题y ' = x 2 + 100 y 2 , y ( 0 ) = 0.取步长 h = 0.1 ，计算到 x = 0.3 （保留到小数点后 4 位）2 ．用改进的欧拉法和梯形法解初值问题y ' = x 2 + x ? y, y ( 0 ) = 0.取步长 h = 0.1 ，计算到 x = 0.5 ，并与准确解 y = ?e ? x + x 2 ? x + 1 相比较．3 ．用梯形法解初值问题? ' ? y + y = 0, ? ? y ( 0 ) = 1. ?证明其近似解为? 2?h? yn = ? ? , ? 2+h?n并证明当 h → 0 时，它收敛于原初值问题的准确解 y = ?e ? x ．4 ．利用欧拉方法计算积分∫在点 x = 0.5,1,1.5, 2 的近似值．x0et25 ．取 h = 0.2 ，用四阶经典龙格―库塔方法求解下列初值问题：1)? ' ? y = x + y, ? ? y ( 0 ) = 1; ?0 & x & 1;2)? y ' = 3 y / (1 + x ) , ? ? ? y ( 0 ) = 1. ?0 & x & 1;6 ．证明对任意参数 t ，下列龙格―库塔公式是二阶的： h ? yn +1 = yn + ( K 2 + K 3 ) , ? 2 ? ? K1 = f ( xn , yn ) , ? ? K 2 = f ( xn + th, yn + thK1 ) , ? ? K 3 = f ( xn + (1 ? t ) h, yn + (1 ? t ) hK1 ) . ? 7 ．证明中点公式 ( 3.10 ) 是二阶的，并求其绝对稳定区间． 8 ．对于初值问题y ' = ?100 ( y ? x 2 ) + 2 x, y ( 0 ) = 1. (1) 用欧拉法求解，步长 h 取什么范围的值，才能使计算稳定． (2) 若用四阶龙格―库塔法计算，步长 h 如何选取？ (3) 若用梯形公式计算，步长 h 有无限制． 9． 分别用二阶显式阿当姆斯方法和二阶隐式阿当姆斯方法解下列初值问题：y ' = 1 ? y, y ( 0 ) = 0.取 h = 0.2, y0 = 0, y1 = 0.181 ，计算 y (1.0 ) 并与准确解 y = 1 ? e ? x 相比较．10 ．证明解 y ' = f ( x, y ) 的下列差分公式yn +1 =1 h ( yn + yn+1 ) + ( 4 y 'n +1 ? y 'n + 3 y 'n ?1 ) 2 4是二阶的，并求出截断误差的主项．11 ．试证明线性二步法yn + 2 + ( b ? 1) yn +1 ? byn =h ?( b + 3) f n + 2 + ( 3b + 1) f n ? ? 4?当 b ≠ ?1 时方法为二阶，当 b = ?1 时方法为三阶．12 ．求方程?u ' = ?10u + 9v, ? ? ' ?v = 10u ? 11v. ?的刚性比，用四阶 R ? K 方法求解时，最大步长能取多少？
第九章.常微分方程初值问题数值解法―汇集和整理大量word文档,专业文献,应用文书,考试资料,教学教材,办公文档,教程攻略,文档搜索下载下载,拥有海量中文文档库,关注高价值的实用信息,我们一直在努力,争取提供更多下载资源。