4 views (last 30 days)

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

James Tursa
on 8 Jan 2020

Edited: James Tursa
on 8 Jan 2020

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.

James Tursa
on 8 Jan 2020

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.