数值分析RungeKutta方法.pptx

3 龙格龙格-库塔法库塔法/*Runge-Kutta Method*/考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：斜率斜率一定取一定取K1 K2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗吗？单步递推法的单步递推法的基本思想基本思想是从是从(xi,yi)点出发，以点出发，以某一斜某一斜率率沿直线达到沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所点。欧拉法及其各种变形所能达到的最高精度为能达到的最高精度为2阶阶。建立高精度的单步递推格式。建立高精度的单步递推格式。第1页/共24页首先希望能确定系数首先希望能确定系数 1、2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 Step 1:将将 K2 在在(xi,yi)点作点作 Taylor 展开展开将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+Step 2:将将 K2 代入第代入第1式，得到式，得到第2页/共24页Step 3:将将 yi+1 与与 y(xi+1)在在 xi 点的点的泰勒泰勒展开作比较展开作比较要求要求 ，则必须有：，则必须有：这里有这里有 个未知个未知数，数，个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格-库库塔格式塔格式。注意到，注意到，就是改进的欧拉法。就是改进的欧拉法。Q:为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？第3页/共24页其中其中 i (i=1,m)，i (i=2,m)和和 ij(i=2,m;j=1,i 1)均为待均为待定系数，确定这些系数定系数，确定这些系数的步骤与前面相似。的步骤与前面相似。).,(.),(),(),(.1122112321313312122122111 +=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法/*Classical Runge-Kutta Method*/：第4页/共24页注：注：龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度642每步须算Ki 的个数 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。深入研究龙格深入研究龙格-库塔法请看库塔法请看此处此处！第5页/共24页7.2 RungeKutta方法 7.2.1 构造高阶单步法的直接方法由Taylor公式：当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程：(7.2-1)其局部截断误差为：第6页/共24页当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程：(7.2-1)其局部截断误差为：即(7.2-1)为p阶方式，上述方式称为Taylor方式。注：利用Tayler公式构造，不实用，高阶导数f(i)不易计算。第7页/共24页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 的高阶导数。第8页/共24页2.RK公式 (7.2-4)其中Kj为y=y(x)在xi+ajh(0 aj 1)处的斜率预测值。aj，bjs，cj为特定常数。第9页/共24页3.常数的确定确定的原则是使精度尽可能高。以二阶为例：(7.2-5)(希望y(xi+1)yi+1=O(hp)的阶数p尽可能高)首先：第10页/共24页另一方面：将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)（希望）第11页/共24页希望：ei+1=y(xi+1)yi+1=O(h3).则应：特例：a2=1 c1=c2=1/2，b21=1，得2阶R-K公式：改进欧拉公式。第12页/共24页希望：ei+1=y(xi+1)yi+1=O(h3).则应：特例：c1=0 c2=1，a2=1/2，b21=1/2，得：(7.2-7)称为中点公式。第13页/共24页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，y0第14页/共24页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=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。第22页/共24页例2(P481)用变步长的标准4阶R-K方法求初值问题：的数值解，要求精度为第23页/共24页感谢您的观看！第24页/共24页