5 views (last 30 days)

Show older comments

Hello,

I am currently working on the Lotka-Volterra model for modelling hare-lynx populations over time.

I have a working function to model this with ode45, that is not a problem.

These are the differential equations that I am using:

% growth rate of the prey population (snowshoe hare)

dxdt(1)=k1*x(1)-k2*x(1)*x(2);

% growth rate of the predator populations (canadian lynx)

dxdt(2)=k2*x(1)*x(2)-k3*x(2);

But then, I also have historical data of number of lynx and hares from 1900 to 1920. I was asked to estimate a new k1 value (different from my initial model, where I know all, k1, k2 and k3), then estimate the steady-state values for X (hare) and Y (lynx) from the historical data. And then, based on the new k1 and the calculated steady-states, calculate new k2 and k3.

So I have the new k1, which is 1.3610. But I don't know how to estimate the steady-states here from just looking at a graph (this was suggested by one of the tutors on a forum), or without knowing the new k2 and k3 yet. If I knew the k2 and k3, that would be fairly easy, as I would only need to calculate X and Y values for which my differential equations are zero. I read that there is an option of plotting nullclines to find steady-states, but I don't know how to apply them with only knowing k1 and having the historical values of number of prey and predators.

The plots above are the historical data that I have imported and plotted. I would very much appreciate if someone could tell me how I can estimate the steady-states from these.

BW,

K

David Goodmanson
on 16 Jul 2021

Edited: David Goodmanson
on 16 Jul 2021

Hi Kamila,

Here is one fairly simple way to find the parameters. The two equations are

dhare/dt = k1*hare - k2*hare*lynx

dlynx/dt = -k3*lynx + k2*hare*lynx

The steady state fixed point occurs when both time derivatives are zero in which case the fixed point is at [hare0,lynx0] and you can solve and get the standard result

hare0 = k3/k2 lynx0 = k1/k2

For Lotka-Volterra, the oscillations form a closed single-loop orbit. Fortunately for the Hudson Bay data someone chose to provide 21 years so there is an almost exactly closed two-loop curve. (The loops would be identical for perfect L-V behavior but of course this is real data).

Suppose 'hare' is a row vector of population by year. I added the first point to the end of the array with hare = [hare hare(1)] (similarly for lynx) so now if you plot hare vs lynx you will get a closed curve with two loops.

So far so good. Now suppose you time integrate both sides of the first L-V equation around an integral number of orbits, going from tstart to tfin. On the left hand side, integrating the derivative you get hare(tstart)-hare(tfin) = 0 since the start position and the end position are the same. On the right hand side you get

0 = k1*Ihare - k2*Iboth where Ihare = Int hare dt and Iboth = Int hare*lynx dt

so

Iboth/Ihare = k1/k2 = lynx0

and similarly for the second equation

Iboth/Ilynx = k3/k2 = hare0

For the integrations you can use t = 0:21 since you have one-year intervals and use the trapz function ,

Ihare = trapz(t,hare) (similary for lynx); Iboth = trapz(t,hare.*lynx)

I got lynx0 = 2.0883e+04 hare0 = 3.5292e+04

which, given the look of the data, ought to be rounded to something like 21000 and 35000. Once you have those, with a new k1 you can use the eqns above to solve for k2 and k3.

A long time ago there was a comic strip about Fat Freddy's Cat. One day Fat Freddy's Cat went out toward the woods and met a big kind of cat he had never seen before. He said, what kind of cat are you? And the cat said, "I'm a lynx". To which FFC said, "Oh, I always wondered where links sausage came from".

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!