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

《常微分方程组的四阶rungekutta龙格库塔法matlab实现》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、常微分方程组的四阶Runge-Kutta方法1.问题：1.1若用普通方法-----仅适用于两个方程组成的方程组编程实现：创建M文件：functionR=rk4(f,g,a,b,xa,ya,N)%UNTITLED2Summaryofthisfunctiongoeshere%Detailedexplanationgoeshere%x'=f(t,x,y)y'=g(t,x,y)%N为迭代次数%h为步长%ya,xa为初值f=@(t,x,y)(2*x-0.02*x*y);g=@(t,x,y)(0.0002*x*y-0.8*y);h=

2、(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)+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)

3、+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=[T'X'Y'];end情况一：对于x0=3000，y0=120控制台中输入：>>rk4('f','g',0,10,300

4、0,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数据:情况二：对于x0=5000，y0=100命令行中输入：>>rk4('f','g'

5、,0,10,5000,100,10)运行结果：ans=1.0e+003*05.00000.10000.00104.18830.11440.00203.29780.10720.00303.34680.09220.00404.20200.08760.00504.88070.09950.00604.20900.11260.00703.38740.10690.00803.40110.09340.00904.15680.08890.01004.77530.0991数据：结论：无论取得初值是哪一组，捕食者与被捕食者的数量总是一个增长

6、另一个减少，并且是以T=5为周期交替增长或减少的。这说明了捕食者与被捕食者是对立的，一个增长则另一个必定减少，从而达到食物链相对稳定。1.2改进方法---------一般方法代码：创建M文件：functionx=rk4(f,x0,h,a,b,N)%x0初值%(a,b)为迭代区间%N迭代次数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/

7、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.7120

8、5.50334.98663.19302.76653.65435.25820.12000.09260.07740.08860.11930.11950.09510.07990.0884Columns10through124.99423.35412.86890.11570.11850.0970结果分析：与1.1的结果相同，