Lotka-Volterra model Euler's Method

[Robben]

Homework Statement

*I am not sure if this should be in the computer science section or here?

I am trying to graph the densities, of the Lotka-Volterra prey and predator model, as a function of time, i.e. $p(t)$ vs $t$ and $q(t)$ vs $t$. Also, the phase space, i.e. $p$ vs $q$, but I am not sure how to do that?

Let $p$ be the prey density and $q$ is the predator density, thus: $$\frac{dp}{dt} = ap\left(1-\frac{p}{K}\right)-\frac{bpq}{1+bp}$$ $$\frac{dq}{dt}=mq\left(1-\frac{q}{kp}\right).$$

Where $bpq$ is the interaction rate between the species, $\frac{bpq}{1+bp}$ is the effective rate of eating prey, $m$ is the mortality rate of the predators, $K$ and $k$ are the carrying capacitance of each population. Let $a=0.2, \ m=0.1, \ K=500,\ k = 0.2.$

Homework Equations

$$p_i = p_{i-1}+\text{p-slope}_{i-1}\Delta t$$
$$q_i = q_{i-1}+\text{q-slope}_{i-1}\Delta t.$$

The Attempt at a Solution

We are also given $a=0.2, \ m=0.1, \ K=500,\ k = 0.2, b = 0.1, p(0) = 10, q(0) = 5.$ I chose $\Delta t = h = 0.1.$

Thus:

$$\begin{align} p_i& = p_{i-1} + \left(0.2~ p_{i-1}\left(1 - \dfrac{p_{i-1}}{500}\right) -\dfrac{0.1~ p_{i-1}~ q_{i-1}}{1 + 0.1~ p_{i-1}} \right)(0.1) \\ q_i& = q_{i-1} + \left(0.1~q_{i-1}\left(1 - \dfrac{q_{i-1}}{0.2~ p_{i-1}}\right) \right)(0.1) \end{align}$$

Which gives the iteration:

$p_0 = 10, \ q_0 = 5$
$p_1 = 9.946, \ q_1 = 4.925$
$p_2 = 9.89538, \ q_2 = 4.85231$
$p_3 = 9.84803, \ q_3 = 4.78187$
$\ldots$

My java code:

public class Lotka {
    public static void main(String[] args) {
        double[] p = new double[1000];
        double[] q = new double[1000];
        p[0] = 10;
        q[0] = 5;
        for(int i = 1; i < 1000; i++) {
            p[i] = p[i-1] + (0.2*p[i-1]*(1-p[i-1]/500)-(0.2*p[i-1]*q[i-1])/(1+0.2*p[i-1]))*0.001;
            q[i] = q[i-1] + (0.1*q[i-1]*(1-q[i-1]/(0.2*p[i-1])))*0.001;
            System.out.println("" + p[i] + " " + q[i]);  
        }   
    }   
}

 

[soarce]

The phase space $(p,q)$ will containt trajectories having $t$ as parameter, i.e. you form pairs of $p$ and $q$ at equal time $t$ and then you plot them. Did you plot $p(t)$ and $q(t)$ ?

 

[Robben]

The phase space $(p,q)$ will containt trajectories having $t$ as parameter, i.e. you form pairs of $p$ and $q$ at equal time $t$ and then you plot them. Did you plot $p(t)$ and $q(t)$ ?

So, I should make the parameter $t$ vary instead of letting $t = 0.1$? So, something like $0.1$ through $0.5$?

 

[soarce]

Let's run $t$ from $t_{min}$ and $t_{max}$, from your program you will obtain three columns of data $t_i$, $p_i$ and $q_i$. The curves in the phase space are given by $(p_i,q_i)$.

 

[tms]

The first thing you should do is replace all the numbers in your code with named variables. Not only does this make your code more readable and maintainable, but you will find that you have used different values in your code from what you say they should be.

You should also reformat your equations in the code to add spaces, again for readability; as it is now, they are completely unreadable.

Also, since you are only using the last value of p and q, there is no reason to store all of them. Just assign the initial values to variables like p_old and q_old, and then at the end of each iteration store the new values in those variable.

 

[Robben]

Let's run $t$ from $t_{min}$ and $t_{max}$, from your program you will obtain three columns of data $t_i$, $p_i$ and $q_i$. The curves in the phase space are given by $(p_i,q_i)$.

I see, thank you!

You should also reformat your equations in the code to add spaces, again for readability; as it is now, they are completely unreadable.

Also, since you are only using the last value of p and q, there is no reason to store all of them. Just assign the initial values to variables like p_old and q_old, and then at the end of each iteration store the new values in those variable.

What do you mean by add spaces?

Also, I need all values of i because I am going to plot them in gnuplot.

 

