222
8 Solving Ordinary Differential Equations
The corresponding differential equation becomes
N
= r(N)N .
The reader is strongly encouraged to repeat the steps in the derivation of the Forward
Euler scheme and establish that we get
N
n+1
= N
n
+ Δt r(N
n )N
n ,
which computes as easy as for a constant r, since r(N n ) is known when computing
N n+1 . Alternatively, one can use the Forward Euler formula for the general problem
u = f (u, t) and use f (u, t) = r(u)u and replace u by N.
The simplest choice of r(N) is a linear function, starting with some growth value
¯
r and declining until the population has reached its maximum, M, according to the
available resources:
r(N) = ¯
r(1 − N/M) .
In the beginning, N M and we will have exponential growth e ¯
rt , but as N
increases, r(N) decreases, and when N reaches M, r(N) = 0 so there is no more
growth and the population remains at N(t) = M. This linear choice of r(N) gives
rise to a model that is called the logistic model. The parameter M is known as the
carrying capacity of the population.
Let us run the logistic model with aid of the ode_FE function. We choose N(0) =
100, Δt = 0.5 month, T = 60 months, r = 0.1, and M = 500. The complete
program, called logistic.py, is basically a call to ode_FE:
from ode_FE import ode_FE
import matplotlib.pyplot as plt
for dt, T in zip((0.5, 20), (60, 100)):
u, t = ode_FE(f=lambda u, t: 0.1*(1 - u/500.)*u, \
U_0=100, dt=dt, T=T)
plt.figure() # Make separate figures for each pass in the loop
plt.plot(t, u, ’b-’)
plt.xlabel(’t’); plt.ylabel(’N(t)’)
plt.savefig(’tmp_{:g}.png’.format(dt))
plt.savefig(’tmp_{:g}.pdf’.format(dt))
Figure 8.9 shows the resulting curve. We see that the population stabilizes around
M = 500 individuals. A corresponding exponential growth would reach N 0 e rt =
100e 0.1·60 ≈ 40, 300 individuals!
What happens if we use “large” Δt values here? We may set Δt = 20 and
T = 100. Now the solution, seen in Fig. 8.10, oscillates and is hence qualitatively
wrong, because one can prove that the exact solution of the differential equation is
monotone.
8 Solving Ordinary Differential Equations
The corresponding differential equation becomes
N
= r(N)N .
The reader is strongly encouraged to repeat the steps in the derivation of the Forward
Euler scheme and establish that we get
N
n+1
= N
n
+ Δt r(N
n )N
n ,
which computes as easy as for a constant r, since r(N n ) is known when computing
N n+1 . Alternatively, one can use the Forward Euler formula for the general problem
u = f (u, t) and use f (u, t) = r(u)u and replace u by N.
The simplest choice of r(N) is a linear function, starting with some growth value
¯
r and declining until the population has reached its maximum, M, according to the
available resources:
r(N) = ¯
r(1 − N/M) .
In the beginning, N M and we will have exponential growth e ¯
rt , but as N
increases, r(N) decreases, and when N reaches M, r(N) = 0 so there is no more
growth and the population remains at N(t) = M. This linear choice of r(N) gives
rise to a model that is called the logistic model. The parameter M is known as the
carrying capacity of the population.
Let us run the logistic model with aid of the ode_FE function. We choose N(0) =
100, Δt = 0.5 month, T = 60 months, r = 0.1, and M = 500. The complete
program, called logistic.py, is basically a call to ode_FE:
from ode_FE import ode_FE
import matplotlib.pyplot as plt
for dt, T in zip((0.5, 20), (60, 100)):
u, t = ode_FE(f=lambda u, t: 0.1*(1 - u/500.)*u, \
U_0=100, dt=dt, T=T)
plt.figure() # Make separate figures for each pass in the loop
plt.plot(t, u, ’b-’)
plt.xlabel(’t’); plt.ylabel(’N(t)’)
plt.savefig(’tmp_{:g}.png’.format(dt))
plt.savefig(’tmp_{:g}.pdf’.format(dt))
Figure 8.9 shows the resulting curve. We see that the population stabilizes around
M = 500 individuals. A corresponding exponential growth would reach N 0 e rt =
100e 0.1·60 ≈ 40, 300 individuals!
What happens if we use “large” Δt values here? We may set Δt = 20 and
T = 100. Now the solution, seen in Fig. 8.10, oscillates and is hence qualitatively
wrong, because one can prove that the exact solution of the differential equation is
monotone.