Runge kutta 4 with two ODE's - function inside a function

Hi there I am using a Runge Kutta-4 method in order to solve a predator prey model. My problem that is I do not understand how to implement the W inside the y1 in the right way.

the differential equations are:

dy1/dx = (a - alpha*y2)*y1 - W(t)*y1,

dy2/dx (-c+gamma*y1)*y2

Can please somebody help me???

% parameters
a = 3;
c = 3;
alpha = 2*10^(-3);
gamma = 7*10^(-3);

% define time
ti = 0;
tf = 5;
delta_t = 0.5;

t = ti:delta_t:tf;
h = delta_t;
N = ceil(tf/h);

% define the seasonal fishing load factor
W = zeros(11,1);

for n =1:N
W(n,1)=fishing_load_factor(t(n));
end


% Define functions
fy1 = @(t,y1,y2,W) (a-alpha*y2)*y1-W*y1;
fy2 = @(t,y1,y2) (-c+gamma*y1)*y2;

%initial conditions
y1(1)=600;
y2(1)=1000;
t(1)=0;

% Update loop
for i = 1:N
    %Update time
    t(1+i)=t(i)+h;
    
    % Update y1 and y2
    K11 = fy1(t(i)      ,y1(i)        , y2(i), W(i));
    K12 = fy2(t(i)      ,y1(i)        , y2(i));
    
    K21 = fy1(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K11, W(i));
    K22 = fy1(t(i)+h/2  ,y1(i)+h/2*K12, y2(i)+h/2*K12);
    
    K31 = fy1(t(i)+h/2  ,y1(i)+h/2*K21, y2(i)+h/2*K21, W(i));
    K32 = fy1(t(i)+h/2  ,y1(i)+h/2*K22, y2(i)+h/2*K22);
 
    K41 = fy1(t(i)+2    ,y1(i)+h*K31   , y2(i)+h*K31, W(i));
    K42 = fy1(t(i)+2    ,y1(i)+h*K32  , y2(i)+h*K32);
    
    y1(1+i) = y1(i)+h/6*(K11 + 2*K21 + K21*K31 + K41);
    y2(1+i) = y2(i)+h/6*(K12 + 2*K22 + K22*K32 + K42);
end

where W is a function

function W = fishing_load_factor(t)

if 0 <= t & t < 3/12
    W=0;
elseif 3/12 <= t & t <8/12
    W=2;
elseif 8/12 <= t & t < 1
    W=0;
elseif 1 <= t
    W=fishing_load_factor(t-1)
end

Answer 1

[James Tursa]

Since W is a function of time, it needs to match the time that you are using in each particular line of code. E.g., take this line:

K11 = fy1(t(i)      ,y1(i)        , y2(i), W(i));

The time for this line is t(i), so the W will be fishing_load_factor(t(i)). E.g.,

K11 = fy1(t(i)      ,y1(i)        , y2(i), fishing_load_factor(t(i)));

For the above lines you get the same result, but now look at this line:

K21 = fy1(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K11, W(i));

Here the time is different, but you didn't adjust your W accordingly. It should be this instead:

K21 = fy1(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K11, fishing_load_factor(t(i)+h/2));

Similarly for the other lines of code involving W.

In addition, the times of these lines are incorrect:

    K41 = fy1(t(i)+2    ,y1(i)+h*K31   , y2(i)+h*K31, W(i));
    K42 = fy1(t(i)+2    ,y1(i)+h*K32  , y2(i)+h*K32);

The times should be t(i)+h instead of t(i)+2.

Then, the factors for the K31 and K32 terms in these lines are incorrect:

    y1(1+i) = y1(i)+h/6*(K11 + 2*K21 + K21*K31 + K41);
    y2(1+i) = y2(i)+h/6*(K12 + 2*K22 + K22*K32 + K42);

Instead of K21*K31 and K22*K32 they should be 2*K31 and 2*K32 respectively.

FInally, you are using the derivative calculations with the wrong states. y1 derivatives should be calculated with the fy1 function and applied to the y1 state, and y2 derivatives should be calculated with the fy2 function and applied to the y2 state. But look what you have done with these lines:

    K21 = fy1(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K11, W(i));
    K22 = fy1(t(i)+h/2  ,y1(i)+h/2*K12, y2(i)+h/2*K12);

K11 came from a fy1 calculation, so it should only be used with the y1 state, but you are using it with the y2 state in this expression:

    y2(i)+h/2*K11

That expression should be using the derivative that came from the fy2 function, namely K12:

    y2(i)+h/2*K12

And then your K22 calculation, which will be used with the y2 state, uses the fy1 function instead of the fy2 function. These two lines should be this instead:

    K21 = fy1(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K12, fishing_load_factor(t(i)+h/2));
    K22 = fy2(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K12);

So, LOTS of these types of errors in your code that need to get corrected. Just follow this rule:

y1 derivatives always use the fy1 function and y1 always gets propagated with the Kx1 derivatives

y2 derivatives always use the fy2 function and y2 always gets propagated with the Kx2 derivatives

I would add that it looks like your h is way too big to pick up the change points in your W(t) function properly. You will need to shrink h by a fair amount to make sure you pick up those W(t) change points adequately.

