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:NW(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
Best Answer