常微分方程组的四阶RungeKutta龙格库塔法matlab实现

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的数口.捕食者被捕食者模型断占.r1小卬满足J,却=Arf-工付巩方-t=Ur$一.网讣关于捕食者做捕食名模型的背景知识可参见同仁,李承治常微分方程载程?第一章南用举例设=C.O2.C=O.DQD2.Z=U8初蛤值取,rU=布加护=120h3I5C.1iK.jfD1001.1假设用普通方法-仅适用于两个方程组成的方程组编程实现：创立M文件：functionR=rk4(f,g,a,b,xa,ya,N)%UNTITLED2Summaryofthisfunctiongoeshere%Detailedexplanati

3、ongoeshere%x=f(t,x,y)y=g(t,x,y)%的迭代次数%的步长%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)运行结果：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)运行结果：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数据:时刻t017Am34缭子数5000418833297.S33-16.K4202狐史数-1001144107.292.287.65689104880-42093387434014156.84775.399.5112.6106.993.48&999.1结论：无论取得初值是哪一组,捕食者与被捕食者的数量总是一个增长另一个减少,并且是以T=5为周期交替增长或减少的.这说明了捕食者与被捕食者是对立的,一个增长那么另一个必定减少,从而到达食物链相对稳定.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)-+个函数向过代码：创立M文件：functionx=rk4(f,x0,h,a,b,N)%x0U值%(a,b)为迭代区间%港代次数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以第一组数据为例：对于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)运行结果：ans=1.0e+003*Columns1through93.00002.66373.71205.50334.98663.19302.76653.65435.25820