[tms]

What do you mean by add spaces?

I mean space out the equations so they are more readable. Instead of

q[i] = q[i-1] + (0.1*q[i-1]*(1-q[i-1]/(0.2*p[i-1])))*0.001;

write something like

q[i] = q[i-1] + (0.1 * q[i-1] * (1 - q[i-1] / (0.2 * p[i-1]))) * 0.001;

It's easier to see how the various terms group that way.

Even better would be to replace the numbers with variables:

q = q_old + (m * q_old * (1 - q_old / (k * p_old))) * delta_t;

Also, I need all values of i because I am going to plot them in gnuplot.

Just redirect the output to a file, and have gnuplot read that file.

 

[Robben]

I mean space out the equations so they are more readable. Instead of

q[i] = q[i-1] + (0.1*q[i-1]*(1-q[i-1]/(0.2*p[i-1])))*0.001;

write something like

q[i] = q[i-1] + (0.1 * q[i-1] * (1 - q[i-1] / (0.2 * p[i-1]))) * 0.001;

It's easier to see how the various terms group that way.

Even better would be to replace the numbers with variables:

q = q_old + (m * q_old * (1 - q_old / (k * p_old))) * delta_t;

Just redirect the output to a file, and have gnuplot read that file.

I see, I have edited my code using your suggestions. Also, I have redirected the output file to gnuplot and got the following graph (prey vs predator) that I inserted and uploaded below. I varied the parameter $b$ from $0.01$ to $0.5$. Is that how the graph suppose to look of the Lotka model because the graph I am seeing on google when I search for Lotka model looks different.

public class Lotka {
    public static void main(String[] args) {
        double[] prey = new double[1000];
        double[] predator = new double[1000];  
        double[] b = new double[101]; //parameter
        prey[0] = 10;
        predator[0] = 5;  
        b[0] = 0.01;   
  
        for (int j = 1; j < 101; j ++) {
            b[j] = b[0] * ((double) j/2);  
        }   
   
        for(int i = 1; i < 100; i++) {
            prey[i] = prey[i-1] + (0.2 * prey[i-1] * (1 - prey[i-1] / 500)
                -(b[i] * prey[i-1] * predator[i-1]) / (1 + b[i] * prey[i-1])) * .001;
            predator[i] = predator[i-1] + (0.1 * predator[i-1] * (1 - predator[i-1]
                /(0.2 * prey[i-1]))) * 0.001;  
            System.out.println("" + predator[i] + " " + prey[i]);
        }  
    }      
}

 

[tms]

It's better, but there are still problems. For one thing, all the parameters, not just b, should be named. (I'm sure Java has a way to name constants, but I don't remember what it is. In C it would be "const double b;" or whatever the variable is. At any rate, if you can't make them constants, make them variables; it is much easier to read and modify that way.)

I don't understand why you made b vary. The parameters should be constant for any particular trial.

In your description, you had $\Delta t$ be 0.1, but in your code it is 0.001. The larger value seems more appropriate, given how little change there is from one iteration to another. With a hundred time more time, you may see a better plot.

You still don't need to save the results of each iteration, since you aren't using them.

 

[Robben]

It's better, but there are still problems. For one thing, all the parameters, not just b, should be named. (I'm sure Java has a way to name constants, but I don't remember what it is. In C it would be "const double b;" or whatever the variable is. At any rate, if you can't make them constants, make them variables; it is much easier to read and modify that way.)

I don't understand why you made b vary. The parameters should be constant for any particular trial.

In your description, you had $\Delta t$ be 0.1, but in your code it is 0.001. The larger value seems more appropriate, given how little change there is from one iteration to another. With a hundred time more time, you may see a better plot.

You still don't need to save the results of each iteration, since you aren't using them.

I need to vary b because it serves as the interaction rate, so I want to see what happens as we increase $b$. I noticed that the densities die off when $b = 0.5$ which is why I was varying $b$.

I had to change delta t because if I left it at $0.1$ then it wouldn't converge which is why i changed it to $0.001$.

 

[tms]

If you want to vary $b$ like that, you have to run different trials, each with a different value of $b$. In anyone trial, $b$ must stay the same (unless, of course, a varying $b$ is part of your model, which it isn't here).

Regarding $\Delta t$, I don't know what you mean by "converge". If you want to keep it at it's current value, then you will need many more iterations per trial.

 

[Robben]

If you want to vary $b$ like that, you have to run different trials, each with a different value of $b$. In anyone trial, $b$ must stay the same (unless, of course, a varying $b$ is part of your model, which it isn't here).

Regarding $\Delta t$, I don't know what you mean by "converge". If you want to keep it at it's current value, then you will need many more iterations per trial.

I see, thank you!