資源描述：
《數(shù)值分析教教案15》由會員上傳分享，免費在線閱讀，更多相關內容在工程資料-天天文庫。
1、3.5龍格-庫塔法龍格-庫塔方法的基本原理是：用方程p+i=X+〃g（W，X）[y^o）=yo中函數(shù)g（九y）在前一節(jié)點心上取值g（兀，％）的線性組合構造一個表示的近似值公式，從而避免了求X+i時用&（兀刃的高階導數(shù)。該方法有多種推導方法，這里用數(shù)值積分法推導。為此，先將微分方程略加變形得出：dy=g（x,y）dxf對該式兩邊在相鄰兩節(jié)點旺和￡+1之間求積分，移項得出：sg(x,y)dxX：(3-15)采用不同的方法計算式（3-15）中定積分，便可得出數(shù)值積分的不同近似結果。如果用矩形求積公式計算，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)右邊的函數(shù)是丘和K2的線性組合，而人和K?是把兀門X和兀+i=xi+h的值代入函數(shù)g(x,y)得出的。這樣計算X+i時不再像用Taylor公式那樣求g(“*J的導數(shù)。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實現(xiàn)三階龍格-庫塔法計算公式（還有休恩法）為：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、微分方程的一般表達式的右端函數(shù)%h：積分步長%a：自變量取值下限%b：自變量取值上限%yO：函數(shù)初值%varvec：常微分方程的變量組%subs(S)表示：用數(shù)值替代所有的符號變量。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函數(shù)運行時需要調用下列函數(shù)：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、格-庫塔法求下面常微分方程的數(shù)值解。0<%<1包=2兀-3y+2
