计算方法5.常微分方程的数值解法概要.ppt

《计算方法5.常微分方程的数值解法概要.ppt》由会员上传分享，免费在线阅读，更多相关内容在PPT专区-天天文库。

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

2、hxfKyxfKKKhyyiiiiii++==++=+llStep2:将K2代入第1式，得到10/1/20212Step3:将yi+1与y(xi+1)在xi点的泰勒展开作比较要求，则必须有：这里有个未知数，个方程。32存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。注意到，就是改进的欧拉法。Q:为获得更高的精度，应该如何进一步推广？10/1/20213构造高阶单步法的直接方法由Taylor公式：当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程：其局部截断

3、误差为：差分方程为p阶方式，上述方式称为Taylor方式。注：利用Tayler公式构造，不实用，高阶导数f(i)不易计算。10/1/20214RungeKutta方法1.基本思想因为其中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/1/202152.R

4、K公式其中Kj为y=y(x)在xi+ajh(0aj1)处的斜率预测值。aj，bjs，cj为特定常数。10/1/202163.常数的确定确定的原则是使精度尽可能高。以二阶为例：(希望y(xi+1)–yi+1=O(hp)的阶数p尽可能高)首先：10/1/20217另一方面：将K2在(xi，yi)处展开。K2=f(xi,yi)+a2hfx'(xi,yi)+b21hK1fy'(xi,yi)+O(h2)可得：yi+1=yi+hc1f(xi,yi)+hc2f(xi,yi)+h2c2[a2fx'(xi,yi)+b21K1fy'(x

5、i,yi)]+O(h3)=yi+h(c1+c2)f(xi,yi)+c2a2h2[fx'(xi,yi)+(b21/a2)f(xi,yi)fy'(xi,yi)]+O(h3)（希望）10/1/20218希望：ei+1=y(xi+1)–yi+1=O(h3).则应：特例：a2=1c1=c2=1/2，b21=1，得2阶R-K公式：改进欧拉公式。10/1/20219希望：ei+1=y(xi+1)–yi+1=O(h3).则应：特例：c1=0c2=1，a2=1/2，b21=1/2，得：称为中点公式。10/1/2021104.最常用的

6、R-K公式——标准4阶R-K公式算法：输入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+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，y010/1/202111Matlab代码：functionRunge_Kutta4(a,b,h,y0)n=(b-a)/h;x0=a;fori=1:nK1=f(x0,y0)K2=f(x0+h

7、/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;输入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+h/2,y0+h*K2/2)K4=f(x0+h,y0+h*K3)x0=x0+hy0=y0+h*(K1+2*K2+2*

8、K3+K4)/6输出x0，y010/1/202112例：用标准4阶R-K公式求：的数值解。取h=0.2，并与标准解y=2ex–x–1比较。解：因为f(x，y)=x+y，从由得：10/1/202113xi+1K1K2K3K4yi+1y(xi+1)0.211.21.221.4441.2428001.2428