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

第9章 常微分方程初值问题数值解法 9.2 简单的数值方法与基本概念 9.2.2 梯形方法 9.2.3 单步法的局部截断误差与阶 9.3 龙格—库塔方法 9.3.2 二阶显式R-K方法 9.3.3 三阶与四阶显式R-K方法 　　对欧拉法E(hλ)=1+hλ，其绝对稳定域为|1+hλ|1,绝对稳定区间为-2λ0,在例5中λ=-100,-2-100h0,即0h2/100=0.02为稳定区间，在例4中取h=0.025，故它是不稳定的，当取h=0.005时它是稳定的. 对二阶R-K方法，解模型方程(4.1)可得到 故 绝对稳定域由|E(hλ)|1得到，于是可得绝对稳定区间为-2hλ0，即0h2/λ. 对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)处做泰勒展开，由于f(x, y )是二元函数，故要用二元泰勒展开，各项展开式为 将以上结果代入(3.7)，则有 要使公式(3.6)具有p=2阶，必须使 即 (3.9)的解是不唯一的. 可令c2=a≠0，则得 这样得到的公式称为二阶R-K方法. 则由此可以看出在改进的欧拉公式中相当于取(xn,yn), (xn+1,yn+1)两点处斜率的平均值，近似代替平均斜率，其精度比欧拉公式提高了. 如取a=1/2，则c1= c2=1/2, λ2=μ21=1. 这就是改进的欧拉公式(3.1). 称为中点公式(变形的欧拉公式)，相当于数值积分的中矩形公式.也可以表示为 如取a=1，则c1=0, c2=1, λ2=μ21=1/2. 得计算公式 对r=2的R-K公式(3.6)能否使局部误差提高到O(h4)? 为此 需把K2多展开一项，从(3.8)的 看到展开式中的项 是不能通过选择参数削掉的，实际上要使 h3 的项为零，需增加3个方程，要确定4个参数c1, c2, λ2, μ21，这是不可能的. 故r＝2的显式R-K方法的阶只能是p=2，而不能得到三阶公式. 要得到三阶显式R-K方法，必须r=3. 此时计算(3.4), (3.5)的公式表示为 其中c1, c2, c3及λ2, μ21, λ3, μ31, μ32均为待定常数，公式(3.11)的局部截断误差为 只要K1, K2将按二元泰勒展开，使Tn+1＝O(h4)，可得待定参数满足方程 这是8个未知数6个方程的方程组，解不是唯一的. 可以得到很多公式. 满足条件(3.12)的公式(3.11)统称为三阶R-K公式. 下面只给出其中一个常见的公式. 此公式称为库塔三阶方法. 继续上述过程，经过较复杂的数学演算，可以导出各种四阶R-K公式，下列经典公式是其中常用的一个： 四阶R-K方法的每一步需要计算四次函数值f，可以证明其局部截断误差为O(h5).　(例3见书p349) 然而值得指出的是，龙格-库塔方法的推导基于泰勒展开方法，因而它要求所求的解具有较好的光滑性质. 反之，如果解的光滑性差，那么，使用龙格-库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法. 实际计算时，我们应当针对问题的具体特点选择合适的算法. 9.3.4 变步长的龙格-库塔方法 单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加了. 步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累. 因此同积分的数值计算一样，微分方程的数值解法也有个选择步长的问题. 在选择步长时，需要考虑两个问题： 1. 怎样衡量和检验计算结果的精度？ 2. 如何依据所获得的精度处理步长？ 　　我们考察四阶R-K公式(3.13) ，从节点xn出发，先以h为步长求出一个近似值，记为 ，由于公式的局部截断误差为O(h5)，故 然后将步长折半，即取为步长 ,从xn跨两步到xn+1，再求得一个近似值 ，每跨一步的局部截断误差是 ，因此有 比较(3.14)式和(3.15)式我们看到，步长折半后，误差大约减少到1/16，即有 由此易得下列事后估计式 这样，我们可以通过检查步长，折半前后两次计算结果的偏差 来判定所选的步长是否合适，具体地说