Fourth-order Runge–Kutta method (RK4) collapses after a few iterations

Question

I'm trying to solve:

x' = 60*x - 0.2*x*y;
y' = 0.01*x*y - 100* y;

using the fourth-order Runge-Kutta algorithm.

Starting points: x(0) = 8000, y(0) = 300 range: [0,15]

Here's the complete function:

function [xx yy time r] = rk4_m(x,y,step)
A = 0;
B = 15;

h = step;
iteration=0;
t = tic;

xh2 = x;
yh2 = y;


rr = zeros(floor(15/step)-1,1);
xx = zeros(floor(15/step)-1,1);
yy = zeros(floor(15/step)-1,1);
AA = zeros(1, floor(15/step)-1);

while( A < B)


    A = A+h;
    iteration = iteration + 1;

    xx(iteration) = x;
    yy(iteration) = y;
    AA(iteration) = A;
    [x y] = rkstep(x,y,h);


    for h2=0:1
        [xh2 yh2] = rkstep(xh2,yh2,h/2);
    end
    r(iteration)=abs(y-yh2);



end
time = toc(t);

xlabel('Range');
ylabel('Value');     
hold on
plot(AA,xx,'b');
plot(AA,yy,'g');
plot(AA,r,'r');
fprintf('Solution:\n');
fprintf('x: %f\n', x);
fprintf('y: %f\n', y);
fprintf('A: %f\n', A);
fprintf('Time: %f\n', time);

end

function [xnext, ynext] = rkstep(xcur, ycur, h)
    kx1 = f_prim_x(xcur,ycur);
    ky1 = f_prim_y(xcur,ycur);

    kx2 = f_prim_x(xcur+0.5*h,ycur+0.5*h*kx1);
    kx3 = f_prim_x(xcur+0.5*h,ycur+0.5*h*kx2);
    kx4 = f_prim_x(xcur+h,ycur+h*kx3);

    ky2 = f_prim_y(xcur+0.5*h*ky1,ycur+0.5*h);
    ky3 = f_prim_y(xcur+0.5*h*ky2,ycur+0.5*h);
    ky4 = f_prim_y(xcur+h*ky2,ycur+h);

    xnext = xcur + (1/6)*h*(kx1 + 2*kx2 + 2*kx3 + kx4);
    ynext = ycur + (1/6)*h*(ky1 + 2*ky2 + 2*ky3 + ky4);
end

function [fx] = f_prim_x(x,y)
fx = 60*x - 0.2*x*y;
end

function [fy] = f_prim_y(x,y)
fy = 0.01*x*y - 100*y;
end

And I'm running it by executing: [xx yy time] = rk4_m(8000,300,10)

The problem is that everything collapses after 2-3 iterations returning useless results. What am I doing wrong? Or is just this method not appropriate for this kind equation?

The semicolons are intentionally omitted.

---

Looks like I didn't pay attention to actual h size. It works now! Thanks!

Answer 1

Looks like some form of the Lotka-Volterra equation?

I'm not sure if if your initial condition is [300;8000] or [8000;300] (you specify it both ways above), but regardless, you have an oscillatory system that you're trying to integrate with a large fixed time step that is (much) greater than the period of oscillation. This is why your error explodes. If you try increasing n (say, 1e6), you'll find that eventually you'll get a stable solution (assuming that your Runge-Kutta implementation is otherwise correct).

Is there a reason why you're not using Matlab's builtin ODE solvers, e.g. ode45 or ode15s?

function ode45demo

[t,y]=odeode45(@f,[0 15],[300;8000]);

figure;
plot(t,y);

function ydot=f(t,y)
ydot(1,1) = 60*y(1) - 0.2*y(1)*y(2);
ydot(2,1) = 0.01*y(1)*y(2) - 100*y(2);

You'll find that adaptive step size solvers are much more efficient for these types of oscillatory problems. Because your system has such a high frequency and seems rather stiff, I suggest that you also look at what ode15s gives and/or adjust the 'AbsTol' and 'RelTol' options with odeset.