Notes from Chapter 1 of Hirsch, Smale, Devaney (HSD)
The most common differential equation is
\[ x'(t) = ax(t), \] which gives us the relationship between a function \(x(t)\) and its derivative \(x'(t)\). Finding a solution means finding a function (or a general class of functions) that satisfies this relationship. Here, the solution is \[ x(t)= k\exp(at). \] Why? If you take the derivative, then \(x'(t) = ak\exp(at) = ax(t)\). Thus, we satisfy the original ODE.
This is actually a general solution because \(k\) is not determined. To determine \(k\), we must specify the initial condition, that is \(x(t_0) = k\exp(at_{0}) = u_0\). Rewriting this, we get \(k = u_{0}\exp(-at_0)\). So, to determine \(k\), we need to know both \(u_0\) and \(t_0\). Often, we just set \(t_{0}=0\), resulting in \(k=u_0\). Then, the initial value problem is written as: \[ x'=ax, \quad x(t_{0}=0)=u_0. \] Generally, any solution must (1) satisfy the differential equation and (2) satisfy the initial value.
The constant \(a\) is a parameter of the equation and it affects how the solution looks. As an example, if \(a=0\), then our solution is \(x(t)=k\), a constant function. Not that interesting. Otherwise, we will get the plots below:
import matplotlib.pyplot as plt
import numpy as np
t = np.linspace(-2, 2, 100)
def f(t, k, a):
return k * np.exp(a*t)
fig, axs = plt.subplots(1,2, figsize=(8,4), sharey=True, sharex=True)
for k in [-2, -1, 0, 1, 2]:
axs[0].set_title("$a = 1$")
axs[0].plot(t, f(t, k, a=1))
for k in [-2, -1, 0, 1, 2]:
axs[1].plot(t, f(t, k, a=-1))
axs[1].set_title("$a = -1$")
plt.suptitle("Solutions to $\\dot{x} = ax$ for different values of $k$")
plt.show()
The equilibrium point is where \(x'=0\). It is a point where the dynamics come to a rest. In this example, it happens to occur at \(x=0\). And, when \(a>0\), solutions move away from the equilibrium point, making it a source. When \(a<0\), solutions move towards it, making the equilibrium point a sink. When \(a=0\), the solution is a constant line \(x(t) = k\) (not plotted).
The system is structurally stable when \(a\neq 0\), because for any \(b\) with the same sign as \(a\), the solutions will behave similarly (diverge or converge to equilibrium). However, at \(a=0\), the smallest change in \(a\), leads to a very different solution. Thus, there is a bifurcation at \(a=0\).
If \(a>0\), we can use the solution of \(x'=ax\) as a population model. The assumption is that the rate of growth (\(x'\)) depends on the original population (\(x\)). That is, population grows faster as there is more population. However, this is unbounded! So, we might want to add in something which makes the growth rate negative when the population hits \(N\): \[ x' = ax(1 - \frac{x}{N}) \] The parameter \(a\) then gives rate of growth when \(x\) is small, and \(N\) gives a “carrying capacity”. When \(x\) is small, the ODE behaves as \(x' \approx ax\). Only when \(x\) is large does the curbing kick in! For simplicity, we can set \(N=1\), such that the “ideal” size is 1. This simplifies down to
\[ x' = ax(1-x) \] This is a nonlinear first-order ODE. Nonlinear because we depend on a nonlinear function of \(x\). First-order because we only work with the first derivative.
From the above derivation, the solution to \(x'(t) = ax(t)(1-x(t))\) is
\[ x(t) = \frac{x(0)\exp(at)}{1-x(0)+x(0)\exp(at)}. \]
Below, we plot the solution and the rate of population growth for different values of \(x(0)\) and \(a\). In the first row, we vary the size of the initial population \(x(0)\) while keeping \(a\) constant. In the second row, we vary the growth rate \(a\) while keeping \(x(0)\) constant.
def x(t, x0, a):
return (x0 * np.exp(a*t)) / (1 - x0 + x0 * np.exp(a*t))
def f(t, x0, a):
return a*x(t, x0, a)*(1 - x(t, x0, a))
t = np.linspace(0, 3, 100)
fig, axs = plt.subplots(2,2, figsize=(10,6), sharex=True)
for x0 in [0.1, 0.3, 0.9, 1.5]:
a = 2
axs[0,0].set_title(f"Solution when $a = {a}$")
axs[0,0].plot(t, x(t, x0, a=a), label=f"$x(0)={x0}$")
axs[0,0].set_xlabel("$t$")
axs[0,0].set_ylabel("$x(t)$")
axs[0,0].legend()
axs[0,1].set_title(f"Rate of population growth when $a = {a}$")
axs[0,1].plot(t, f(t, x0, a=a), label=f"$x(0)={x0}$")
axs[0,1].set_xlabel("$t$")
axs[0,1].set_ylabel("$x'(t)$")
axs[0,1].legend()
for a in [0.5, 1, 2]:
x0 = 0.4
axs[1,0].set_title(f"Solution when $x(0) = {x0}$")
axs[1,0].plot(t, x(t, x0=x0, a=a), label=f"$a={a}$")
axs[1,0].legend()
axs[1,1].set_title(f"Rate of population growth when $x(0) = {x0}$")
axs[1,1].plot(t, f(t, x0=x0, a=a), label=f"$a={a}$")
plt.show()
There is another useful plot, the “phase portrait”, which shows the rate of change (\(x'\)) as a function of the population size (\(x\)). This simply shows the ODE relationship without solving it. Regardless of the value of \(a\), there are two equilibrium points: at \(x=0\) and at \(x=1\). Recall that equilibrium points refer to when the dynamics come to a rest (\(x'=0\)).
When combined with the plots above, we can see that when \(a>0\), the equilibrium at \(x=0\) is a source (solutions move away from it) and the equilibrium at \(x=1\) is a sink (solutions move towards it).
x_vals = np.linspace(-0.1, 1.2, 100)
plt.title(f"Phase portraits across $a$")
plt.xlabel("$x$")
plt.ylabel("$x'$")
plt.axhline(0, color='black', lw=0.5, ls='--')
plt.axvline(0, color='black', lw=0.5, ls='--')
for a in [0.5, 1, 2]:
plt.plot(x_vals, a*x_vals*(1 - x_vals), label=f"$a={a}$")
plt.legend()
plt.show()
The ODE we have been working with is:
\[ x' = ax(1-x) \]
Let’s change a couple things. For simplicity, let’s set \(a=1\) and introduce a constant “harvesting” factor \(h\)
\[ x' = x(1-x) - h, \quad \text{s.t. } h \geq 0. \]
Again, we could solve this via integration. See below!
Just intuitively, thinking about the phase portrait (\(x'\) vs. \(x\)) all we have done is added a vertical shift parameter. So, it is not surprising that this affects where the equilibria occur.
x_vals = np.linspace(-0.1, 1.2, 100)
a = 1
plt.title(f"Phase portraits across $h$")
plt.xlabel("$x$")
plt.ylabel("$x'$")
plt.axhline(0, color='black', lw=0.5, ls='--')
plt.axvline(0, color='black', lw=0.5, ls='--')
for h in [0, 0.25, 0.5]:
plt.plot(x_vals, x_vals*(1 - x_vals) - h, label=f"$h={h}$")
plt.legend()
plt.show()
In the case that \(h < \frac{1}{4}\), we have two equilibria. When \(h=\frac{1}{4}\), we encounter a single equilibrium. Any value of \(h\) beyond that results in no equilibrium point for the ODE. Most importantly, \(h=\frac{1}{4}\) represents a bifurcation point: its a regime change from two equilibria to no equilibria.
Let’s be a bit more realistic with our model:
\[ x' = ax(1-x) - h(1+\sin(2\pi t)), \]
where we reintroduce \(a\) and add periodicity to our harvesting model. Moreover, this harvesting is time-dependenting making this a nonautonomous differential equation.
Unfortunately, this is not something we can solve analytically. So, we turn to qualitative observations. Interestingly, because of the periodicity, we can solve for all points \(0 \leq t \leq 1\), and this would give us solutions for all points in time.
Let’s say we had a hypothetical map \(p\) which takes an initial condition \(x(0)\) and gives us the value of \(x(1)\). That is, \(p(x(0)) = x(1)\). If we composed this function twice, we would get \(p(p(x(0))) = x(2)\). If we composed it \(n\) times, we would get \(p^{n}(x(0)) = x(n)\). This map is called the Poincaré map. However, these are usually really hard to compute.
A flow is a function \(\phi(t, x_0)\), or \(\phi_t(x_0)\), which gives the solution at time \(t\) given an initial condition \(x_0\). The notation just emphasizes that the solution depends on both time and the initial condition.