资源描述：
《数值分析教教案15》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。
1、3.5龙格-库塔法龙格-库塔方法的基本原理是：用方程p+i=X+〃g（W，X）[y^o）=yo中函数g（九y）在前一节点心上取值g（兀，％）的线性组合构造一个表示的近似值公式，从而避免了求X+i时用&（兀刃的高阶导数。该方法有多种推导方法，这里用数值积分法推导。为此，先将微分方程略加变形得出：dy=g（x,y）dxf对该式两边在相邻两节点旺和£+1之间求积分，移项得出：sg(x,y)dxX：(3-15)采用不同的方法计算式（3-15）中定积分，便可得出数值积分的不同近似结果。如果用矩形求积公式计算，JF+1g（兀刃况¥=（兀+i—兀）g（兀，X）=/zg（X,」J代入式（3-
2、15）就可得出X+l=X+人？（兀厂必）,这个结果和由Taylor公式得出的Euler公式结果完全一样。3.5.1二阶龙格-库塔公式若用梯形求积公式计算式（3-15）中的积分，则有:F+11g(x,y)dx^-h[g(xi.yJxi2J+g(Gl,x+l)]令K]=g(%y)上式中g(兀+iJ+J里的无+i=兀+人为+1用欧拉公式代换，则可得出X+l=X+/2g(£』J=X+/2K]，把勺+i和X+1代入上式则有:1g(兀,y)dx=t〃[K]+g(xi+h,yi+hK})]令：K?=gg+hdi+flKJ,则得出式(3-15)的一个近似结果:1yM=y^-h(K}+K2)(
3、3-16)这就是二阶龙-格库塔公式，它的局部截断差为°(力彳)。式(3-16)右边的函数是丘和K2的线性组合，而人和K?是把兀门X和兀+i=xi+h的值代入函数g(x,y)得出的。这样计算X+i时不再像用Taylor公式那样求g(“*J的导数。3.5.2三阶龙格-库塔公式若用抛物线(Simpson)求积公式计算式(3-15)中的积分，则有:3+i1Ig(ly)dx--灿g(“y)+4g(和⑵X+1/2)+,X+1)1Jxi6=x{+h.h=XM_兀，兀+1/2=兀+力/2,兀+]X+1/2和兀+1都是未知的，可以用欧拉格式估算X+I/2二X+1处（兀・J）和X+l=兀+力g（
4、X」）o类似二阶龙格-库塔公式的推导，并令:Ki二”/hh“、K2=g(乞+亍％+空(),K3=g(xt+h.yt-hK、+2hK2)把它们代入式（3-15）则得出三阶龙格-库塔公式，它的局部截断误差为曲）。h九r+护+%+()(3-17)3.5.3三阶龙格-库塔公式的MATLAB实现三阶龙格-库塔法计算公式（还有休恩法）为：Ki二gCwJ“/hh“、K2=g(无+〒”+qKJK3=g{xL+h.yi-hK、+2hK2)h—+泸+4D三阶龙格■库塔公式的MATLAB程序代码如下所示:functiony=DELGKT3_kuta（f,h,a,b,yO,varvec）%f:一阶常
5、微分方程的一般表达式的右端函数%h：积分步长%a：自变量取值下限%b：自变量取值上限%yO：函数初值%varvec：常微分方程的变量组%subs(S)表示：用数值替代所有的符号变量。formatlong;N=(b-a)/h;y=zeros(N+l,l);y(l)=yO;x=a:h:b;var二findsym(f);fori=2:N+lKI=Funval(f,varvec,[x(i-1)y(i-l)]);K2=Funval(f,varvec,[x(i-1)+h/2y(i・l)+Kl*h/2]);K3=Funval(f,varvec,[x(i-1)+hy(i-l)・h*Kl+K2
6、*2*h]);y(i)二y(i-l)+h*(Kl+4*K2+K3)/6;endformatshort;DELGKT3_kuta函数运行时需要调用下列函数：functionfv=Funval(f,varvec,varval)var=findsym(f);讦length(var)<4ifvar(1)==varvec(l)fv=subs(f,varvec(1),varval(1));elsefv=subs(f,varvec(2),varval(2));endelsefv二subs(f,varvec,varval);end【例3-9】三阶龙格-库塔求解一阶常微分方程应用实例。用三阶龙
7、格-库塔法求下面常微分方程的数值解。0<%<1包=2兀-3y+2
