资源预览内容
第1页 / 共24页
第2页 / 共24页
第3页 / 共24页
第4页 / 共24页
第5页 / 共24页
第6页 / 共24页
第7页 / 共24页
第8页 / 共24页
第9页 / 共24页
第10页 / 共24页
亲,该文档总共24页,到这儿已超出免费预览范围,如果喜欢就下载吧!
资源描述
1 第三节 Runge-Kutta方法 2 2009, Henan Polytechnic University 23 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 3 龙格 - 库塔法 /* Runge-Kutta Method */ 考察改进的欧拉法,可以将其改写为: 斜率 一定取K1 K2 的平均值吗? 步长一定是一个h 吗 ? 单步递推法的基本思想是从 ( xi , yi ) 点出发,以某一斜 率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所 能达到的最高精度为2阶。 建立高精度的单步递推格式。 3 2009, Henan Polytechnic University 33 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 首先希望能确定系数 1、2、p,使得到的算法格式有2阶 精度,即在 的前提假设下,使得 Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开 将改进欧拉法推广为: ),( ),( 12 1 22111 phKyphxfK yxfK KKhyy ii ii ii += = += + Step 2: 将 K2 代入第1式,得到 4 2009, Henan Polytechnic University 43 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较 要求 ,则必须有: 这里有 个未知 数, 个方程。 3 2 存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库 塔格式。注意到, 就是改进的欧拉法。 Q: 为获得更高的精度,应该如何进一步推广? 5 2009, Henan Polytechnic University 53 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 其中i ( i = 1, , m ),i ( i = 2, , m ) 和 ij ( i = 2, , m; j = 1, , i1 ) 均为待定 系数,确定这些系数的 步骤与前面相似。 ).,( . ),( ),( ),( . 112211 23213133 12122 1 22111 + += += += = += mm mmmmim ii ii ii mmii hKhKhKyhxfK hKhKyhxfK hKyhxfK yxfK KKKhyy 最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ : 6 2009, Henan Polytechnic University 63 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 注: 龙格-库塔法的主要运算在于计算 Ki 的值,即计算 f 的 值。Butcher 于1965年给出了计算量与可达到的最高精 度阶数的关系: 753 可达到的最高精度 642 每步须算Ki 的个数 由于龙格-库塔法的导出基于泰勒展开,故精度主要受 解函数的光滑性影响。对于光滑性不太好的解,最好 采用低阶算法而将步长h 取小。 深入研究龙格- 库塔法请看此处! 7 2009, Henan Polytechnic University 73 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 7.2 RungeKutta方法 7.2.1 构造高阶单步法的直接方法 由Taylor公式: 当h充分小时,略去Taylor公式余项,并以yi、yi+1分别代替y(xi)、y(xi+1), 得到差分方程: (7.2-1) 其局部截断误差为: 8 2009, Henan Polytechnic University 83 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 当h充分小时,略去Taylor公式余项,并以yi、yi+1分别代替y(xi)、y(xi+1), 得到差分方程: (7.2-1) 其局部截断误差为: 即(7.2-1)为p阶方式,上述方式称为Taylor方式。 注:利用Tayler公式构造,不实用,高阶导数f (i)不易计算。 9 2009, Henan Polytechnic University 93 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 7.2.2 RungeKutta方法 1. 基本思想 因为 = y(xi) + hf (,y() = y(xi) + hK 其中K = f (,y()称为y(x)在xi,xi+1上的平均斜率。 若取 K1 = f (xi,y(xi) Euler公式 取 K2 = f (xi+1,y(xi+1) 向后Euler公式 一阶精度 取 梯形公式 二阶精度 猜想:若能多预测几个点的斜率,再取其加权平均作为K,可望得到较高精 度的数值解,从而避免求f 的高阶导数。 10 2009, Henan Polytechnic University 103 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 2. RK公式 (7.2-4) 其中Kj为y = y(x)在xi + ajh (0 aj 1)处的斜率预测值。 aj,bjs,cj为特定常数。 11 2009, Henan Polytechnic University 113 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 3. 常数的确定 确定的原则是使精度尽可能高。 以二阶为例: (7.2-5) (希望y(xi+1) yi+1 = O(hp)的阶数p尽可能高) 首先: 12 2009, Henan Polytechnic University 123 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 另一方面: 将K2在(xi,yi)处展开。 K2 = f (xi,yi) + a2hf x(xi,yi) + b21hK1 f y(xi,yi) + O(h2). 代入(7.2-5)得: 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/a12) f (xi,yi) f y(xi,yi)+O(h3) (希望) 13 2009, Henan Polytechnic University 133 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 希望:ei+1 = y(xi+1) yi+1 = O(h3). 则应: 特例:a2 = 1 c1 = c2 = 1/2,b21 = 1,得2阶R-K公式: 改进欧拉公式。 14 2009, Henan Polytechnic University 143 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 希望:ei+1 = y(xi+1) yi+1 = O(h3). 则应: 特例:c1 = 0 c2 = 1,a2 = 1/2,b21 = 1/2,得: (7.2-7) 称为中点公式。 15 2009, Henan Polytechnic University 153 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 4. 最常用的R-K公式 标准4阶R-K公式 (7.2-8) 算法: 输入a,b,n,y0 h=(b-a)n,x0 = a for i = 1, i=n, i+ 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 y0 = y0 + h*(K1+2*K2+2*K3+K4)/6 输出x0,y0 16 2009, Henan Polytechnic University 163 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 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; 输入a,b,n,y0 h=(b-a)n,x0 = a for i = 1, i 时,说明步长h / 2仍然偏大,须将步长减半,继续计算; 当 / 2p时,说明已有 ,步长h / 2偏小,应取 步长h。 24 2009, Henan Polytechnic University 243 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 例2(P481) 用变步长的标准4阶R-K方法求初值问题: 的数值解,要求精度为
网站客服QQ:2055934822
金锄头文库版权所有
经营许可证:蜀ICP备13022795号 | 川公网安备 51140202000112号