Math216 Introduction to Differential Equations

Deniz Bilman, Department of Mathematics, University of Michigan

Lecture 6: Autonomous equations and population dynamics
Return of the qualitative analysis

We covered several equations today in the setting that $y=\phi(t)$ is the population of a species.

The first model for the change of this population assumed that the rate of change of the population is proportional to the current size of the population. In other words, $\frac{dy}{dt}\sim y$. This yields the following initial value problem:

$$ \frac{dy}{dt}=r y, ~y(0)=y_0, $$

where $r$ is a proportionality constant. The solution to this IVP is given by

$$ \phi(t) = y_0 e^{rt}. $$

Therefore, for any initial condition $y_0>0$ (positivity due to the fact that $y$ is the population of a species), the corresponding solution

  • grows exponentially as time elapses (as $t$ increases) if $r \gt 0$
  • decays to $0$ exponentially as time elapses (as $t$ increases) if $r \lt 0$.

The second model had the hypothesis that the porportionality constant $r$ used above in fact depends on the current value of $y$ as well:

$$ \frac{dy}{dt}=h(y) y. $$

The rate function $h(y)$ should satisfy:

  • $h(y)\approx r$ for small values of $y$,
  • $h(y)$ decreases as $y$ increases (population grows at a slower rate as it gets larger, maybe due to lack of resources?),
  • $h(y)<0$ if $y$ is larger than a certain value (if the population becomes too big, it starts to decay).

The simplest function that satisfies all of the assumptions above is

$$ h(y) = r- ay, $$

for some constant $a \gt 0$. Using this, we have the following ODE

$$ \frac{dy}{dt} = (r-ay)y, $$

which can also be rewritten in the form

$$ \frac{dy}{dt} = r\left(1-\frac{y}{K}\right)y, $$

where $K=\frac{r}{a} \gt 0$. This equation is called the logistic equation. At a glance, we can see that it has two equilibrium solutions: $y=\phi_1(t)=0$ and $y=\phi_2(t)=K$. Recall that these solutions are the zeros of the function on the right hand side of the ODE:

$$ f(y)=r\left(1-\frac{y}{K}\right)y $$
  • Since $f'(0) \gt 0$, $y=\phi_1(t)$ is an unstable equilibrium solution (phase portrait drawn in class).
  • $f'(K) \lt 0$, therefore $y=\phi_2(t)$ is an asymptotically stable equilibrium solution.

Therefore, given any initial condition $y(0)=y_0>0$, the corresponding solution has the saturation property:

$$ y(t)\to K,\text{ as }t\to+\infty. $$

For this reason, the constant $K$ is called the saturation level, or environment carrying capacity.

Example

Consider $y'=\frac{1}{2}(1-\frac{y}{10})$. Observe that $K=10$. The integral curves can be seen in the plot below.


In [1]:
using PlotlyJS;
# # using Plots
# # plotly()
const timespan = 0:0.1:20.0;
const initialdata = 0.0:1:15;
const r = 0.5;
const K = 10.0;
tracesODE1 = GenericTrace[];
for q0 in initialdata
  trace = scatter(x=timespan, y=q0*K./(q0+(K-q0)*exp(-r*timespan)), name = "y(0)=$(q0)");
  push!(tracesODE1, trace);
end
plot(tracesODE1, Layout(title="Integral Curves for the solution of the logistic equation", xaxis=attr(title="t"),yaxis=attr(title="y(t)")))


Plotly javascript loaded.

To load again call

init_notebook(true)

Out[1]:

Remark: Note that if $K$ is taken to be large, the nonlinear term (the $y^2$ term) in the logistic equation has a small coefficient:

$$ y' = r\left(1-\frac{y}{K}\right)y = ry -\frac{1}{K}y^2. $$

Neglecting the nonlinear term results in the ODE $y'=ry$, whose solution grow exponentially in time since $r>0$. All of the solutions are unbounded in time in absence of the nonlinear term in the ODE. However, if even a very weak nonlinearity is introduced in the equation, the solutions exhibit a dramatically different behavior: they saturate at $y=K$, no matter how big $K$ is. The solution of the logistic equation is bounded for any initial value.

Now consider the following variant of the logistic equation, which is obtained by introducing a minus sign in the right hand side:

$$ \frac{dy}{dt}=-r\left(1-\frac{y}{T}\right)y, $$

where $r>0$ and $T>0$ are constants. Not note that the zeros of $f(y)=-r\left(1-\frac{y}{T}\right)y$ are $y=\phi_1(t)=0$ and $y=\phi_2(t)=T$, just as was the case for the logistic equation. However, the nature of these equilibrium solutions is different (reversed compared to the logistic equation.)

  • Since $f'(0) \lt 0$, $y=\phi_1(t)=0$ is an asymptotically stable equilibrium solution (phase portrait drawn in class).
  • $f'(T) \gt 0$, therefore $y=\phi_2(t)=T$ is an unstable equilibrium solution.

These imply that for initial conditions $0 \lt y_0 \lt T$, the corresponding solution has the property that $y(t)\to 0$ as $t\to \infty$. The solution gets away from $y=T$ and approaches $y=0$ as time elapses. On the other hand, if the initial value $y_0$ is greater than $T$, then the solution $y(t)$ grows exponentially in time. For this reason, the equilibrium solution $y=T$ is called the threshold level. A population whose size is initially less than $T$ shrinks to $0$ as $t\to\infty$. A population whose size is initially greater than $T$ keeps growing in time.

