Chapter 1

Notes from Chapter 1 of Hirsch, Smale, Devaney (HSD).

Simple Example of an ODE

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)= ke^{at}. \] Why? If you take the first derivative of the solution, then \(x'(t) = ake^{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 also specify the initial condition, that is \(x(t_0) = ke^{at_0} = u_0\). We are claiming that at \(t_0\), the value of the function is \(u_0\).

Rewriting this, we get \(k = u_{0}e^{-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'(t)=ax(t), \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:

t = np.linspace(-2, 2, 100)

def sol(t, k, a):
    return k * np.exp(a*t)
    
fig, axs = plt.subplots(1,3, figsize=(10,4), sharey=True, sharex=True)

for k in [-2, -1, 0, 1, 2]:
    axs[0].set_title("$a = 1$")
    axs[0].plot(t, sol(t, k, a=1))

for k in [-2, -1, 0, 1, 2]:
    axs[1].plot(t, sol(t, k, a=0))
    axs[1].set_title("$a = 0$")
    
for k in [-2, -1, 0, 1, 2]:
    axs[2].plot(t, sol(t, k, a=-1))
    axs[2].set_title("$a = -1$")

fig.suptitle("Solutions to $\\dot{x} = ax$ for different values of $k$")
fig.tight_layout()
plt.show()

An equilibrium is where \(x'=0\), i.e. where the dynamics come to a rest, denoted \(x_\ast\). For \(x'=ax\), if \(a\neq 0\) then the only equilibrium is \(x_\ast=0\). 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\), we have \(x'=0\) for all \(x\), so every point is an equilibrium and every solution is a constant line \(x(t)=k\).

The system is structurally stable when \(a\neq 0\), because for any \(\tilde a\) 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\).

Population Model

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 are more people. However, this is unbounded! So, we might want to add in something which makes the growth rate negative when the population hits a target \(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.

Note 1: Solution

We use separation to solve this

\[ \int \frac{dx}{x(1-x)}=\int a \,dt. \] Then, use the method of partial fractions

\[ \int \left(\frac{1}{x} + \frac{1}{1-x}\right) dx = \int a \, dt. \] Integrating \[ \begin{aligned} \int \left(\frac{1}{x} + \frac{1}{1-x}\right) dx &= \int a \, dt \\ \ln(x) - \ln(1-x)+K &= at \\ \ln\left(\frac{x}{1-x}\right)+K &= at \\ \ln\left(\frac{x}{1-x}\right)&= at -K\\ \left(\frac{x}{1-x}\right)&= \exp(at -K)\\ x&= \frac{\exp (at -K)}{1+\exp(at -K)}\\ x&= \frac{C\exp(at)}{1+C\exp(at)}\\ \end{aligned} \] If we algebraically solve for \(C\), we get \[ \begin{aligned} C\exp(at) &= \frac{x(t)}{1-x(t)} \\ C &= \frac{x(0)}{1-x(0)} \end{aligned} \]

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, 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}$")

fig.tight_layout()
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. If \(a\neq 0\), there are two equilibrium points: at \(x=0\) and at \(x=1\) (if \(a=0\), then \(x'=0\) for all \(x\) so every point is an equilibrium). 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()

Harvesting and Bifurcations

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!

Note 2: Solution

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: it’s a regime change from two equilibria to no equilibria.