Analysing Lotka and Volterra equation (ODE)

[yaya10]

Hello everyone,


I am working on a project on basically about ODE and phase plane

and I am working on this paper by Hanski

http://www.arctic-predators.uit.no/biblio_IPYappl/HanskiNature93 mustelid predators.pdf

How do i find numerical solutions which methods should I use ?

after finding the solution numerically, how can we tell whether the numbers in the table fit

the data described in table 1?

I am sorry i know you would have difficulty to understand me because i am not a native speaker:)

regards,

yaya

 

[jackmell]

Because of the conditonal, I recommend you write a custom Euler method which first doesn't include the conditional, then includes it. But that depends on how good you are. If you're not use to solving DEs numerically, then I recommend you build-up to this one. First code the Euler for y'=2y-1. Perfect that first. Then code a coupled system:

x'=2x-1.2 xy
y'=-y+1.2 xy

Get that one perfected. Then try modifying that code to include a conditional like:

$$x'=2x-1.2 xy$$
$$y'=\begin{cases}-y+1.2 xy & y>y_c \\ -y & y\leq y_c \end{cases}$$

Not sure what will happen there. But something like that. Then jump to your system first without the conditional, then with it.

 

[yaya10]

Thank you jackmell

Again , I have a question sorry to ask too many questions

The question is that when data has some different behaviour like

Chaos behaviour and limit cycles

what can we say about Chaos and limit cycle

I need to analyse the data but actually i get stuck with these terms which I do not know their meanings..!

Does Chaos mean the system is not stable

and Limit cycle the system will stabilise ??

Your explanation and help would highly appreciated

Tanks again

 

[HallsofIvy]

I would NOT recommend "Euler", that is much too inaccurate. A third order Runge-Kutta solver would be far superior.

 

[jackmell]

Chaos and limit cycles huh? You gotta' look those up dude. Chaos I don't believe means unstable although that may depend on how you define "stable". Consider the Lorenz attractor which is chaotic yet bounded. There are some properties of chaos that you should look up such as dense periodic points, ergodic behavior, and sensitivity. Do at least Wikipedia on them right? Then of course we can't talk about chaos without mentioning the magnus opus, "Chaos" by Peitegen. You got that book? If no, try and find it and go directly to chapter 10, "Deterministic Chaos" which reviews the Lorenz system and other chaotic systems. Work the examples there. Will take some time and you'll get a good understanding of what chaos is. Next get "Differential Equations" by Blanchard and Hall which has good examples on limit cycles or anothe text which covers them nicely. Limit cycles are exactly that: the system approaches a cyclic behavior at some limit, usually the long-term behavior approaches a limit. The Van der Pol equation first comes to mind when I think of limit cycles. Chaos however is when the system reaches a bounded behavior which is not a cycle nor is it at an equlibirum point. But perhaps you can define it better with some study. I think you should look at Van der Pol and the Lorenz system before studying yours. If you have time, do all that work first, then come back to your problem and it will be much easier for you to study and understand.

So my best advice to you is to put your project away for now and spend at least several days just researching chaos and limit cycles and actually working with them numerically and graphically, then come back to your problem.

 

[epenguin]

I could not open you link.
However it is not necessary to use the Euler method.
It is quite easy to integrate the Volterra-Lotka equations analytically.
They give cycles (start at any point and you will come back to it) around a centre - not limit cycles.
The more complicated problem with the conditional can be solved too.
You can solve the equations for the first condition (pure Lotka-Volterra). And you can solve it for the second condition even easier. Now how to put them together?

In the example of the condition of Jackmell, if y at the centre is less than y_c , then the solutions near enough the centre will be just what they are without the conditional.
Having equations of the solutions you should be able to work out at what initial conditions will your solution hit the conditional.

Beyond that? Take any general starting point (x₀, y₀) that lies on the conditional. From your known Lotka-Volterra solution work out at what x, y, the solution again hits the conditional. (In the example for given initial (x₀, y_c) at what other x on the solution cycle does y = y_c?). Do the same thing for solution of the equation of the second condition. If you find there is a starting point that for which both solutions hit the conditional at the same point then you're on a limit cycle. There may be one, or none, or even more than one such cycle. You will probably be able to work out mathematically whether for starting points of the limit cycle the trajectories spiral in towards the cycle (stable cycle) or away from it (unstable).