Example

Consider $y'=-2(1-\frac{y}{10})$. Observe that $T=10$. The integral curves can be seen in the plot below. We only plot solutions corresponding to initial conditions $0\leq y_0 \leq 10$, which are the solutions that decay to $0$ in time.


In [2]:
using PlotlyJS;
# # using Plots
# # plotly()
const timespan = 0:0.1:10;
const initialdata = 0.0:1.0:10.0;
const r = 2.0;
const T = 10.0;
tracesODE1 = GenericTrace[];
for q0 in initialdata
  trace = scatter(x=timespan, y=q0*T./(q0+(T-q0)*exp(r*timespan)), name = "y(0)=$(q0)");
  push!(tracesODE1, trace);
end
plot(tracesODE1, Layout(title="Integral Curves", xaxis=attr(title="t"),yaxis=attr(title="y(t)")))


WARNING: redefining constant timespan
WARNING: redefining constant initialdata
WARNING: redefining constant r
Out[2]:

We now plot a solution that corresponds to $y_0 \gt 10$, we take $y_0=11$. Observe that there is something more than exponential growth is happenning to this solution. It seems to develop a vertical asymptote at some finite value of $t$. We can see that it happens near $t\approx 1.198$. You can mouse-over and read the data points.


In [3]:
using PlotlyJS;
# # using Plots
# # plotly()
const timespan = 0:0.001:1.198;
const initialdata = 10.0:1.0:11.0;
const r = 2.0;
const T = 10.0;
tracesODE1 = GenericTrace[];
for q0 in initialdata
  trace = scatter(x=timespan, y=q0*T./(q0+(T-q0)*exp(r*timespan)), name = "y(0)=$(q0)");
  push!(tracesODE1, trace);
end
plot(tracesODE1, Layout(title="Integral curve that blows up in finite time", xaxis=attr(title="t"),yaxis=attr(title="y(t)")))


WARNING: redefining constant timespan
WARNING: redefining constant initialdata
Out[3]:

Indeed, we can solve the IVP

$$ \frac{dy}{dt} = r\left(1-\frac{y}{T} \right)y,~$y(0)=y_0$ $$

under the assumption that $y_0 \gt T$ and obtain the solution

$$ y(t) = \frac{y_0 T}{y_0 + (T-y_0)e^{rt}}. $$

Now, given the values of the positive constants (parameters) $r$ and $T$ we can see that the denominator of this formula becomes 0 for a finite value of $t$:

$$ y_0 + (T-y_0)e^{rt^*} = 0 \iff t^* = \frac{1}{r}\log\frac{y_0}{y_0-T}. $$

Therefore the solution corresponding to the initial condition $y_0\gt T$ has a vertical asymptote at

$$ t = \frac{1}{r}\log\frac{y_0}{y_0-T}. $$

This means that the population becomes infinite in finite amount of time. Going back to the last example, if we plug in the values $r=2$ and $T=10$, we can compute $t^*$ using this formula and verify that the result agrees with what we have seen in the plot $t^*\approx 1.198$.


In [4]:
r = 2.0;
T = 10.0;
y0 = 11.0;
tstar = (1/r)*log(y0/(y0-T));
print(tstar)


1.1989476363991853

Although the threshold mechanism can be widely seen in real life scenarios, the fact that the corresponding solution becomes infinite in finite time is unrealistic. We would like to have the threshold level so that small disturbances gets washed out with time, whereas significant disturbances change the state of a quantity and cause it to saturate at some other state. In other words, we would like to combine the saturation property caused by the dynamics of the logistic equation with the threshold property caused by the dynamics of the variant of the logistic equation. This can be done by considering the ODE

$$ \frac{dy}{dt} = -r\left(1-\frac{y}{T}\right)\left(1-\frac{y}{K}\right)y, $$

where $r$, $T$, and $K$ are positive constants and $T \lt K $. Now the nonlinearity is cubic, and you should peform the phase portrait analysis. The equilibrium solutions are $y=\phi_1(t)=0$, $y=\phi_2(t)=T$, and $y=\phi_3(t)=K$. Below is an example with parameters $r=1$, $T=8$, and $K=12$. We can see both the threshold and the saturation mechanisms in action.


In [5]:
using ODE;
using PlotlyJS;
const timespan = 0:0.1:20.0;
tracesODE1 = GenericTrace[];
function rhsODE(t,y)
  return -1.0*(1-y/8.0)*(1-y/12.0)*y
end;
const initialdata=0.0:1.0:15.0;
for j in 1:length(initialdata)
  y0 = initialdata[j]
  t,y = ode45(rhsODE,y0,timespan)
  trace = scatter(x=t, y=y, name = "y(0)=$(initialdata[j])");
  push!(tracesODE1, trace);
end

plot(tracesODE1, Layout(title="Integral Curves for logistic growth with a threshold", xaxis=attr(title="t"),yaxis=attr(title="y(t)")))


WARNING: redefining constant timespan
WARNING: redefining constant initialdata
Out[5]:

Acknowledgements

I would like to thank the following people for correcting a mistake in these notes:

  • Ahmed Nafis Arafat

In [ ]: