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

常微分方程组的四阶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=zeros(1,N+1);

Y=zeros(1,N+1);

T=a:h:b;

X(1)=xa;

Y(1)=ya;

for j=1:N

f1=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(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.1200

0.0010 2.6637 0.0926

0.0020 3.7120 0.0774

0.0030 5.5033 0.0886

0.0040 4.9866 0.1193

0.0050 3.1930 0.1195

0.0060 2.7665 0.0951

0.0070 3.6543 0.0799

0.0080 5.2582 0.0884

0.0090 4.9942 0.1157

0.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.1000

0.0010 4.1883 0.1144

0.0020 3.2978 0.1072

0.0030 3.3468 0.0922

0.0040 4.2020 0.0876

0.0050 4.8807 0.0995

0.0060 4.2090 0.1126

0.0070 3.3874 0.1069

0.0080 3.4011 0.0934

0.0090 4.1568 0.0889

0.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:h:b;

t=zeros(1:N+1);

x(:,1)=x0;

for i=1:N+1

L1=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 9

3.0000 2.6637 3.7120 5.5033

4.9866 3.1930 2.7665 3.6543

5.2582 0.1200 0.0926 0.0774 0.0886 0.1193 0.1195 0.0951 0.0799 0.0884 Columns 10 through 12

4.9942 3.3541 2.8689

0.1157 0.1185 0.0970

结果分析：与1.1 的结果相同，当然这个方法还适合于一般方程组。