数值分析Runge-Kutta方法.ppt

《数值分析Runge-Kutta方法.ppt》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、第三节Runge-Kutta方法§3龙格-库塔法/*Runge-KuttaMethod*/考察改进的欧拉法，可以将其改写为：斜率一定取K1K2的平均值吗？步长一定是一个h吗？单步递推法的基本思想是从(xi,yi)点出发，以某一斜率沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。建立高精度的单步递推格式。首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得Step1:将K2在(xi,yi)点作Taylor展开将改进欧拉法推广为：),(),(][12122111phKyphxfKyx

2、fKKKhyyiiiiii++==++=+llStep2:将K2代入第1式，得到Step3:将yi+1与y(xi+1)在xi点的泰勒展开作比较要求，则必须有：这里有个未知数，个方程。32存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。注意到，就是改进的欧拉法。Q:为获得更高的精度，应该如何进一步推广？其中i(i=1,…,m)，i(i=2,…,m)和ij(i=2,…,m;j=1,…,i1)均为待定系数，确定这些系数的步骤与前面相似。)...,(......),(),(),(]...[1122112321313312122122

3、111--++++++=+++=++==++++=mmmmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyybbbabbaballl最常用为四级4阶经典龙格-库塔法/*ClassicalRunge-KuttaMethod*/：注：龙格-库塔法的主要运算在于计算Ki的值，即计算f的值。Butcher于1965年给出了计算量与可达到的最高精度阶数的关系：753可达到的最高精度642每步须算Ki的个数由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好

4、采用低阶算法而将步长h取小。深入研究龙格-库塔法请看此处！7.2RungeKutta方法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.

5、2.2RungeKutta方法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.RK公式(7.2-4)其中Kj为y=y(x)在xi+ajh(0aj1)处的斜率预测值。aj，bjs，cj为特定常数。3.

6、常数的确定确定的原则是使精度尽可能高。以二阶为例：(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)+h2c2[a2fx'(xi，yi)+b21K1fy'(xi，yi)]+O(h3)=yi+h(c1+c2)f(xi，yi)+c2a2h2[fx'(xi，yi)+(b21/a12)f(xi，yi)fy

7、'(xi，yi)]+O(h3)（希望）希望：ei+1=y(xi+1)–yi+1=O(h3).则应：特例：a2=1c1=c2=1/2，b21=1，得2阶R-K公式：改进欧拉公式。希望：ei+1=y(xi+1)–yi+1=O(h3).则应：特例：c1=0c2=1，a2=1/2，b21=1/2，得：(7.2-7)称为中点公式。4.最常用的R-K公式——标准4阶R-K公式(7.2-8)算法：输入a，b，n，y0h=(b-a)n，x0=afori=1,i<=n,i++K1=f(x0,y0)K2=f(x0+h/2,y0+h*K1/2)K3=f(x0

8、+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，y0M