5. Differential and difference forms¶
Difference equation model:
Differential equation model:
Recall that the delta sign (\(\Delta\)) means change in or the difference. Compare the difference equation with the differential form, which uses the terminology \(dN\) and \(dt\). These represent infinitesimally small time steps, corresponding to our common-sense perception of time as divisible ever more finely without limit. In differential equations populations change smoothly rather than in finite steps—growth approximating that of organisms that can reproduce at any time, such as bacterial or human populations.
It turns out that differential equations are harder for computers to solve than difference equations. Computers cannot make infinitely fine time steps, but have to approximate by using very small time steps instead. On the other hand, difference equations can be harder to solve mathematically.
r <- 1; s <- -0.001; N <- 1; t <- 0; dt <- 1
population <- data.frame(time = integer(), size = integer())
time = seq(from = t, to = 20, by = dt)
for (t in 1:(1 + length(time))) {
population[t, ] = c(t, N)
dN = (r + s * N) * N; N = N + dN;
if(N < 0) N = 0;
}
population
| time | size | |
|---|---|---|
| <dbl> | <dbl> | |
| 1 | 1 | 1.000000 |
| 2 | 2 | 1.999000 |
| 3 | 3 | 3.994004 |
| 4 | 4 | 7.972056 |
| 5 | 5 | 15.880558 |
| 6 | 6 | 31.508924 |
| 7 | 7 | 62.025036 |
| 8 | 8 | 120.202967 |
| 9 | 9 | 225.957181 |
| 10 | 10 | 400.857715 |
| 11 | 11 | 641.028522 |
| 12 | 12 | 871.139478 |
| 13 | 13 | 983.394966 |
| 14 | 14 | 999.724273 |
| 15 | 15 | 999.999924 |
| 16 | 16 | 1000.000000 |
| 17 | 17 | 1000.000000 |
| 18 | 18 | 1000.000000 |
| 19 | 19 | 1000.000000 |
| 20 | 20 | 1000.000000 |
| 21 | 21 | 1000.000000 |
| 22 | 22 | 1000.000000 |
Above is computer code for a difference equation presented earlier, which leveled off at 1000, but with an additional if() statement to catch negative population size. If the population is far above its carrying capacity, the calculation could show such a strong decline that the next year’s population would be negative—meaning that the population would die out completely. The highlighted addition just avoids projecting negative populations, which are mathematically possible but biologically meaningless. Below is similar code for the corresponding differential equation, with differences highlighted.
N <- 1; t <- 0
dt = 1 / (20 * 24 * 60)
time = seq(from = 0, to = 20, by = dt)
population <- data.frame(time = integer(), exponential = integer())
for(t in 1:(1 + length(time))){
population[t, ] = c(t, N)
dN = (r + s * N) * N * dt; N = N + dN
if(N < 0) N = 0
}
head(population,100)
Warning
Feel free to run the above loop on your own machines, but note that we have defined our time vector to contain every single second from the beginning of day 0 to the end of day 20. It will take some time for R calculate all the time steps.
This intends to model infinitely small time steps. Of course it cannot do that exactly, but must settle for very small time steps. Instead of \(dt = 1\), for example, representing one year, it is set here to about one second, dividing the time period (20 days) by 24 hours, creating 20 * 24 * 60 time steps. This is hardly infinitely small, but for populations of bacteria and humans it is close enough for practical purposes. Still, it is important to check for negative populations in case the time step is not small enough.
How small is close enough to infinitely small? is the question. To find out, you can set the time step to something small and run the code, which will produce a set of population values through time. Then set the step smaller still and run the code again. It will run more slowly because it is calculating more steps, but if essentially the same answer appears—if the answer “converges”—then you can make the step larger again, speeding the calculation. With a few trials you can find a time step that is small enough to give accurate answers but large enough to allow your code to run reasonably fast.
Fig. 5.1 Figure 5.1 Differential logistic growth (maroon) compared with discrete (green). No dots appear on the differential form, since that represents infinitesimal time steps, whereas the difference form has a dot at each point calculated.¶
Figure 5.1 shows the results of running the differential equation version of the program (the second one above, in maroon) versus the difference equation version (the first above, in green). The differential equation has the same parameters and general shape, but the population approaches its carrying capacity more quickly. Because the differential time steps produce offspring earlier—not waiting for the end of the step—offspring are available to reproduce earlier, and so forth.
This particular method for differential equations is called “Euler’s method” (pronounced “Oiler’s”), a basic approach not often used in the twentieth century because computers not so long ago were millions of times slower than they are now. Today this method is often fast enough, and is desirable because of its relative simplicity.
5.1. A mathematical view¶
Differential equations can be amenable to mathematical analysis. To repeat, here is the differential population model.
It turns out there is something simple about infinity, and when time steps are infinitely small the methods of calculus developed over the centuries can solve this differential equation exactly, mathematically. If you apply a symbolic mathematics computer package, or the methods for integration of functions developed in calculus, you can find the population value N for any future time t. This is called the “solution” to the differential equation.
Most differential equations cannot be solved this way but, fortunately, the basic equations of ecology can. This solution becomes useful in projecting forward or otherwise understanding the behavior of a population. If you know the starting N, s, and r, you can plug them into the formula to find the population size at every time in the future, without stepping through the differential equation.
5.2. Exponential solution¶
To think about this mathematically, first set \(s\) to zero, meaning no density dependence. The differential equation then reduces to \(\frac{dN}{dt} = r\), and if you replace \(s\) in the equation above with 0 you get this:
A measure often used for exponential growth, and that we will apply later in this book, is “doubling time”—the time that must elapse for the population to double. For exponential growth, this is always the same, no matter how large or small the population. For exponential growth, the equation above is:
\(N_0\) is the starting population at time 0, \(N(t)\) is the population at any time \(t\), and \(r\) is the constant growth rate—the “intrinsic rate of natural increase.” How much time, \(\tau\), will elapse before the population doubles? At some time \(t\), the population will be \(N(t)\), and at the later time \(t+\tau\), the population will be \(N(t+\tau)\). The question to be answered is this: for what \(τ\) will the ratio of those two populations be 2?
Substituting the right-hand side of the exponential growth equation from above gives:
The factor \(N_0\) cancels out, and taking the natural logarithms of both sides gives:
Since the log of a ratio is the difference of the logs, this yields:
Since logarithms and exponentials are inverse processes—each one undoes the other—the natural logarithm of \(e^{x}\) is simply \(x\). That gives:
And finally, the doubling time, \(\tau\), is:
In other words, the doubling time for exponential growth, where \(r\) is positive and \(s\) is 0, is just the natural logarithm of 2 (0.69314718…) divided by the growth rate \(r\).
5.3. Logistic solution¶
Recall that the carrying capacity is \(-\frac{s}{r}\), also called \(K\). So, wherever \(-\frac{s}{r}\) appears, substitue \(K\) as follows:
This is the solution given in textbooks for logistic growth. There are slight variations in ways it is written, but they are equivalent.
5.4. Orthologistic solution¶
Finally, let \(s\) be positive. This creates a vertical asymptote and orthologistic growth. The position in time of the vertical asymptote is the “singularity” mentioned earlier. The interesting question is, when \(s\) is positive, what is the time of the singularity? That is, when will the population grow beyond all bounds in this model?
What must happen to the denominator for the population to grow to unbounded values? It has to get closer and closer to 0, for then the \(N(t)\) will grow closer and closer to infinity. So to find the singularity, you only have to set the denominator to 0, and then solve for the time \(t\). You can go through the intermediate steps in the algebra below, or use a mathematical equation solver on your computer to do it for you.
Recall the equation from above for solving the differential equation:
Setting the denominator to 0 will lead along this algebraic path:
Multiply through by \(\frac{r}{s}e^{rt}\) to obtain:
Next take the logarithms of both sides, remembering that \(\ln e^{x} = x\):
Finally, divide through by \(r\) to find the time of the singularity:
In the 1960s, Heinz von Foerster wrote about this in the journal Science. Though the consequences he suggested were deadly serious, his work was not taken very seriously at the time, perhaps in part because the time was so far away (about a human lifetime), but perhaps also because he put the date of the singularity on Friday the 13th, 2026, his 115th birthday. In the title of his paper he called this “doomsday”, when the human population would have demolished itself.
Von Foerster used a more complicated model than the \(r+sN\) model we are using, but it led to the same result. Some of the ideas were picked up by Paul Ehrlich and others, and became the late-1960s concept of the “population bomb”—which was taken seriously by many.
Fig. 5.2 The Population Bomb¶