计算方法4常微分方程数值解法.pptVIP

然后可得到有限差分公式 上式可写成 写成矩阵形式 其中 上述方程组有四个未知量，只有三个方程，有无穷多组解。 取任意一组解便得一种二阶龙－库公式。 当c1=c2=1/2, a2=b21=1时二阶Runge-Kutta公式为 yn+1=yn+k1/2+k2/2 k1=hf(xn,yn) k2=hf(xn+h,yn+k1) 此即改进Euler法 取c2=0 ，c2=1，a2=1/2，b21=1/2 yn+1=yn+k2 k1=hf(xn,yn) k2=hf(xn+h/2,yn+k1/2) 此为中点法或变形的 Euler公式 三阶龙格－库塔法是用三个值k1,k2,k3的加权平均来近似k*， 即有 yn+1=yn+c1k1+c2k2+c3k3 k1=hf(xn,yn) k2=hf(xn+a2h,yn+b21k1) k3=hf(xn+a3h,yn+b31k1+b32k2) 要使其具有三阶精度，必须使局部截断误差为O(h4) 类似二阶龙格－库塔法的推导，c1,c2,c3,a2,a3,b21,b31,b32应满足 c1+c2+c3=1 a2=b21 a3=b31+b32 c2a2+c3a3=1/2 c2a22+c3a32=1/3 c3a32=1/6 由该方程组任意解可得三阶龙格-库塔公式 例：Kutta公式 kn+1=yn+(k1+4k2+k3)/6 k1=hf(xn,yn) k2=hf(xn+h/2,yn+k1/2) k3=hf(xn+h,yn-k1+2k2) 类似可推出四阶龙格-库塔公式，常用的有 例：经典Runge-Kutta法 yn+1=yn+(k1+2k2+2k3+k4)/6 k1=hf(xn,yn) k2=hf(xn+h/2,yn+k1/2) k3=hf(xn+h/2,yn+k2/2) k4=hf(xn+h,yn+k3) 局部截断误差 O(h5) 还有: Gill公式及m (m4)阶龙格－库塔法。 m4时：计算量太大，精确度不一定提高，有时会降低。 Gill公式 节省存储单元 控制舍入误差 对于经典的四阶Runge-Kutta法给出如下算法： [算法4.2]求解： dy/dx=f(x,y) a≤x≤b y (a)=y0 Step 1: 输入a,b,y0 及N Step 2: (b-a)/N=h,a=x,y0=y Step 3: 输出 (x,y) Step 4: For I=1 T0 N hf(x,y)=k1 hf(x+h/2,y+ k1/2)= k2 hf(x+h/2,y+k2/2)=k3 hf(x+h,y＋k3)=k4 y+(k1+2k2+2k3+k4)/6=y x+h=x 输出(x,y) Step 5 : 停止 [例4.3]用四阶经典Runge－Kutta方法解初值问题： (1)求 ， (2)求 ， 自适应龙格－库塔法 用户提出问题I ： 问题：①：如何判断|y(xn)-yn|ε ∵精度值y(xn)未知。 ②：如何取h=? 解①：如用p阶龙格－库塔法计算，局部截断误差为O(hp+1) y’=f(x,y) y(x0)=y0 要求误差ε＝10-8 求数值解 xn h/2 xn+1 h 如 xn-----------------?xn+1 令 yn=y(xn) yn+1(h) 则 y(xn+1)-yn+1(h)≈chp+1 步长折半xn?xn+h/2?xn+1分两步计算y(xn+1)的近似值yn+1(h/2)。 则 y(xn+1)-yn+1(h/2)≈2c(h/2)p+1 定理：对于问题I若用P阶龙格－库塔法计算y(xn+1)在步长折半前后的近似值分别为yn+1(h), yn+1(h/2)则有误差公式 注：10 误差的事后估计法 20 停机准则：△ε (可保证|y(xn+1)-yn+1(h/2)|ε) 解②： ⑴ h取大，局部截断误差chp+1大，不精确 ⑵ h取小，运算量大（步多），舍入误差积累大 解决策略：变步长龙格－库塔法 If（△ε） 将步长折半反复计算，直至△ε为止， 取h为最后一次的步长， yn+1为最后一次计算的结果。 Else if (△