《常微分方程组的四阶RungeKutta龙格库塔法matlab实现》由会员分享，可在线阅读，更多相关《常微分方程组的四阶RungeKutta龙格库塔法matlab实现（5页珍藏版）》请在金锄头文库上搜索。
1、 常微分方程组的四阶Runge-Kutta方法1.问题：1.1若用普通方法-仅适用于两个方程组成的方程组编程实现：创建M 文件：function R = rk4(f,g,a,b,xa,ya,N)%UNTITLED2 Summary of this function goes here% Detailed explanation goes here%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=(b-a)/N;T=zeros(1,N+1);X=zero
2、s(1,N+1);Y=zeros(1,N+1);T=a:h:b;X(1)=xa;Y(1)=ya;for j=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)+h*g2/2);g3=feval(g,T(j)+h/2,X(j)+h/2*f2,Y(j)+h/2*g2);f4=feval(f,T(
3、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,3000,120,10)运行结果：ans =1.0e+003 *0 3.0000 0.12000.0010 2.6637 0.09260.0020 3.7120 0.07740.0030 5.5033 0.08860.0040 4
4、.9866 0.11930.0050 3.1930 0.11950.0060 2.7665 0.09510.0070 3.6543 0.07990.0080 5.2582 0.08840.0090 4.9942 0.11570.0100 3.3541 0.1185数据:情况二：对于x0=5000，y0=100命令行中输入： rk4(f,g,0,10,5000,100,10)运行结果：ans =1.0e+003 *0 5.0000 0.10000.0010 4.1883 0.11440.0020 3.2978 0.10720.0030 3.3468 0.09220.0040 4.2020 0.0
5、8760.0050 4.8807 0.09950.0060 4.2090 0.11260.0070 3.3874 0.10690.0080 3.4011 0.09340.0090 4.1568 0.08890.0100 4.7753 0.0991数据：结论：无论取得初值是哪一组，捕食者与被捕食者的数量总是一个增长另一个减少，并且是以T=5 为周期交替增长或减少的。这说明了捕食者与被捕食者是对立的，一个增长则另一个必定减少，从而达到食物链相对稳定。1.2改进方法-一般方法代码：创建M 文件：function x=rk4(f,x0,h,a,b,N)%x0初值%(a,b)为迭代区间%N迭代次数t=a
6、:h:b;t=zeros(1:N+1);x(:,1)=x0;for i=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 *Columns 1 through 93.0000 2.6637 3.7120 5.5033 4.9866 3.1930 2.7665 3.6543 5.25820.1200 0.0926 0.0774 0.0886 0.1193 0.1195 0.0951 0.0799 0.0884Columns 10 through 124.9942 3.3541 2.86890.1157 0.1185 0.0970结果分析：与1.1 的结果相同，当然这个方法还适合于一般方程组。
