第七章常微分方程数值解法

一、

基本思想7.3龙格-库塔(Runge-Kutta)方法1.从积分的角度考虑.2.从微分的角度考虑.应用微分中值定理，有，利用所给微分方程，有其中K*=f(ξ,

y(ξ))

称为区间[xn,xn+1]上的平均斜率。(7)对差商这样，只要对平均斜率K*提供一种算法，由式按照这种观点考察Euler格式，作为平均斜率K*的近似，精度自然很低。便相应地导出一种计算格式。它简单地取点xn的斜率可将改进的欧拉格式改写成的算术平均值作为平均斜率。该公式可以看作是用xn和xn+1两个点处的斜率和由改进型欧拉公式我们可以猜想，如果在[xn,xn+1]内多预测几个点的斜率，再对他们进行加权平均，可能得到精度更好的平均斜率！再考察改进的Euler格式。二、下面以2阶龙格-库塔方法为例来阐述这种思想考察区间[xn,xn+1]

上的一点xn+α

=xn+αh,0<α≤1,用xn和xn+α

的斜率K1和K2的加权平均作为平均斜率K*的近似值：即取其中和是待定常数。若取k1=f(xn,yn)，则问题在于如何确定xn+α

处的斜率k2和常数λ1和λ2。K*=λ1K1+λ2K2yn+1=yn+h(λ1K1+λ2K2)仿照改进的欧拉方法，用欧拉方法预测y(xn+α)的值，并用它来估计斜率K2：于是得到如下形式的算法：通过适当选取参数λ1,

λ2和α的值，使得公式具有2阶精度！！yn+α

=yn+αhK1K2=f(xn+α,yn+α

)(8)式(8)是为二级二阶Runge-Kutta方法的一般形式由泰勒公式展开，要使公式具有2阶精度，只需方程组有无穷多解：二级二阶R-K方法有无穷多种取取（中点法）（改进的Euler法）三、三级三阶R-K方法形式：（9）我们要确定参数使得（9）式的局部截断误差为O(h4).将K2,K3按二元函数泰勒公式展开，代入(9)中第一式，再将y(xn+1)在xn处作泰勒展开，最后通过比较两个展开式的系数，使局部截断误差为O(h4).得到如下方程组：此为5个方程，但含有7个未知量的方程组，它有无穷多个解。

库塔三阶方法四、四级四阶R-K方法：局部截断误差经典龙格-库塔方法解：例2：用经典的龙格-库塔方法求解下列初值问题。经典的四阶龙格-库塔公式：0.10.20.30.40.51.09541.18321.26491.34161.41420.60.70.80.91.01.48321.54921.61251.6733

1.7321同保留5位的精确值完全一致：0.10.20.30.40.51.09541.18321.26491.34161.41420.60.70.80.91.01.48321.54921.61251.6733

1.7321

0.11.1000

1.0959

1.0954

0.21.1918

1.1841

1.18320.31.2774

1.2662

1.26490.41.3582

1.3434

1.34160.51.4351

1.4164

1.41420.61.5090

1.4860

1.48320.71.5803

1.5525

1.54920.81.6498

1.6165

1.61250.91.7178

1.6782

1.6733

1.01.7848

1.7379

1.7321

7.4单步法的收敛性和稳定性收敛性和稳定性从不同角度描述了数值方法的可靠性.只有既收敛又稳定的方法，才能提供比较可靠的计算结果。一、收敛性定义若一种数值方法对任意固定的xn=x0+nh,当h→0时，有yn→y(xn),则称该数值方法是收敛的。例设考察Euler方法求解此问题的收敛性。该问题的精确解为Euler方法计算公式为从而有当h→0时,考察在固定点x*=xn=x0+nh近似值的变化。开始取步长h=h0,x*=x2=x0+2h再取步长h=h0∕2,x*=x4=x0+4h进一步取步长h=h0∕4,x*=x8=x0+8h,

……x0

x1

x2x0

x4

x0

x8

这样，是固定的。当h→0

