计算方法5.常微分方程的数值解法.ppt

构造高阶单步法的直接方法 由Taylor公式： 当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程： 其局部截断误差为： Runge?Kutta方法 1. 基本思想 因为 其中K?=f(?,y(?))称为y(x)在[xi,xi+1]上的平均斜率。 若取K1=f(xi,y(xi))——Euler公式 取K2=f(xi+1,y(xi+1))——向后Euler公式 一阶精度 取 —— 梯形公式 二阶精度 猜想：若能多预测几个点的斜率，再取其加权平均作为K?，可望得到较高精度的数值解，从而避免求f 的高阶导数。 2. R?K公式 其中Kj为y = y(x)在xi + ajh (0 ? aj ? 1)处的斜率预测值。 aj，bjs，cj为特定常数。 3. 常数的确定 确定的原则是使精度尽可能高。 以二阶为例： (希望y(xi+1) – yi+1 = O(hp)的阶数p尽可能高) 首先： 另一方面： 将K2在(xi，yi)处展开。 K2 = f (xi,yi) + a2hf x(xi,yi) + b21hK1 f y(xi,yi) + O(h2) 可得： yi+1 = yi + hc1 f (xi,yi) + hc2 f (xi,yi) + h2c2[a2 f x(xi,yi) + b21K1 f y(xi,yi)] + O(h3) = yi + h(c1 + c2) f (xi,yi) + c2a2h2[f x(xi,yi) + (b21/a2) f (xi,yi) f y(xi,yi)]+O(h3) （希望） 希望：ei+1 = y(xi+1) – yi+1 = O(h3). 则应： 特例：a2 = 1 ? c1 = c2 = 1/2，b21 = 1，得2阶R-K公式： ??改进欧拉公式。 希望：ei+1 = y(xi+1) – yi+1 = O(h3). 则应： 特例：c1 = 0 ? c2 = 1，a2 = 1/2，b21 = 1/2，得： 称为中点公式。 4. 最常用的R-K公式 —— 标准4阶R-K公式 算法： Matlab代码： 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+2*K3+K4)/6; y0=y1 end; end; function f=f(x,y) f=x+y; end; 例：用标准4阶R-K公式求： 的数值解。取h = 0.2，并与标准解y = 2ex – x – 1比较。 解：因为f(x，y) = x + y，从由得： 1.Runge-Kutta 法的一般形式 2. 2级2阶Runge-Kutta方法 式(9.11) 称L级p阶Runge-Kutta方法(简称R-K法)。 当L＝1就是欧拉法，当L＝2时为改进的欧拉法。 3. 经典Runge-Kutta方法 我们可以构造出一族3级3阶，一族4级4阶和一族5级4 阶等R-K方法。最常用的4级4阶是如下的经典R-K方法： 4．R-K-Fehhlberg 方法 R-K-Fehhlberg方法是在R-K方法的基础上引进误差和步长 控制的办法。即利用5阶R-K法 5. 隐式R-K方法 类似于显式R-K公式，稍加改变，就得到隐式R-K方法。 通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式 R-K法有: 6.变步长方法 在单步法中每一积分步步长实际上是相互独立的，步长的 选择具有了灵活性。根据合理地选择每一积分步的步长， 既保证精度的要求，又可以减少计算量，从而减小舍入误 差。其方便的控制手段是基于误差的事后估计式。 * 5.3 龙格 - 库塔法 /* Runge-K