数值分析RungeKua方法.ppt

2009,Henan Polytechnic University,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,3 Runge-Kutta,方法,第五章 常微分方程数值解法,第三节,Runge-Kutta,方法,3,龙格,-,库塔法,/*,Runge-Kutta,Method*/,考察改进的欧拉法，可以将其改写为：,斜率,一定取,K,1,K,2,的,平均值,吗？,步长一定是一个,h,吗？,单步递推法的,基本思想,是从,(,x,i,y,i,),点出发，以,某一斜,率,沿直线达到,(,x,i,+1,y,i,+1,),点。欧拉法及其各种变形所,能达到的最高精度为,2,阶。,建立高精度的单步递推格式。,首先希望能确定系数,1,、,2,、,p,，使得到的算法格式有,2,阶精度，即在 的前提假设下，使得,Step 1,:,将,K,2,在,(,x,i,y,i,),点作,Taylor,展开,将改进欧拉法推广为：,),(,),(,1,2,1,2,2,1,1,1,p,hK,y,p,h,x,f,K,y,x,f,K,K,K,h,y,y,i,i,i,i,i,i,+,+,=,=,+,+,=,+,l,l,Step 2,:,将,K,2,代入第,1,式，得到,Step 3,:,将,y,i,+1,与,y,(,x,i,+1,),在,x,i,点的泰勒展开作比较,要求 ，则必须有：,这里有 个未知数，个方程。,3,2,存在,无穷多个解,。所有满足上式的格式统称为,2,阶龙格,-,库塔格式,。,注意到，就是改进的欧拉法。,Q:,为获得更高的精度，应该如何进一步推广？,其中,i,(,i,=1,m,),，,i,(,i,=2,m,),和,ij,(,i,=2,m,;,j,=1,i,1,),均为待定系数，确定这些系数的步骤与前面相似。,),.,(,.,.,),(,),(,),(,.,1,1,2,2,1,1,2,32,1,31,3,3,1,21,2,2,1,2,2,1,1,1,-,-,+,+,+,+,+,+,=,+,+,+,=,+,+,=,=,+,+,+,+,=,m,m m,mm,m,i,m,i,i,i,i,i,i,m,m,i,i,hK,hK,hK,y,h,x,f,K,hK,hK,y,h,x,f,K,hK,y,h,x,f,K,y,x,f,K,K,K,K,h,y,y,b,b,b,a,b,b,a,b,a,l,l,l,最常用为四级,4,阶,经典龙格,-,库塔法,/*Classical Runge-Kutta Method*/,：,注：,龙格,-,库塔法,的主要运算在于计算,K,i,的值，即计算,f,的值。,Butcher,于,1965,年给出了计算量与可达到的最高精度阶数的关系：,7,5,3,可达到的最高精度,6,4,2,每步须算,K,i,的个数,由于龙格,-,库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用,低阶算法,而将步长,h,取小,。,深入研究龙格,-,库塔法请看,此处,！,7.2 Runge,Kutta,方法,7.2.1,构造高阶单步法的直接方法,由,Taylor,公式：,当,h,充分小时，略去,Taylor,公式余项，并以,y,i,、,y,i,+1,分别代替,y,(,x,i,),、,y,(,x,i,+1,),，得到差分方程：,(7.2-1),其局部截断误差为：,当,h,充分小时，略去,Taylor,公式余项，并以,y,i,、,y,i,+1,分别代替,y,(,x,i,),、,y,(,x,i,+1,),，得到差分方程：,(7.2-1),其局部截断误差为：,即,(7.2-1),为,p,阶方式，上述方式称为,Taylor,方式。,注：利用,Tayler,公式构造，不实用，高阶导数,f,(,i,),不易计算。,7.2.2 Runge,Kutta,方法,1.,基本思想,因为,=,y,(,x,i,)+,hf,(,，,y,(,),=,y,(,x,i,)+,hK,其中,K,=,f,(,，,y,(,),称为,y,(,x,),在,x,i,，,x,i,+1,上的平均斜率。,若取,K,1,=,f,(,x,i,，,y,(,x,i,),Euler,公式,取,K,2,=,f,(,x,i,+1,，,y,(,x,i,+1,),向后Euler,公式 一阶精度,取,梯形公式 二阶精度,猜想：若能多预测几个点的斜率，再取其加权平均作为,K,，可望得到较高精度的数值解，从而避免求,f,的高阶导数。,2.R,K,公式,(7.2-4),其中,K,j,为,y,=,y,(,x,),在,x,i,+,a,j,h,(0,a,j,1),处的斜率预测值。,a,j,，,b,js,，,c,j,为特定常数。,3.,常数的确定,确定的原则是使精度尽可能高。,以二阶为例：,(7.2-5),(,希望,y,(,x,i,+1,),y,i,+1,=O(,h,p,),的阶数,p,尽可能高,),首先,：,另一方面,：,将,K,2,在,(,x,i,，,y,i,),处展开。,K,2,=,f,(,x,i,，,y,i,)+,a,2,hf,x,(,x,i,，,y,i,)+,b,21,hK,1,f,y,(,x,i,，,y,i,)+O(,h,2,).,代入,(7.2-5),得,：,y,i,+1,=,y,i,+,hc,1,f,(,x,i,，,y,i,)+,hc,2,f,(,x,i,，,y,i,),+,h,2,c,2,a,2,f,x,(,x,i,，,y,i,)+,b,21,K,1,f,y,(,x,i,，,y,i,)+O(,h,3,),=,y,i,+,h,(,c,1,+,c,2,),f,(,x,i,，,y,i,),+,c,2,a,2,h,2,f,x,(,x,i,，,y,i)+(,b,21,/,a,12,),f,(,x,i,，,y,i,),f,y,(,x,i,，,y,i,)+O(,h,3,),（希望）,2024/12/13 周五,13,希望,：,e,i,+1,=,y,(,x,i,+1,),y,i,+1,=O(,h,3,).,则应,：,特例：,a,2,=1,c,1,=,c,2,=1/2,，,b,21,=1,，得,2,阶,R-K,公式：,改进欧拉公式。,希望,：,e,i,+1,=,y,(,x,i,+1,),y,i,+1,=O(,h,3,).,则应,：,特例：,c,1,=0,c,2,=1,，,a,2,=1/2,，,b,21,=1/2,，得：,(7.2-7),称为中点公式。,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,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,仍然偏大，须将步长减半，继续计算；,当,/2,p,时，说明已有 ，步长,h,/2,偏小，应取步长,h,。,例,2(P481),用变步长的标准,4,阶,R-K,方法求初值问题：,的数值解，要求精度为,2024/12/13 周五,26,