For plotting solutions of d.e.'s like this from Euler's or other method there are pnety of 'phase plane plotters' to be found on the web like this one: http://www.math.rutgers.edu/courses/ODE/sherod/phase-local.html

Fun. And the combination with analytical solutions nice.

 

[jackmell]

I didn't have problems clicking on the link and I read through it some more. I mean why didn't he just call them rats and ferrets? That's part of the problem in my opinion and feel he should have written it as:

$$\frac{dR}{dt}=r(t)R\left(1-\frac{R}{K(t)}\right)-\frac{c R F}{R+D}$$

$$\frac{dF}{dt}=\nu F\left(1-\frac{qF}{R}\right)H(R-R_c)+d(R,t)F H(R_c-R)$$

R(t) is the (prey) rat and F(t) is the predator (ferret).

with $H$ being the Heaveside function and:

$$d(R,t)=d_l f(t)H(R-R_c)+d_h g(t) H(R_c-R)$$

$$r(t)=r_0 S_a(t)$$
$$K(t)=K_0 S_b(t)$$

with the extra functions describing the dynamics the author is trying to model.

 

[epenguin]

OK I have been able to open and see that now. My comments were about the simple Lotka-Voterra equation: the paper you cite is a more complicated system.

You don't have to solve equations analytically, that is for fun and instruction with toy models. You can use the appelet I gave instead also for simple L-V. And you can use it for the more complicated eqs. 1, 2, 3 in the paper, also condensed into Jackmell's first two equations.

The appelet is suitable when you take time out of the equation, i.e. plot N against P, using for eqs.1 and 2

dN/dP = [rN(1 - N/K) - cPN/(N + D)]/vP(1 - qP/n)

and similarly you can get a different equation for dN/dP for N< N_crit combining eqs. 1 and 3. (I think you can solve that one analytically, not sure about 1 with 2.)

I think an important point to realize is that equations like Lotka-Volterra or the above in two variables, as the text mentions, cannot give rise to chaotic behaviour but rather, as the text mentions, can only give rise to convergence to s stable point or else a limit cycle (in special cases like simple L-V pure cycles around a centre). This arises essentially because the paths of d.e. solutions cannot cross, and if you ever do get back to a point where you were before then you have a cycle, all tied in with the fact the model is deterministic.

So I was a bit puzzled at them getting chaos out of the system until I read the paper carefully They say as well as the population variables they also include an independent time variable (expressed symbolically in Jackmell's last two equations). So you have a 3-variable system in which I am not surprised you can have chaos. I am not surprised either that you can have cycles. My understanding is that if you have a system e.g. spiralling towards a stable point, and the time of one or two turns of a spiral is a year then you can understand how an annual input (in their case a periodic change of rate constants) would turn the thing into a cycle. That is called 'entrainment' or at least my rough and perhaps defective understanding of it.

 

[yaya10]

I would NOT recommend "Euler", that is much too inaccurate. A third order Runge-Kutta solver would be far superior.

Thanks for your help

I have asked my supervisor he told me that he will help me with Matlab or Scilab

So, it would become much more easier if i master it :)

Thanks again HallsofIvy

 

[yaya10]

I would NOT recommend "Euler", that is much too inaccurate. A third order Runge-Kutta solver would be far superior.

Thanks for your help

I have asked my supervisor he told me that he will help me with Matlab or Scilab

So, it would become much more easier if i master it :)

Thanks again HallsofIvy

 

[yaya10]

Chaos and limit cycles huh? You gotta' look those up dude. Chaos I don't believe means unstable although that may depend on how you define "stable". Consider the Lorenz attractor which is chaotic yet bounded. There are some properties of chaos that you should look up such as dense periodic points, ergodic behavior, and sensitivity. Do at least Wikipedia on them right? Then of course we can't talk about chaos without mentioning the magnus opus, "Chaos" by Peitegen. You got that book? If no, try and find it and go directly to chapter 10, "Deterministic Chaos" which reviews the Lorenz system and other chaotic systems. Work the examples there. Will take some time and you'll get a good understanding of what chaos is. Next get "Differential Equations" by Blanchard and Hall which has good examples on limit cycles or anothe text which covers them nicely. Limit cycles are exactly that: the system approaches a cyclic behavior at some limit, usually the long-term behavior approaches a limit. The Van der Pol equation first comes to mind when I think of limit cycles. Chaos however is when the system reaches a bounded behavior which is not a cycle nor is it at an equlibirum point. But perhaps you can define it better with some study. I think you should look at Van der Pol and the Lorenz system before studying yours. If you have time, do all that work first, then come back to your problem and it will be much easier for you to study and understand.

