ODE solving using RK4 method (Predator prey)

function Q3
clc
clear all
close all 
 
y = zeros(1,1); 
x = zeros(1,1);
%Initial conditions 
t = 0.0; 
x(1) = 2; 
y(1) = 1;
%Parameters
h = 0.001; 
tfinal = 30.0; 
nsteps=(tfinal-t)/h; 
 
%store intial values 
tout(1) = t; 
xout(1) = x(1); 
yout(1) = y(1); 
 
for i=2:nsteps
    k1 = getf(t,y,x); 
    k2 = getf(t+h/2,y+h/2*k1,x+h/2*k1);
    k3 = getf(t+h/2,y+h/2*k2, x+h/2*k2);
    k4 = getf(t+h,y+h*k3, x+h/2*k3);
    y = y+h/6*(k1+2*k2 + 2*k3 + k4);
    x = x+h/6*(k1+2*k2 + 2*k3 + k4);
 
    t = t+h; 
    xout(i) = x(1); 
    yout(i) = y(1); 
    tout(i) = t; 
end 
 
plot(tout,yout,"x",tout,xout,"-r")
 
function f = getf(t,y,x)
alpha = 1.2;
beta = 0.6;
gamma = 0.8;
delta = 0.3;
f(1) = alpha*x(1)-beta*x(1)*y(1); 
f(2) = -gamma*x(1)+delta*y(1)*x(1);

There is something wrong with this because the solution isnt right

[John D'Errico]

Please do not remove your question, just saying the solution is not right. In fact, it is far more likely tjhat you just do not understand what was said. Of course, nobody can help you now, since we cannot see your question. This insults the person who tried to help you, and it makes your question useless for others who might try to understnd the problem later.

[Diya Saha]

I already got help for the question and it worked so I was really thankful for James. However since this was a hw question i was worried about plagarism and wasn't compfortable to have my full code on the web. Is that really so much of a problem?

[Stephen23]

Original question from Google Cache:

"ODE solving using RK4 method (Predator prey)"

function Q3
clc
clear all
close all 
 
y = zeros(1,1); 
x = zeros(1,1);
%Initial conditions 
t = 0.0; 
x(1) = 2; 
y(1) = 1;
%Parameters
h = 0.001; 
tfinal = 30.0; 
nsteps=(tfinal-t)/h; 
 
%store intial values 
tout(1) = t; 
xout(1) = x(1); 
yout(1) = y(1); 
 
for i=2:nsteps
    k1 = getf(t,y,x); 
    k2 = getf(t+h/2,y+h/2*k1,x+h/2*k1);
    k3 = getf(t+h/2,y+h/2*k2, x+h/2*k2);
    k4 = getf(t+h,y+h*k3, x+h/2*k3);
    y = y+h/6*(k1+2*k2 + 2*k3 + k4);
    x = x+h/6*(k1+2*k2 + 2*k3 + k4);
 
    t = t+h; 
    xout(i) = x(1); 
    yout(i) = y(1); 
    tout(i) = t; 
end 
 
plot(tout,yout,"x",tout,xout,"-r")
 
function f = getf(t,y,x)
alpha = 1.2;
beta = 0.6;
gamma = 0.8;
delta = 0.3;
f(1) = alpha*x(1)-beta*x(1)*y(1); 
f(2) = -gamma*x(1)+delta*y(1)*x(1);

There is something wrong with this because the solution isnt right

[Diya Saha]

@Stephen thankyou so much for doing that.

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[James Tursa]

Your basic problem is that you have two states, x and y, but your function arguments are inconsistent with this. Take this code:

k1 = getf(t,y,x); 

getf( ) returns a 2-element vector. So k1 is a 2-element vector. Yet you turn around and treat it like it is a scalar in the very next line:

k2 = getf(t+h/2,y+h/2*k1,x+h/2*k1);

I.e., you end up passing vectors into getf in those 2nd and 3rd input aguments when they should be scalars.

It is unclear to me which element is supposed to be x and which one is y. Since y is first in your argument list, let's assume f(1) is y and f(2) is x. Then the calculation for k2 should be something like this instead:

k2 = getf(t+h/2,y+h/2*k1(1),x+h/2*k1(2));

So all of your function calls will have to be fixed for this.

Having said all of that, I would advise carrying one state vector variable that contains both y and x, instead of carrying two separate y and x variables. That would simplify your code and also put your functions in a form that could be used by calls to the MATLAB supplied integrators such as ode45( ).

[Diya Saha]

Thankyou so much! The of k1 and all the other k's being a 2-element vector really helped me edit my code and I got the answer. Thakyou again!