数值分析7.2龙格—库塔方法.ppt

*7.3.4 变步长的龙格-库塔方法 单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加了. 步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累. 因此同积分的数值计算一样，微分方程的数值解法同样存在选择步长的问题. 在选择步长时，需要考虑两个问题： 1. 怎样衡量和检验计算结果的精度？ 2. 如何依据所获得的精度处理步长？ 　　我们考察四阶R-K公式(3.13) ，从节点xn出发，先以h为步长求出一个近似值，记为 ，由于公式的局部截断误差为O(h5)，故 然后将步长折半，即取为步长 ,从xn跨两步到xn+1，再求得一个近似值 ，每跨一步的局部截断误差是 ，因此有 比较(3.14)式和(3.15)式我们看到，步长折半后，误差大约减少到1/16，即有 由此易得下列事后估计式 这样，我们可以通过检查步长，折半前后两次计算结果的偏差 来判定所选的步长是否合适，具体地说，将区分以下两种情况处理： 　　这种通过加倍或折半处理步长的方法称为变步长方法.表面上看，为了选择步长，每一步的计算量增加了，但总体考虑往往是合算的. 　　1.对于给定的精度ε,如果Δε,我们反复将步长折半计算,直至Δε为止,这时取最终得到的 作为结果； 　　2.如果Δε，我们将反复将步长作加倍计算，直至Δε为止，这时再将步长折半计算一次，就得到所要的结果. 上页 下页 7.3 龙格—库塔方法 对许多实际问题来说，欧拉公式与改进欧拉公式精度还不能满足要求，为此从另一个角度来分析这两个公式的特点，从而探索一条构造高精度方法的途径. 7.3.1 显式龙格—库塔法的一般形式 上节给出了显式单步法的表达式(2.11), 其局部截断误差为O(hp+1)，对欧拉法Tn+1=O(h2)，即方法为p=1阶，若用改进欧拉法(2.9)： 此时增量函数为 它比欧拉法的?(xn, yn, h)=f(xn, yn), 增加了计算一个函数 f 的值，期望 p=2.若要使得到的公式阶数 p 更大, ? 就必须包含更多的 f 值. 实际上从方程(1.1)等价的积分形式(2.4) ，即 若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加求积节点，为此可将(3.3)的右端用求积公式表示为 一般说来，点数 r 越多，精度越高，上式右端相当于增量函数? (xn, yn, h)，为得到便于计算的显式方法，可类似于改进欧拉法(3.1),(3.2)，将公式表示为 其中 这里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)是其中一种，下面将证明其阶 p=2. 要使公式(3.4),(3.5)具有更高的阶 p，就要增加点数 r. 下面我们只就 r=2推导R-K方法. 并给出 r=3,4 时的常用公式，其推导方法与 r=2 时类似，只是计算较复杂. 补充：二元函数的泰勒公式 一元函数 的泰勒公式: 推广 多元函数泰勒公式 记号： (设下面涉及的偏导数连续): 一般地, 表示 表示 定理： 的某一邻域内有直 到 n + 1 阶连续偏导数 , 为此邻域内任 一点, 则有 其中 ① ② ① 称为f 在点(x0 , y0 )的 n 阶泰勒公式, ②称为其拉格 朗日型余项 . 7.3.2 二阶显式R-K方法 对r=2的R-K方法，由(3.4),(3.5)式可得如下计算公式 这里 c1, c2, λ2, μ21 均为待定常数，我们希望适当选取这些系数，使公式阶数 p 尽量高. 根据局部截断误差定义，推导出(3.6)的局部截断误差为 其中 这里yn=y(xn), yn+1=y(xn+1). 为得到Tn+1的阶 p，要将上式各项在(xn, yn)处做泰勒展开，用二元泰勒展开，各项展开式为 将以上结果代入(3.7)，则有 要使公式(3.6)具有 p=2 阶，必须使 即 (3.9)的解是不唯一的. 可令c2=a≠0，则得 这样得到的公式称为二阶R-K方法. 如取a=1/2，则c1= c2=1/2, λ2=μ21=1. 这就是改进的欧拉公式(3.1). 称为中点公式(变形的欧拉公式)，相当于数值积分的中矩形公式.也可以表示为 如取a=1，则c1=0, c2=1, λ2=μ21=1/