So my best advice to you is to put your project away for now and spend at least several days just researching chaos and limit cycles and actually working with them numerically and graphically, then come back to your problem.

Thanks for your advice. actually I was a bit confusing about Chaos and Limit cycles but now

I know that they are just types of behaviours.
Spending sometime away form my project, has helped my a lot on my project.

Thnaks jackmell

 

[yaya10]

I could not open you link.
However it is not necessary to use the Euler method.
It is quite easy to integrate the Volterra-Lotka equations analytically.
They give cycles (start at any point and you will come back to it) around a centre - not limit cycles.
The more complicated problem with the conditional can be solved too.
You can solve the equations for the first condition (pure Lotka-Volterra). And you can solve it for the second condition even easier. Now how to put them together?

In the example of the condition of Jackmell, if y at the centre is less than y_c , then the solutions near enough the centre will be just what they are without the conditional.
Having equations of the solutions you should be able to work out at what initial conditions will your solution hit the conditional.

Beyond that? Take any general starting point (x₀, y₀) that lies on the conditional. From your known Lotka-Volterra solution work out at what x, y, the solution again hits the conditional. (In the example for given initial (x₀, y_c) at what other x on the solution cycle does y = y_c?). Do the same thing for solution of the equation of the second condition. If you find there is a starting point that for which both solutions hit the conditional at the same point then you're on a limit cycle. There may be one, or none, or even more than one such cycle. You will probably be able to work out mathematically whether for starting points of the limit cycle the trajectories spiral in towards the cycle (stable cycle) or away from it (unstable).

For plotting solutions of d.e.'s like this from Euler's or other method there are pnety of 'phase plane plotters' to be found on the web like this one: http://www.math.rutgers.edu/courses/ODE/sherod/phase-local.html

Fun. And the combination with analytical solutions nice.

thank you so much for your explanation and time


and the link is so cool. After my supervisor helps me with Matlab I will use the link you gave me to compare my solutions.

Thank you again epenguin

By the way, I like your nickname:)

 

[jackmell]

I didn't have problems clicking on the link and I read through it some more. I mean why didn't he just call them rats and ferrets? That's part of the problem in my opinion and feel he should have written it as:

$$\frac{dR}{dt}=r(t)R\left(1-\frac{R}{K(t)}\right)-\frac{c R F}{R+D}$$

$$\frac{dF}{dt}=\nu F\left(1-\frac{qF}{R}\right)H(R-R_c)+d(R,t)F H(R_c-R)$$

R(t) is the (prey) rat and F(t) is the predator (ferret).

with $H$ being the Heaveside function and:

$$d(R,t)=d_l f(t)H(R-R_c)+d_h g(t) H(R_c-R)$$

$$r(t)=r_0 S_a(t)$$
$$K(t)=K_0 S_b(t)$$

with the extra functions describing the dynamics the author is trying to model.

I'd like to make some improvements on these equations based on what the author is trying to model. I'll use a square-wave generator, S(t), with period one to model summer and winter. Summer is 1 and winter is -1. Then with H(t) being the Heaviside function, I can write (still rats and ferrets as I see it):

$$\frac{dR}{dt}=r(t)R\left(1-\frac{R}{K(t)}\right)-\frac{c R F}{R+D}$$

$$\frac{dF}{dt}=H(S(t))\left[\nu F\left(1-\frac{qF}{R}\right)H(R-R_c)+d_lF H(R_c-R)\right]+H(-S(t))\bigg[d_l H(R-R_c)+d_h F H(R_c-r)\bigg]$$

$$r(t)=r_a H(S(t))+H(-S(t)) r_b$$
$$K(t)=k_a H(S(t))+H(-S(t))k_b$$


I haven't run them so I'm not 100% sure they're doing what the author wants.

 

[yaya10]

Thank you so much for helping me

I have a meeting with my supervisor and he is going to give me some codes to run them on MATLAB and scilab. I hope they are easy to run because he wants me to analyse the data using one of them. I am happy to share the final results of my project with you guys, if you would lke t.

Thank you all so much.

Yaya

 

[jackmell]

Hi yaha,

