常微分方程組的四階RungeKutta龍格庫塔法matlab實現(xiàn)

1、常微分方程組的四階Runge-Kutta方法1.問題:卜加以兩個方程組成的常微分方程組為例說明之,考慮初值問時-7T=f,%價*=|虱工2)tviih丁(hj=y(to)=蜘其四冊RungmKutUi方法可表小如卜:+2/2+2/i+/*Hk=強十gi+2例+%其中f/i=flk：工心叫c：ffi二上以,jA-thhhhhhh=/Pt+/HR+x/i馱+彳的=觀士+于也十./,除+n；BRA玄、tnth73=,%十歹以+4f?、期h+斗鈍,ff3=磯+不,1+不為,?卜一可禽卜Jul1f4=/ijfe+九工4hh羿+h驕,S4=!tfe+h,Xk+Vs,tfk+M器.匹1r推自名一班輔食型Ip

2、ffidator-pppywklel電口,n分別表示兔子被捕食者與狐貍捕食者在時刻1的數(shù)口.捕食者被捕食者模型斷占.r1小卬滿足J,卻=Arf-工付鞏方-t=Ur$一.網(wǎng)訃關(guān)于捕食者做捕食名模型的背景知識可參見同仁,李承治常微分方程載程?第一章南用舉例設(shè)=C.O2.C=O.DQD2.Z=U8初蛤值取,rU=布加護(hù)=120h3I5C.1iK.jfD1001.1假設(shè)用普通方法-僅適用于兩個方程組成的方程組編程實現(xiàn)：創(chuàng)立M文件：functionR=rk4(f,g,a,b,xa,ya,N)%UNTITLED2Summaryofthisfunctiongoeshere%Detailedexplanati

3、ongoeshere%x=f(t,x,y)y=g(t,x,y)%的迭代次數(shù)%的步長%ya,xa為初值f=(t,x,y)(2*x-0.02*x*y);g=(t,x,y)(0.0002*x*y-0.8*y);h=(b-a)/N;T=zeros(1,N+1);X=zeros(1,N+1);Y=zeros(1,N+1);T=a:h:b;X(1)=xa;Y(1)=ya;forj=1:Nf1=feval(f,T(j),X(j),Y(j);g1=feval(g,T(j),X(j),Y(j);f2=feval(f,T(j)+h/2,X(j)+h/2*f1,Y(j)+g1/2);g2=feval(g,T(j)+

4、h/2,X(j)+h/2*f1,Y(j)+h/2*g1);f3=feval(f,T(j)+h/2,X(j)+h/2*f2,Y(j)+h*g2/2);g3=feval(g,T(j)+h/2,X(j)+h/2*f2,Y(j)+h/2*g2);f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3);g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3);X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6;Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6;R=TXY;end情況一：對于x0=3000,y0=120限制臺中輸入

5、：rk4(f,g,0,10,3000,120,10)運行結(jié)果：ans=1.0e+003*03.00000.12000.00102.66370.09260.00203.71200.07740.00305.50330.08860.00404.98660.11930.00503.19300.11950.00602.76650.09510.00703.65430.07990.00805.25820.08840.00904.99420.11570.01003.35410.1185時刻t01Al34兔子教2663.7371205503349866孤獨教1209工677.488.6119356789103.

6、1930256653.65435.25824994233541119.595.179.9帆41157118.5情況二：對于x0=5000,y0=100命令行中輸入：rk4(f,g,0,10,5000,100,10)運行結(jié)果：ans=1.0e+003*05.00000.10000.00104.18830.11440.00203.29780.10720.00303.34680.09220.00404.20210.08760.00504.88070.09950.00604.20900.11260.00703.38740.10690.00803.40110.09340.00904.15680.0889

7、0.01004.77530.0991數(shù)據(jù):時刻t017Am34繚子數(shù)5000418833297.S33-16.K4202狐史數(shù)-1001144107.292.287.65689104880-42093387434014156.84775.399.5112.6106.993.48&999.1結(jié)論：無論取得初值是哪一組,捕食者與被捕食者的數(shù)量總是一個增長另一個減少,并且是以T=5為周期交替增長或減少的.這說明了捕食者與被捕食者是對立的,一個增長那么另一個必定減少,從而到達(dá)食物鏈相對穩(wěn)定.1.2改良方法-般方法試給舟般的階律俄分方程綱r譏、-J=Jt.l,-“口.)dtdx.2i-,叫,二2,工m)

8、with-m、dX即d,二/(小DX=(xLx2xii).X(tO)=XO：f也是M如)=邛*2,：如)_埠91-941l工作(0)-+個函數(shù)向過代碼：創(chuàng)立M文件：functionx=rk4(f,x0,h,a,b,N)%x0U值%(a,b)為迭代區(qū)間%港代次數(shù)t=a:h:b;t=zeros(1:N+1);x(:,1)=x0;fori=1:N+1L1=f(t(i),x(:,i);L2=f(t(i)+h/2,x(:,i)+(h/2)*L1);L3=f(t(i)+h/2,x(:,i)+(h/2)*L2);L4=f(t(i)+h,x(:,i)+h*L3);x(:,i+1)=x(:,i)+(h/6)*(L1+2*L2+2*L3+L4);end以第一組數(shù)據(jù)為例：對于x0=3000,y0=120限制臺中輸入：f=(t,x)2*x(1)-0.02*x(1)*x(2),0.0002*x(1)*x(2)-0.8*x(2);x0=3000,120;N=10;a=0;b=10;rk4(f,x0,a,b,10)運行結(jié)果：ans=1.0e+003*Columns1through93.00002.66373.71205.50334.98663.19302.76653.65435.25820