(即n→∞)时对于此例，Euler方法是收敛的。再考察隐式Euler方法。此问题的隐式Euler格式为这时有从而数值解这时当h→0时仍然有yn→y(xn).一种方法的收敛性不能只针对某一个具体实例的成立与否，而是必须考虑对所有初值问题都收敛。下面讨论一般单步法的收敛性。首先，关于求解初值问题的所有显式单步方法都可以写成如下形式其中h为步长，称为增量函数,它依赖于方程右端函数f,且不同的单步方法，φ的形式不同。(1)(2)Euler方法：改进的Euler方法：单步法式（2）在节点xn+1的局部截断误差定义为注意，局部截断误差en+1是在yn=y(xn)精确成立的假定下，从xn出发，用式（2）向前推进一步到xn+1而产生的误差，并不是在点xn+1处，由于使用式（2）求解所引进的全部误差。(3)事实上，用式（2）计算，每推进一步都会产生误差，因此近似值yn包含有前n步的累计误差。把从x0处的初值y0出发，用单步法（2）推进n+1步，到点xn+1所得到的近似值yn+1与其精确值y(xn+1)的偏差En+1=y(xn+1)-yn+1称为整体截断误差。由上述定义可知，一种单步法是否收敛，就是看其整体截断误差En=

y(xn)-yn在h→0

(n→∞)时是否有En→0定义(4)为了估计单步法式（2）的整体截断误差En，我们需要下面的引理。引理设α>0,

β≥0为实数，且实数序列{ηn}满足递推关系则有证利用不等式eα>1+α,用数学归纳法（略）。定理若初值问题（1）的一个单步方法（2）的局部截断误差为O(hp+1)(p>=1)，且（2）式中的增量函数φ(x,y,h)满足对y的Lipschitz条件，即存在L>0使得对一切有则单步法（2）的整体截断误差为En=O(h

p

).从而单步法收敛。(5)证依据局部截断误差的定义，有用上式减去（2）式，得根据定理中的假设条件，有再利用引理的结论可得到（C为正常数）假定初值是准确的，即E0=

y(a)-y0=0则有，从而有当h→0时，En→0.依据这个定理，判断一个单步法（2）的收敛性，就归结为验证其增量函数φ(x,y,h)是否对y满足Lipschitz条件。推论1Euler方法是收敛的。推论2

改进的Euler方法是收敛的。推论3各级Runge-Kutta方法是收敛的。在以上讨论中，(2)是p阶的方法，即en+1=O(hp+1),一般都是指p至少等于1，即p+1>=2.若将en+1按变量h在h=0作Taylor展开，得到很明显，en+1=O(hp+1)，而p>=1的充要条件是二、相容性(6)定义若方法（2）增量函数满足则称单步法（2）与初值问题（1）相容。据以上讨论，与初值问题(1)相容的方法至少是一阶的。

对于与初值问题相容的单步方法(2)，若φ(x,y,h)

满足对y的Lipschitz条件，则由前面的定理，此方法是收敛的。(7)若把（2）写成(8)并令x

=

xn

(=

x0+nh)固定，h→0.因为`yn→y(xn),(8)式可写成即(9)这意味着h→0时，数值计算格式(8)趋于微分方程(9).这就是“相容”一词的意义。不难验证，已讨论过的各种方法都是与所对应的初值问题相容的三、稳定性定义

设用某一数值方法计算yn时，所得到的实际计算结果为ŷn,且由误差δn=yn-ŷn引起以后各点处ym(m>n)的误差为δm,如果总有|δm|≤|δn|，则称该数值方法是绝对稳定的。对于微分方程若Reλ>0，它是不稳定的（这里指微分方程本身的性质）若Reλ<0，它是稳定的.但用不同的数值方法来求解时，还有数值稳定与数值不稳定的两种可能。所以我们考虑Reλ<0的情形，对于实数λ，考虑λ<0。当λ为实数时，Euler方法稳定条件：-2<λh<0改进的Euler方法