计算方法5常微分方程的数值解法概要课件

1、5.3 龙格 - 库塔法 /* Runge-Kutta Method */ 考察改进的欧拉法，可以将其改写为：斜率一定取K1 K2 的平均值吗？步长一定是一个h 吗？单步递推法的基本思想是从 ( xi , yi ) 点出发，以某一斜率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所能达到的最高精度为2阶。建立高精度的单步递推格式。10/15/202215.3 龙格 - 库塔法 /* Runge-Kutta M首先希望能确定系数 1、2、p，使得到的算法格式有2阶精度，即在 的前提假设下，使得 Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开将改

2、进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+llStep 2: 将 K2 代入第1式，得到10/15/20222首先希望能确定系数 1、2、p，使得到的算法格式有2阶精Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较要求 ，则必须有：这里有 个未知数， 个方程。32存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库塔格式。注意到， 就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？10/15/20223Step 3: 将 yi+1 与 y( xi+1 ) 在 x构造高阶单步法的直接方

3、法由Taylor公式：当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程： 其局部截断误差为：差分方程为p阶方式，上述方式称为Taylor方式。注：利用Tayler公式构造，不实用，高阶导数f (i)不易计算。10/15/20224构造高阶单步法的直接方法差分方程为p阶方式，上述方式称为Ta RungeKutta方法1. 基本思想因为 其中K=f(,y()称为y(x)在xi,xi+1上的平均斜率。若取K1=f(xi,y(xi)Euler公式取K2=f(xi+1,y(xi+1)向后Euler公式 一阶精度取 梯形公式 二阶精度猜想：若能多预

4、测几个点的斜率，再取其加权平均作为K，可望得到较高精度的数值解，从而避免求f 的高阶导数。10/15/20225 RungeKutta方法10/11/202252. RK公式 其中Kj为y = y(x)在xi + ajh (0 aj 1)处的斜率预测值。aj，bjs，cj为特定常数。10/15/202262. RK公式10/11/202263. 常数的确定确定的原则是使精度尽可能高。以二阶为例： (希望y(xi+1) yi+1 = O(hp)的阶数p尽可能高)首先：10/15/202273. 常数的确定10/11/20227另一方面：将K2在(xi，yi)处展开。K2 = f (xi,yi)

5、+ a2hf x(xi,yi) + b21hK1 f y(xi,yi) + O(h2)可得：yi+1 = yi + hc1 f (xi,yi) + hc2 f (xi,yi) + h2c2a2 f x(xi,yi) + b21K1 f y(xi,yi) + O(h3) = yi + h(c1 + c2) f (xi,yi) + c2a2h2f x(xi,yi) + (b21/a2) f (xi,yi) f y(xi,yi)+O(h3)（希望）10/15/20228另一方面：10/11/20228希望：ei+1 = y(xi+1) yi+1 = O(h3). 则应：特例：a2 = 1 c1 =

6、c2 = 1/2，b21 = 1，得2阶R-K公式： 改进欧拉公式。10/15/20229希望：ei+1 = y(xi+1) yi+1 = O(h希望：ei+1 = y(xi+1) yi+1 = O(h3). 则应：特例：c1 = 0 c2 = 1，a2 = 1/2，b21 = 1/2，得： 称为中点公式。10/15/202210希望：ei+1 = y(xi+1) yi+1 = O(h4. 最常用的R-K公式 标准4阶R-K公式 算法：输入a，b，n，y0h=(b-a)n，x0 = afor i = 1, i=n, i+K1 = f(x0, y0)K2 = f(x0+h/2, y0+h*K1/

7、2)K3 = f(x0+h/2, y0+h*K2/2)K4 = f(x0+h, y0+h*K3)x0 = x0+hy0 = y0 + h*(K1+2*K2+2*K3+K4)/6输出x0，y010/15/2022114. 最常用的R-K公式 标准4阶R-K公式输入a，bMatlab代码：function Runge_Kutta4(a,b,h,y0)n=(b-a)/h;x0=a;for i=1:n K1=f(x0,y0) K2=f(x0+h/2,y0+h*K1/2) K3=f(x0+h/2,y0+h*K2/2) K4=f(x0+h,y0+h*K3) x0=x0+h y1=y0+h*(K1+2*K2

8、+2*K3+K4)/6; y0=y1end;end;function f=f(x,y) f=x+y;end;输入a，b，n，y0h=(b-a)n，x0 = afor i = 1, i ,反复将步长折半进行计算，直至 为止,这时取最终得到的作为结果；如果 为止，这时再将步长折半一次，就得到所要的结果。这种通过加倍或折半处理步长的计算方法称为变步长方法。 注：推荐使用精度好计算量低的变步长方法。用四阶显式R-K方法做变步长方法是实践中较好的方法！ 10/15/2022256.变步长方法在单步法中每一积分步步长实际上是相互独立的，步例 分别用改进的欧拉格式和四阶龙格库塔格式解初值问题（取步长h=0.2）：10/15/202226例 分别用改进的欧拉格式和四阶龙格库塔格式解初值问题节点 改进欧拉法 四阶龙格库塔法 准确解 0 1 1 1 0.2 1.186667 1.183229 1.1832160.4 1.348312 1.341667 1.3416410