数值分析Runge-Kutta方法.ppt

《数值分析Runge-Kutta方法.ppt》由会员分享，可在线阅读，更多相关《数值分析Runge-Kutta方法.ppt（24页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、1第三节第三节第三节第三节 Runge-KuttaRunge-KuttaRunge-KuttaRunge-Kutta方法方法方法方法2 2009,Henan Polytechnic University2 3 3 3 3 Runge-KuttaRunge-KuttaRunge-KuttaRunge-Kutta方法方法方法方法第五章第五章第五章第五章 常微分方程数值解法常微分方程数值解法常微分方程数值解法常微分方程数值解法3 龙格龙格-库塔法库塔法/*Runge-Kutta Method*/考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：斜率斜率一定取一定取K1 K2 的的平

2、均值平均值吗？吗？步长一定是一个步长一定是一个h 吗吗？单步递推法的单步递推法的基本思想基本思想是从是从(xi,yi)点出发，以点出发，以某一斜某一斜率率沿直线达到沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所点。欧拉法及其各种变形所能达到的最高精度为能达到的最高精度为2阶阶。建立高精度的单步递推格式。建立高精度的单步递推格式。3 2009,Henan Polytechnic University3 3 3 3 3 Runge-KuttaRunge-KuttaRunge-KuttaRunge-Kutta方法方法方法方法第五章第五章第五章第五章 常微分方程数值解法常微分方程数值解法常微

3、分方程数值解法常微分方程数值解法首先希望能确定系数首先希望能确定系数 1、2、p，使得到的算法格式有使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 Step 1:将将 K2 在在(xi,yi)点作点作 Taylor 展开展开将将改进欧拉法推广为：改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+Step 2:将将 K2 代入第代入第1式，得到式，得到4 2009,Henan Polytechnic University4 3 3 3 3 Runge-KuttaRunge-KuttaRunge-KuttaR

4、unge-Kutta方法方法方法方法第五章第五章第五章第五章 常微分方程数值解法常微分方程数值解法常微分方程数值解法常微分方程数值解法Step 3:将将 yi+1 与与 y(xi+1)在在 xi 点点的的泰勒泰勒展开作比较展开作比较要求要求 ，则必须有：，则必须有：这里有这里有 个未知个未知数，数，个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格-库库塔格式塔格式。注意到，注意到，就是改进的欧拉法。就是改进的欧拉法。Q:为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？5 2009,Henan Poly

5、technic University5 3 3 3 3 Runge-KuttaRunge-KuttaRunge-KuttaRunge-Kutta方法方法方法方法第五章第五章第五章第五章 常微分方程数值解法常微分方程数值解法常微分方程数值解法常微分方程数值解法其中其中 i (i=1,m)，i (i=2,m)和和 ij(i=2,m;j=1,i 1)均为待均为待定系数，确定这些系数定系数，确定这些系数的步骤与前面相似。的步骤与前面相似。).,(.),(),(),(.1122112321313312122122111 +=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhK

6、yhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法/*Classical Runge-Kutta Method*/：6 2009,Henan Polytechnic University6 3 3 3 3 Runge-KuttaRunge-KuttaRunge-KuttaRunge-Kutta方法方法方法方法第五章第五章第五章第五章 常微分方程数值解法常微分方程数值解法常微分方程数值解法常微分方程数值解法注：注：龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于

7、1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度642每步须算Ki 的个数 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。深入研究龙格深入研究龙格-库塔法请看库塔法请看此处此处！7 2009,Henan Polytechnic University7 3 3 3 3 Runge-KuttaRunge-KuttaRunge-KuttaR

8、unge-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 University8 3 3 3 3 Runge-KuttaRunge-KuttaRunge-KuttaRunge-Kutta方法方法方法方法第五章第五章第五章第五章 常微分方程数值解

9、法常微分方程数值解法常微分方程数值解法常微分方程数值解法当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 University9 3 3 3 3 Runge-KuttaRunge-KuttaRunge-KuttaRunge-Kutta方法方法方法方法第五章第五章第五章第五章 常微分方程数值解法常微分方程数值解法常微分方程数值解法常微分方

10、程数值解法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 University10 3 3 3 3 Runge-KuttaRunge-KuttaRunge-KuttaRunge-Kutta方法

11、方法方法方法第五章第五章第五章第五章 常微分方程数值解法常微分方程数值解法常微分方程数值解法常微分方程数值解法2.RK公式 (7.2-4)其中Kj为y=y(x)在xi+ajh(0 aj 1)处的斜率预测值。aj，bjs，cj为特定常数。11 2009,Henan Polytechnic University11 3 3 3 3 Runge-KuttaRunge-KuttaRunge-KuttaRunge-Kutta方法方法方法方法第五章第五章第五章第五章 常微分方程数值解法常微分方程数值解法常微分方程数值解法常微分方程数值解法3.常数的确定确定的原则是使精度尽可能高。以二阶为例：(7.2-5)

12、(希望y(xi+1)yi+1=O(hp)的阶数p尽可能高)首先：12 2009,Henan Polytechnic University12 3 3 3 3 Runge-KuttaRunge-KuttaRunge-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)+h2c2a

13、2 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 University13 3 3 3 3 Runge-KuttaRunge-KuttaRunge-KuttaRunge-Kutta方法方法方法方法第五章第五章第五章第五章 常微分方程数值解法常微分方程数值解法常微分方程数值解法常微分方程数值解法希望：ei+1=y(xi+1)yi+1=O(h3).则应：特例：a2=1 c1=c2

14、=1/2，b21=1，得2阶R-K公式：改进欧拉公式。14 2009,Henan Polytechnic University14 3 3 3 3 Runge-KuttaRunge-KuttaRunge-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 University15 3 3 3 3 Runge

15、-KuttaRunge-KuttaRunge-KuttaRunge-Kutta方法方法方法方法第五章第五章第五章第五章 常微分方程数值解法常微分方程数值解法常微分方程数值解法常微分方程数值解法4.最常用的R-K公式 标准4阶R-K公式 (7.2-8)算法：输入输入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/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，y016 2009,Henan Pol

16、ytechnic University16 3 3 3 3 Runge-KuttaRunge-KuttaRunge-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+K

17、4)/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 时，说明步长h/2仍然偏大，须将步长减半，继续计算；当 /2p时，说明已有 ，步长h/2偏小，应取步长h。24 2009,Henan Polytechnic University24 3 3 3 3 Runge-KuttaRunge-KuttaRunge-KuttaRunge-Kutta方法方法方法方法第五章第五章第五章第五章 常微分方程数值解法常微分方程数值解法常微分方程数值解法常微分方程数值解法例2(P481)用变步长的标准4阶R-K方法求初值问题：的数值解，要求精度为