[Wytse Petrie]

Well thank you very much for your clear explaination. It was easy to correct the script with your help. This is the end result, which is correct :)

clear all 
close all
clc

%% sub-assignment 1
%chosen for RK4, it is the only who fullfill the three criteria that are
%suggested in the question. RK4 is the only explicitit method that may be
%stable for systems which have imaginairy eigenvalues.

%% sub-assignment 2
% parameter
a = 3;
c = 3;
alpha = 2*10^(-3);
gamma = 7*10^(-3);

% define time
ti = 0;
tf = 5;
delta_t = 0.5;

t = ti:delta_t:tf;
h = delta_t;
N = ceil(tf/h);

% define the seasonal fishing load factor
W = zeros(11,1);

for n =1:N
W(n,1)=fishing_load_factor(t(n));
end


% Define functions
fy1 = @(t,y1,y2,W) (a-alpha*y2)*y1-W*y1;
fy2 = @(t,y1,y2) (-c+gamma*y1)*y2;

%initial conditions
y1(1)=600;
y2(1)=1000;
t(1)=0;

% Update loop
for i = 1:N
    %Update time
    t(1+i)=t(i)+h;
    
    % Update y1 and y2
    K11 = fy1(t(i)      ,y1(i)        , y2(i), fishing_load_factor(t(i)));
    K12 = fy2(t(i)      ,y1(i)        , y2(i));
    
    K21 = fy1(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K12, fishing_load_factor(t(i)+h/2));
    K22 = fy2(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K12);
    
    K31 = fy1(t(i)+h/2  ,y1(i)+h/2*K21, y2(i)+h/2*K22, fishing_load_factor(t(i)+h/2));
    K32 = fy2(t(i)+h/2  ,y1(i)+h/2*K21, y2(i)+h/2*K22);
 
    K41 = fy1(t(i)+h    ,y1(i)+h*K31   , y2(i)+h*K32, fishing_load_factor(t(i)+h));
    K42 = fy2(t(i)+h    ,y1(i)+h*K31  , y2(i)+h*K32);
    
    y1(1+i) = y1(i)+h/6*(K11 + 2*K21 + 2*K31 + K41);
    y2(1+i) = y2(i)+h/6*(K12 + 2*K22 + 2*K32 + K42);
    
    
end

[Wytse Petrie]

I'm now calculating the global error truncation of this Ringe Kutta 4 method. Therefore i need to calculate the find y1 and y2 with dx*2. I now got the following code, but it won't work.

clear all
close all
clc

% parameter
a = 3;
c = 3;
alpha = 2*10^(-3);
gamma = 7*10^(-3);

% define time
ti = 0;
tf = 5;
delta_t0 = 1;

t0 = ti:delta_t0:tf;
h0 = delta_t0;
N0 = ceil(tf/h0)+1;

% define the seasonal fishing load factor
W = zeros(N0,1);

for n =1:N0
W(n,1)=fishing_load_factor(t0(n));
end


% Define functions
fy1 = @(t,y1,y2,W) (a-alpha*y2)*y1-W*y1;
fy2 = @(t,y1,y2) (-c+gamma*y1)*y2;

%initial conditions
y1(1)=600;
y2(1)=1000;
t0(1)=0;

% Update loop
for i = 1:N0-1
    %Update time
    t0(1+i)=t0(i)+h0;
    
    % Update y1 and y2
    K11 = fy1(t0(i)      ,y1(i)        , y2(i), fishing_load_factor(t0(i)));
    K12 = fy2(t0(i)      ,y1(i)        , y2(i));
    
    K21 = fy1(t0(i)+h0/2  ,y1(i)+h0/2*K11, y2(i)+h0/2*K12, fishing_load_factor(t0(i)+h0/2));
    K22 = fy2(t0(i)+h0/2  ,y1(i)+h0/2*K11, y2(i)+h0/2*K12);
    
    K31 = fy1(t0(i)+h0/2  ,y1(i)+h0/2*K21, y2(i)+h0/2*K22, fishing_load_factor(t0(i)+h0/2));
    K32 = fy2(t0(i)+h0/2  ,y1(i)+h0/2*K21, y2(i)+h0/2*K22);
 
    K41 = fy1(t0(i)+h0    ,y1(i)+h0*K31   , y2(i)+h0*K32, fishing_load_factor(t0(i)+h0));
    K42 = fy2(t0(i)+h0    ,y1(i)+h0*K31   , y2(i)+h0*K32);
    
    y1(1+i) = y1(i)+h0/6*(K11 + 2*K21 + 2*K31 + K41);
    y2(1+i) = y2(i)+h0/6*(K12 + 2*K22 + 2*K32 + K42);
    
    
end

[James Tursa]

In both of your new codes, the h values being used are way too big. You are using h values of 0.5 and 1.0, but look at your fishing_load_factor( ) function which has change points at 3/12 and 8/12. Neither of your h values will be picking up those change points properly because they are way too big. I would try dropping your h value to much less than the 3/12 and 8/12 values, maybe something like 0.01 or so. As it is, your results are meaningless with regards to that fishing_load_factor( ) function.