1

I am trying to model the Lotka-Volterra Predator-Prey system, using the coupled DE's:

dy(1)/dt = rx(1-x/k) - ay(1)y(2) % prey population

dy(2)/dt = aby(1)y(2) - dy(2) % predator population

Here is the code I have :

% Solves equations using numerical ODE solver 45 (nonstiff runge kutta)

runtime = 1000; % Duration time of simulation in seconds.

%Parameter values used in simulation 

r = 0.5; % exponentional growth rate of prey in absence of predator

a = 0.01;% conversion efficiency of predator

b = 0.02; % attack rate

d = 0.10; % death rate

k = 750; % carrying capacity

y0 = [10, 10] % initial conditions y(1)= 10, y(2) = 10

deq1=@(t,y) [r.*y(1)*(1-(y(1)./k))- a.*y(1)*y(2); a.*b.*y(1).*y(2)-    d.*y(2)];

[t,sol] = ode45(deq1,[0 runtime],y0);

How can I change my code so that time series plots (y(1) vs. time) show oscillations, rather than steady states? It seems like something is going wrong in the integration step, as the plots are structured the way I want them but the behavior of the functions is not what I expected.

Steady state time series plots

Steady state time series plots

Phase plot reaching equilibrium

Phase plot reaching equilibrium

Community
  • 1
  • 1
  • Looks like the equilibrium point is not a center but a saddle point. Compute the signature of the Jacobian at that point. – Lutz Lehmann Oct 12 '16 at 09:02
  • So, would error in my code could be due to my parameter values or do you think I have used ode45 incorrectly? – Charlie Robinson Oct 13 '16 at 04:09
  • See https://stackoverflow.com/q/35071393/3088138 for an example that actually oscillates and moves towards the equilibrium. Thus it may be that your equilibrium is still a center, but the dampening is so strong that only 3/4 of an oscillation occur. – Lutz Lehmann Oct 13 '16 at 18:35

0 Answers0