If you want or can, could you please post your code here. I'd like to see how you guys coded the various conditionals and the other parameters. For example, here's mine in Mathematica based on the equations I posted above. Not sure if they're correct however. I used a square-wave generator for the summer-winter cycle but I'm not sure about what some of the paramater values should be.

rsummer = 5.4; 
rwinter = 2.7; 
ksummer = 100; 
kwinter = 25; 
c = 600; 
dlow = -0.1; 
dhigh = -5; 
thed = 15; 
q = 50; 
v = 2.8; 
rcrit = 2*thed; 
dsummer = -5; 
dwinter = -0.1; 
ratstart = 300; 
ferretstart = 25; 

r[t_] := rsummer*HeavisideTheta[SquareWave[t]] + rwinter*HeavisideTheta[-SquareWave[t]]; 

k[t_] := ksummer*HeavisideTheta[SquareWave[t]] + kwinter*HeavisideTheta[-SquareWave[t]]; 

rateqn = Derivative[1][rat][t] == r[t]*rat[t]*(1 - rat[t]/k[t]) - (c*rat[t]*ferret[t])/(rat[t] + thed);

ferreteqn = Derivative[1][ferret][t] == HeavisideTheta[SquareWave[t]]*
     (v*ferret[t]*(1 - (q*ferret[t])/rat[t])*HeavisideTheta[rat[t] - rcrit] + 
      dlow*ferret[t]*HeavisideTheta[rcrit - rat[t]]) + HeavisideTheta[SquareWave[t]]*
     (dlow*ferret[t]*HeavisideTheta[rat[t] - rcrit] + dhigh*ferret[t]*HeavisideTheta[rcrit - rat[t]]);

mysol = NDSolve[{rateqn, ferreteqn, rat[0] == ratstart, ferret[0] == ferretstart}, {rat, ferret}, {t, 0, 160}, 
   MaxSteps -> 100000];

Plot[{rat[t], ferret[t]} /. mysol, {t, 0, 105}, PlotRange -> All]

 

[epenguin]

I am happy to share the final results of my project with you guys, if you would lke t.

We would be interested when you are ready.

 

[yaya10]

Hi yaha,

If you want or can, could you please post your code here. I'd like to see how you guys coded the various conditionals and the other parameters. For example, here's mine in Mathematica based on the equations I posted above. Not sure if they're correct however. I used a square-wave generator for the summer-winter cycle but I'm not sure about what some of the paramater values should be.

rsummer = 5.4; 
rwinter = 2.7; 
ksummer = 100; 
kwinter = 25; 
c = 600; 
dlow = -0.1; 
dhigh = -5; 
thed = 15; 
q = 50; 
v = 2.8; 
rcrit = 2*thed; 
dsummer = -5; 
dwinter = -0.1; 
ratstart = 300; 
ferretstart = 25; 

r[t_] := rsummer*HeavisideTheta[SquareWave[t]] + rwinter*HeavisideTheta[-SquareWave[t]]; 

k[t_] := ksummer*HeavisideTheta[SquareWave[t]] + kwinter*HeavisideTheta[-SquareWave[t]]; 

rateqn = Derivative[1][rat][t] == r[t]*rat[t]*(1 - rat[t]/k[t]) - (c*rat[t]*ferret[t])/(rat[t] + thed);

ferreteqn = Derivative[1][ferret][t] == HeavisideTheta[SquareWave[t]]*
     (v*ferret[t]*(1 - (q*ferret[t])/rat[t])*HeavisideTheta[rat[t] - rcrit] + 
      dlow*ferret[t]*HeavisideTheta[rcrit - rat[t]]) + HeavisideTheta[SquareWave[t]]*
     (dlow*ferret[t]*HeavisideTheta[rat[t] - rcrit] + dhigh*ferret[t]*HeavisideTheta[rcrit - rat[t]]);

mysol = NDSolve[{rateqn, ferreteqn, rat[0] == ratstart, ferret[0] == ferretstart}, {rat, ferret}, {t, 0, 160}, 
   MaxSteps -> 100000];

Plot[{rat[t], ferret[t]} /. mysol, {t, 0, 105}, PlotRange -> All]

Hi jackmell,


My supervisor has not given me the codes yet.

Once he does, I will definitely share them with you.


Cheers,

Yaya

 

[yaya10]

We would be interested when you are ready.

Thats cool :)

yep when I ready, I will post them here or somewhere else.

:)


Cheers