数值分析Runge-Kutta方法PPT课件

第三节Runge-Kutta方法,3龙格-库塔法/*Runge-KuttaMethod*/,斜率一定取K1K2的平均值吗？,步长一定是一个h吗？,单步递推法的基本思想是从(xi,yi)点出发，以某一斜率沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。,建立高精度的单步递推格式。,首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得,Step1:将K2在(xi,yi)点作Taylor展开,Step2:将K2代入第1式，得到,Step3:将yi+1与y(xi+1)在xi点的泰勒展开作比较,要求，则必须有：,这里有个未知数，个方程。,3,2,存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。,注意到，就是改进的欧拉法。,Q:为获得更高的精度，应该如何进一步推广？,其中i(i=1,m)，i(i=2,m)和ij(i=2,m;j=1,i1)均为待定系数，确定这些系数的步骤与前面相似。,最常用为四级4阶经典龙格-库塔法/*ClassicalRunge-KuttaMethod*/：,由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h取小。,深入研究龙格-库塔法请看此处！,7.2RungeKutta方法7.2.1构造高阶单步法的直接方法由Taylor公式：当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程：(7.2-1)其局部截断误差为：,当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程：(7.2-1)其局部截断误差为：即(7.2-1)为p阶方式，上述方式称为Taylor方式。注：利用Tayler公式构造，不实用，高阶导数f(i)不易计算。,7.2.2RungeKutta方法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的高阶导数。,2.RK公式(7.2-4)其中Kj为y=y(x)在xi+ajh(0aj1)处的斜率预测值。aj，bjs，cj为特定常数。,3.常数的确定确定的原则是使精度尽可能高。以二阶为例：(7.2-5)(希望y(xi+1)yi+1=O(hp)的阶数p尽可能高)首先：,另一方面：将K2在(xi，yi)处展开。K2=f(xi，yi)+a2hfx(xi，yi)+b21hK1fy(xi，yi)+O(h2).代入(7.2-5)得：yi+1=yi+hc1f(xi，yi)+hc2f(xi，yi)+h2c2a2fx(xi，yi)+b21K1fy(xi，yi)+O(h3)=yi+h(c1+c2)f(xi，yi)+c2a2h2fx(xi，yi)+(b21/a12)f(xi，yi)fy(xi，yi)+O(h3)（希望）,希望：ei+1=y(xi+1)yi+1=O(h3).则应：特例：a2=1c1=c2=1/2，b21=1，得2阶R-K公式：改进欧拉公式。,希望：ei+1=y(xi+1)yi+1=O(h3).则应：特例：c1=0c2=1，a2=1/2，b21=1/2，得：(7.2-7)称为中点公式。,4.最常用的R-K公式标准4阶R-K公式(7.2-8)算法：,Matlab代码：functionRunge_Kutta4(a,b,h,y0)n=(b-a)/h;x0=a;fori=1:nK1=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+hy1=y0+h*(K1+2*K2+2*K3+K4)/6;y0=y1end;end;functionf=f(x,y)f=x+y;end;,例1：（P478）用标准4阶R-K公式求：的数值解。取h=0.2，并与标准解y=2exx1比较。解：因为f(x，y)=x+y，从而由(7.2-8)得：,例1：（P478）用标准4阶R-K公式求：的数值解。取h=0.2，并与标准解y=2exx1比较。解：,注：步长h的选择(P479)使用数值解法求解初值问题，选择步长h是一个重要问题。从每一步看，h小局部截断误差小，整体截断误差也就小；但从整个区间看，h小则节点多，这不仅使计算工作量增大，而且也使舍入误差的累积严重，导致数值不稳定。所以步长h的大小要适当。在实际计算中，常常采用事后估计误差及自动选择步长的办法来保证节点xi处的数值解yi(i=1，2，N)满足所要求的精度。,设所用的一步法是p阶方法，其整体截断误差有渐近公式y(xi+1)yi+1=Mhp+O(hp+1)其中M是与h无关的常数。考虑用Richardson外推法提高数值解的精度。事实上，从节点xi出发，先以h为步长经一步计算求出数值解，则整体截断误差如上式所示为(7.2-9)然后以h/2为步长，从节点xi出发经过两步计算求出数值解，则整体截断误差为(7.2-10)用2p乘(7.2-10)式减去(7.2-9)式，得,用2p乘(7.2-10)式减去(7.2-9)式，得故有：(7.2-11)(7.2-12)或：(7.2-13),注：(1)(7.2-11)表明若取则整体截断误差为hp+1，即精度提高一阶。(7.2-12)表明数值解的整体截断误差近似为(7.2-13)表明数值解的整体截断误