1.2.1 Self-competition

One improvement would be to demand that we include a notion of self-competition within the population. This could, for example, model competition for food or territory. Mathematically, we need a decay term which is small for small \(x\), allowing the population to grow, but dominates the growth term when \(x\) gets larger, thus restricting its growth.

The simplest example of such a model is the logistic equation. Originally introduced by Pierre François Verhulst in 1838, the equation is nonlinear and takes the form \[\begin{equation} \label{logistic} \mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = a x\left(1-\frac{x}{K}\right), \end{equation}\] where \(a>0\) is the usual growth term and, as we shall see, \(K\) is the limiting population or carrying capacity.

Can we guess what this model looks like? See that there is an equilibrium (\(\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = 0\)) at \(x=K\). The term \((1-x/K)\) is negative if \(x>K\) and positive if \(x<K\) so we might expect it to either decay towards \(K\) from above and up towards \(K\) if from below.

This sort of analysis turns out to be quite common for population models because it is not always possible to solve them analytically. Thankfully this time we can, however.

1.2.2 Solutions to the logistic equation

Let’s go ahead and now work out the solutions of our model. We can separate [logistic]: \[\begin{equation} \int\frac{1}{x(1 - x/K)}\,\mathrm{d}x= at + C. \end{equation}\] We can integrate the first integral using partial fractions (remember those!), \[\begin{equation} \frac{1}{x(1 - x/K)} = \frac{1}{x} + \frac{1}{K\left(1-x/K\right)}, \end{equation}\] and so \[\begin{equation} \log(x) - \log(1-x/K) =a t+C \implies \frac{x}{1-x/K} = A\mathrm{e}^{at}, \end{equation}\] so that finally, \[\begin{equation} \label{logsol} x(t) = \frac{A \mathrm{e}^{a t}}{1+ \frac{A}{K}\mathrm{e}^{at}}. \end{equation}\]

1.3 demonstrates the behaviour of this equation.

What about that limiting population we promised? See that as \(t \to \infty\), \(x(t) \to K\): a fact independent of the initial condition.

1.2.3 Equilibria of the logistic equation

Rather than worry about how the population changes, we might only really care where it will end up, given sufficient time. We have already spotted that the system tends to a state where \(\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}=0\) at \(x=K\). Is this the only possibility?

No! If we set \(\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}=0\) then we have \[\begin{equation} a x\left(1-\frac{x}{K}\right) = 0 \implies x=0, K. \end{equation}\] So there is also an unchanging state of \(x=0\) where there is no population. This makes sense, of course: no population means no reproduction. At this stage we make our first definition:

Equilibrium:

An equilibrium of a system is one in which all time derivatives are zero. For example, consider the system \[\begin{equation} \mathchoice{\frac{\mathrm{d}^2 u}{\mathrm{d}t^2}}{\mathrm{d}^2 u/\mathrm{d}t^2}{\mathrm{d}^2 u/\mathrm{d}t^2}{\mathrm{d}^2 u/\mathrm{d}t^2} + \mathchoice{\frac{\mathrm{d}u}{\mathrm{d}t}}{\mathrm{d}u/\mathrm{d}t}{\mathrm{d}u/\mathrm{d}t}{\mathrm{d}u/\mathrm{d}t} = u^2+v^2,\qquad \mathchoice{\frac{\mathrm{d}^{3} v}{\mathrm{d}t^{3}}}{\mathrm{d}^{3} v/\mathrm{d}t^{3}}{\mathrm{d}^{3} v/\mathrm{d}t^{3}}{\mathrm{d}^{3} v/\mathrm{d}t^{3}} = uv. \end{equation}\] The equations of equilibrium are \[\begin{equation} 0 = u^2+v^2,\qquad 0 = uv \end{equation}\] (the only solution to which is \(u=v=0\)). The definition of equilibrium is often (in a dynamical systems context) referred to as a steady state or fixed point.

Permissible/feasible equilibrium:

We earlier demanded that the population is \(\geq 0\). We thus define a permissible or feasible equilibrium as one which satisfies this criteria. It will be important throughout the course that our models have permissible equilibria to be valid. Another idea we will come to discuss is that a good model should not allow a positive initial population to become negative.

1.2.4 A first look at stability

We note that for any positive \(A\), [logsol] will tend to \(K\). We say that the \(x=K\) equilibrium is stable because any small change from \(x=K\) (say \(x= K-\varepsilon\)) will tend back to \(x=K\) if we go forward in time (convince yourself of this by looking again at 1.3). However, if we are at \(x=0\) and there is a sudden small change to \(x=\varepsilon\), perhaps representing a small population migration, if we go forward in time it will grow inexorably towards \(K\). Thus we say that the \(x=0\) equilibrium is unstable.

But did we need to solve the logistic equation, [logistic], to find this? Actually no! Because our differential equation is of the form \(\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = f(x)\), we can use a common technique where we simply plot \(\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}\) against \(x\) (known as plotting the phase space) and make some observations.

Look at 1.4. Our two equilibria are marked. If you start at a value of \(x\) where \(\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}\) is positive, we know that \(x\) increases in time, so you move to the right over time (indicated by the forward-pointing red arrow). And where \(\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}\) is negative, you move to the left (the backward-pointing red arrow). This instinctively tells us that you will, at \(t=\infty\), always end up at \(x=K\) unless you start exactly at \(x=0\). Do you agree?

We will make a more mathematically precise notion of stability/instability of equilibria in 2 and 2.2.

1.4.1 Interactions between predator and prey

Instead of introducing self-competition, a second possibility to avoid unbounded growth is to model a second population, \(y(t)\), which represents a second species. In the case of our greenfly population, we have our antagonists in ladybirds (1.6). Since ladybirds prey on greenflies, the greenfly population, \(x\), will decrease proportionally to \[\begin{equation*} \text{[the number of ladybirds, $y$]} \times \text{[the number of greenflies, $x$],} \end{equation*}\] i.e., the number of interactions of the two species which may lead to a sad little greenfly funeral. This law will therefore be in the form \[\begin{align} \mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = a x - bxy, \end{align}\] with \(b\) the rate at which fatal interactions occur.

But we must then also model the changing ladybird population, \(y(t)\). We assume in the absence of greenflies it will decrease as its food supply has vanished: \[\begin{equation} \mathchoice{\frac{\mathrm{d}y}{\mathrm{d}t}}{\mathrm{d}y/\mathrm{d}t}{\mathrm{d}y/\mathrm{d}t}{\mathrm{d}y/\mathrm{d}t} = -cy. \end{equation}\] (Ask yourself: Why is this proportional to \(y\)?) However, it will also grow proportionally to the availability of food, i.e., interactions of the two species (at some rate \(d\)), so \[\begin{equation} \mathchoice{\frac{\mathrm{d}y}{\mathrm{d}t}}{\mathrm{d}y/\mathrm{d}t}{\mathrm{d}y/\mathrm{d}t}{\mathrm{d}y/\mathrm{d}t} = -cy + d x y. \end{equation}\] So we have a coupled set of ordinary differential equations, \[\begin{align} \label{lotvol} \mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} &= ax - bxy,\\ \mathchoice{\frac{\mathrm{d}y}{\mathrm{d}t}}{\mathrm{d}y/\mathrm{d}t}{\mathrm{d}y/\mathrm{d}t}{\mathrm{d}y/\mathrm{d}t} &= -cy + dxy, \end{align}\] where, once again, \(ax\) represents prey reproduction with unlimited food, \(-cy\) is predator natural death or emigration, and \(bxy\) and \(dxy\) are interaction terms. Note that \(b\) and \(d\) are not necessarily the same: many prey may have to be eaten for the predator population to grow by one.

This system represents both the individual growth/decay of the species (self-interactions) as well as their mutual interaction. This is the so-called Lotka–Volterra (predator–prey) system developed by the Italian mathematician Vito Volterra (1860–1940) in 1926 to explain the fluctuation of fish populations in the Adriatic Sea (more on Lotka in 4). In more modern theories there are multiple species, each with their own interactions, but we will limit ourselves to this simpler but highly instructive classical system.

An example solution is shown for the parameters \((a,b,c,d) = (\frac23,\frac43,1,1)\), \(x(0)=1\), \(y(0)=1\) in 1.7(a). See how the peaks in the greenfly population naturally increase the ladybird food supply: the ladybird population then increases. In turn this leads to the greenfly population dropping as they get eaten, and this decrease in food supply leads to the ladybird population dropping as food becomes competitive.

This periodic behaviour is made clear using a phase plot, as shown in 1.7(b). In this case we have a parametrised plot \((x(t),y(t))\): a geometric plot of the variables of the system (2D here because we have two variables). Closed curves in phase space indicate a periodic relationship between the two parameters.

1.4.2 Parameter reduction of Lotka–Volterra

The parameters \((a,b,c,d)\) play a key role in determining the system’s behaviour. However, they are not all independent. We can work out how many we actually need by doing a nondimensionalisation.

Thus far we’ve dealt directly with the dimensional form of the differential equations, [lotvol], meaning that the parameters in the equations have relevant dimensions (or units) associated with them. While this is perhaps convenient for making direct predictions from measured values of the parameters, working with dimensional parameters can at best keep the equations looking ‘untidy’, and at worst obscure the mathematical structure of the equations or relevant approximations that can be made.

The process of nondimensionalisation is to separate our variables into a nondimensional bit and a dimensional bit: \[\begin{equation} x = \widehat{x}X, \quad y = \widehat{y}Y, \quad t = \widehat{t}T. \end{equation}\] Here, \(x\), \(X\), \(y\) and \(Y\) have units of ‘individuals’, and \(t\) and \(T\) have units of time. Variables with hats indicate they are nondimensional.

If we use this substitution, our equations are \[\begin{align} \mathchoice{\frac{\mathrm{d}\widehat x}{\mathrm{d}\widehat t}}{\mathrm{d}\widehat x/\mathrm{d}\widehat t}{\mathrm{d}\widehat x/\mathrm{d}\widehat t}{\mathrm{d}\widehat x/\mathrm{d}\widehat t}\frac{X}{T} &= a\widehat{x}X - b\widehat{x}\widehat{y}XY,\\ \mathchoice{\frac{\mathrm{d}\widehat y}{\mathrm{d}\widehat t}}{\mathrm{d}\widehat y/\mathrm{d}\widehat t}{\mathrm{d}\widehat y/\mathrm{d}\widehat t}{\mathrm{d}\widehat y/\mathrm{d}\widehat t}\frac{Y}{T} &= -c\widehat{y}Y + d\widehat{x}\widehat{y}XY. \end{align}\] Rearranging so that there are no dimensions on the left hand sides, we have \[\begin{align} \mathchoice{\frac{\mathrm{d}\widehat x}{\mathrm{d}\widehat t}}{\mathrm{d}\widehat x/\mathrm{d}\widehat t}{\mathrm{d}\widehat x/\mathrm{d}\widehat t}{\mathrm{d}\widehat x/\mathrm{d}\widehat t} &= a\widehat{x}T - b\widehat{x}\widehat{y}YT,\\ \mathchoice{\frac{\mathrm{d}\widehat y}{\mathrm{d}\widehat t}}{\mathrm{d}\widehat y/\mathrm{d}\widehat t}{\mathrm{d}\widehat y/\mathrm{d}\widehat t}{\mathrm{d}\widehat y/\mathrm{d}\widehat t} &= -c\widehat{y}T + d\widehat{x}\widehat{y}XT. \end{align}\] So now we pick \(T\), \(X\) and \(Y\) to remove as many of our parameters as possible. If we choose \[\begin{equation} \label{TYX-scaling} T = \frac{1}{a}, \quad Y = \frac{a}{b}, \quad X = \frac{c}{d}, \end{equation}\] then the system can be written as \[\begin{align} \label{reducedlotvol} \mathchoice{\frac{\mathrm{d}\widehat{x}}{\mathrm{d}\widehat{t}}}{\mathrm{d}\widehat{x}/\mathrm{d}\widehat{t}}{\mathrm{d}\widehat{x}/\mathrm{d}\widehat{t}}{\mathrm{d}\widehat{x}/\mathrm{d}\widehat{t}} &= \widehat{x} - \widehat{x}\widehat{y},\\ \mathchoice{\frac{\mathrm{d}\widehat{y}}{\mathrm{d}\widehat{t}}}{\mathrm{d}\widehat{y}/\mathrm{d}\widehat{t}}{\mathrm{d}\widehat{y}/\mathrm{d}\widehat{t}}{\mathrm{d}\widehat{y}/\mathrm{d}\widehat{t}} &= \gamma(-\widehat{y} +\widehat{x}\widehat{y}), \end{align}\] where \[\begin{equation} \gamma = c/a. \end{equation}\] We see that nondimensionalisation has ‘tidied up’ our equation – we went from four parameters to one. Why might this matter? Well, if you fix a set of parameters you get a solution. Select another set of parameters and you get another solution. So (for given initial conditions) there are as many solutions as there are parameters. If we have decided that the four parameters are positive real numbers, then there is a four-dimensional space of solutions: quite a lot if we want to map out all the system’s behaviour! In fact, we have shown that the parameters relate to each other and that there is, in fact, only a one-dimensional space of solutions, a much simpler search.

The scalings show that the solutions just relate by constant stretching factors. For example, if I choose \(\gamma=1\), I can then set \(c=2\), thus \(a=2\). In fact, the parameters \(b\) and \(d\) are even redundant and we can just choose them (but this is rare, Lotka–Volterra is a slightly odd system). The ratio \(c/d\) just stretches any \(x\) solution (stretches its range/amplitude); \(a/b\) stretches \(y\); and the scaling along \(t\) changes the period of the solutions. That is to say, for a given \(\gamma\) we have a main solution which can be simply scaled to get other solutions, without solving the system again.

1.4.3 Solutions to Lotka–Volterra

We can solve this system using separation of variables. Dividing the two equations (and dropping hats) we obtain \[\begin{equation} \label{sepvareq} \mathchoice{\frac{\mathrm{d}y}{\mathrm{d}x}}{\mathrm{d}y/\mathrm{d}x}{\mathrm{d}y/\mathrm{d}x}{\mathrm{d}y/\mathrm{d}x}= \frac{\gamma y}{x}\left(\frac{-1+x}{1- y}\right)\implies \int \frac{1-y}{y} \,\mathrm{d}y= -\gamma \int \frac{1-x}{x} \,\mathrm{d}x. \end{equation}\] Integrating both sides of [sepvareq] we obtain \[\begin{equation} \label{dynsol} \log y-y =- \gamma(\log x -x) + C, \end{equation}\] where the constant \(C\) can be set by some initial condition, \((x(0),y(0))\). Unfortunately it is not possible to write this relationship in explicit form. This gives us the phase curves determined by the value of the constant \(C\). Parametrising this curve then gives the solutions \(x(t)\) and \(y(t)\), i.e. we could choose some behaviour for \(x(t)\) and [dynsol] will then determine the behaviour of \(y\). We will find that this kind of solution is common to such systems.

1.4.4 Equilibria of Lotka–Volterra

Looking at [reducedlotvol] we can see there are two possible equilibria where \(\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = \mathchoice{\frac{\mathrm{d}y}{\mathrm{d}t}}{\mathrm{d}y/\mathrm{d}t}{\mathrm{d}y/\mathrm{d}t}{\mathrm{d}y/\mathrm{d}t}=0\) \(\forall t\): \[\begin{equation} \label{equlibria} (x,y) = (0,0) \quad \mbox{and} \quad (x,y) = (1,1). \end{equation}\] (Dynamical systems fans would call these fixed points or steady state solutions). The \((0,0)\) solution corresponds to both populations being extinct! The second corresponds to the nonzero population densities at which the population sizes will remain fixed.

In 1.8(a), we see the varying behaviour of the closed phase curves of the system. All curves encircle the equilibrium at \((1,1)\) and as the initial conditions get closer to the equilibrium value, the radius of the curve decreases. In 1.8(b), we see the dramatic variety of morphology the parametric curves can exhibit. When the pair \((x(0),y(0))\) are initially close to the equilibrium, the curves have low-amplitude sinusoidal shape, while if \(y(0)\) is initially small, the curves have extremely sharp gradients and dramatic rates of change at the maxima.

1.4.5 Stability of Lotka–Volterra?

We have our dynamic solutions, [dynsol], and the fixed point equilibria, [equlibria]. A number of questions beg to be asked at this point:

1. Can one or both of the species die out if they are both nonzero at some time \(t\)?

2. Can an oscillating pair of populations relax to nonzero fixed values, i.e., do the populations ever settle?

An immediate observation in regards to (i) is that [dynsol] only allows \(x=0\) when \(y=0\) and vice versa, so they would have to become extinct simultaneously. The existence of periodic solutions as shown in 1.8 seems to suggest neither (i) or (ii) can occur, because the system repeats itself cyclically. A solution which decays into equilibrium would have to have a phase space diagram which spiralled inwards. In fact, we have a precise means of determining the answer to such questions which we will discuss in 2.

2.1.1 Equations of motion

Let’s consider a rigid pendulum: a bead of mass \(m\) attached to a rigid rod of length \(\ell\), depicted in 2.1(a). Three forces act on the bead: its weight, \(mg\); tension in the rod, \(T\); and friction, which opposes the motion of the bead.

The distance travelled by the bead is the arclength, \(s=\ell \theta\). The velocity and acceleration are therefore \[\begin{equation} \mathchoice{\frac{\mathrm{d}s}{\mathrm{d}t}}{\mathrm{d}s/\mathrm{d}t}{\mathrm{d}s/\mathrm{d}t}{\mathrm{d}s/\mathrm{d}t} = \ell\mathchoice{\frac{\mathrm{d}\theta}{\mathrm{d}t}}{\mathrm{d}\theta/\mathrm{d}t}{\mathrm{d}\theta/\mathrm{d}t}{\mathrm{d}\theta/\mathrm{d}t}, \quad \mathchoice{\frac{\mathrm{d}^2 s}{\mathrm{d}t^2}}{\mathrm{d}^2 s/\mathrm{d}t^2}{\mathrm{d}^2 s/\mathrm{d}t^2}{\mathrm{d}^2 s/\mathrm{d}t^2} = \ell\mathchoice{\frac{\mathrm{d}^2 \theta}{\mathrm{d}t^2}}{\mathrm{d}^2 \theta/\mathrm{d}t^2}{\mathrm{d}^2 \theta/\mathrm{d}t^2}{\mathrm{d}^2 \theta/\mathrm{d}t^2}. \end{equation}\] We say friction is proportional to velocity, \(\nu\ell \, \mathchoice{\frac{\mathrm{d}\theta}{\mathrm{d}t}}{\mathrm{d}\theta/\mathrm{d}t}{\mathrm{d}\theta/\mathrm{d}t}{\mathrm{d}\theta/\mathrm{d}t}\), and points in the direction opposing motion.

If we resolve parallel and perpendicular to the rod, given the rod makes an angle \(\theta\) to the vertical, then: \[\begin{align} \text{[parallel]} & & 0 &= T - mg\cos\theta, \\ \text{[perpendicular]} & & m \ell\mathchoice{\frac{\mathrm{d}^2 \theta}{\mathrm{d}t^2}}{\mathrm{d}^2 \theta/\mathrm{d}t^2}{\mathrm{d}^2 \theta/\mathrm{d}t^2}{\mathrm{d}^2 \theta/\mathrm{d}t^2} &= - \nu\ell \mathchoice{\frac{\mathrm{d}\theta}{\mathrm{d}t}}{\mathrm{d}\theta/\mathrm{d}t}{\mathrm{d}\theta/\mathrm{d}t}{\mathrm{d}\theta/\mathrm{d}t} - mg\sin\theta, \end{align}\] the second equation of which can be rearranged to give \[\begin{equation} \label{peneq} \mathchoice{\frac{\mathrm{d}^2 \theta}{\mathrm{d}t^2}}{\mathrm{d}^2 \theta/\mathrm{d}t^2}{\mathrm{d}^2 \theta/\mathrm{d}t^2}{\mathrm{d}^2 \theta/\mathrm{d}t^2} + \frac{\nu}{m} \mathchoice{\frac{\mathrm{d}\theta}{\mathrm{d}t}}{\mathrm{d}\theta/\mathrm{d}t}{\mathrm{d}\theta/\mathrm{d}t}{\mathrm{d}\theta/\mathrm{d}t}+ \frac{g}{\ell}\sin\theta = 0. \end{equation}\] This is a nonlinear ODE; our physical intuition tells us that the pendulum will swing with a decreasing amplitude until it relaxes to \(\theta=0\).

2.1.2 Solutions to the pendulum equation

The damped pendulum system is not an integrable system, i.e., there are no general closed-form solutions. This is generally the case for mechanical systems with friction. On the other hand, numerical solutions are simple to obtain.

An example is shown in 2.2(a): the solution is indeed oscillatory with decreasing amplitude. In (b), we see a pendulum with a very high coefficient of friction (a pendulum in fluid, say) where the damping is so strong that the angle never becomes negative: its motion is killed off on the first swing.

2.1.3 Equilibria of the pendulum equation

Accounting for periodicity, there are two equilibria to [peneq]: \[\begin{equation} \frac{g}{\ell}\sin\theta = 0 \end{equation}\] gives \(\theta(t) = 0\) and \(\theta(t) = \pi\) for all \(t\). The first \(\theta=0\) solution corresponds to a pendulum starting at the bottom of its cycle and not moving.

The second solution is far more interesting. This is when the pendulum rod is vertically upwards as in 2.1(b). The forces in the system are balanced as the rod tension balances gravity and the rotating moments are equal in either direction. But have you tried this? No matter how hard you try, it will eventually fail and the pendulum will start to rotate back to its \(\theta=0\) equilibrium. Why?

Any small variation in the pendulum – a gentle breeze or vibration in the rod – no matter how small, always grows over time. In practice, no system is perfect and such variations always exist. Mathematically we represent small variations by a linear stability analysis.

2.1.4 Order notation, \(\mathcal{O}\)

In asymptotic analysis, we say a function \(g(t,\varepsilon)\) is \(\mathcal{O}(f(t,\varepsilon))\) if \[\begin{equation} \lim_{\varepsilon\to 0} \frac{g(t,\varepsilon)}{f(t,\varepsilon)} = C \end{equation}\] with \(C\) bounded.

So, for example, let’s have \(g(t,\varepsilon) = C \varepsilon\) and \(f(t,\varepsilon) = \varepsilon\). Well, \[\begin{equation} \lim_{\varepsilon\to 0} \frac{g(t,\varepsilon)}{f(t,\varepsilon)} = C, \end{equation}\] so \(C \varepsilon\) is \(\mathcal{O}(\varepsilon)\).

But if \(g(t,\varepsilon) = C \varepsilon^2\) and \(f(t,\varepsilon) = \varepsilon\), \[\begin{equation} \lim_{\varepsilon\to 0} \frac{g(t,\varepsilon)}{f(t,\varepsilon)} = 0. \end{equation}\] So \(C \varepsilon^2\) is also \(\mathcal{O}(\varepsilon)\).

On the other hand, if \(g(t,\varepsilon) = C\) and \(f(t,\varepsilon) = \varepsilon\), then \[\begin{equation} \lim_{\varepsilon\to 0} \frac{g(t,\varepsilon)}{f(t,\varepsilon)} = \infty. \end{equation}\] So \(C\) is in some way much bigger than \(\varepsilon\) when \(\varepsilon\to 0\). We would say \(C\) is \(\mathcal{O}(1)\) but not \(\mathcal{O}(\varepsilon)\).

For non-polynomial functions the way in is to use Taylor expansion: what if \(g(t,\varepsilon) = \sin\varepsilon\) and \(f(t,\varepsilon) = \varepsilon\)? \[\begin{equation} \lim_{\varepsilon\to 0} \frac{g(t,\varepsilon)}{f(t,\varepsilon)} = \lim_{\varepsilon\to 0}\frac{\sin\varepsilon}{\varepsilon} = \lim_{\varepsilon\to 0}\frac{\varepsilon-\varepsilon^3/3!+\cdots}{\varepsilon} = \lim_{\varepsilon\to 0}\left( 1 - \frac{\varepsilon^2}{3!} + \cdots \right) = 1. \end{equation}\] So \(\sin\varepsilon\) is also \(\mathcal{O}(\varepsilon)\).

Long story short: if something is \(\mathcal{O}(\varepsilon^n)\), it is the same size as, or smaller than, \(\varepsilon^n\) when \(\varepsilon\to 0\).

2.1.5 Four steps of a linear stability analysis

The basic steps of a linear stability analysis are as follows

1. Find the system’s equilibria, \(\theta_0\) (we will commonly use the subscript \(0\) to indicate the equilibrium solution).

2. Assume a value which is changed from this equilibrium value by a very small amount, \(\theta = \theta_0 +\varepsilon\theta_1\), with \(\varepsilon\ll 1\). This is supposed to mimic the small vibration in the system. Substitute this into the equation and ignore terms of order \(\varepsilon^2\) and higher. We are left with the behaviour of the system where only small vibrations matter.

3. Solve this system to find out if our small vibration, \(\theta_1\), grows (like it would for \(\theta_0=\pi\) in the pendulum) or decays (as it would for \(\theta_0=0\)).

4. Conclude that the equilibrium is stable if we have decay (small vibrations would disappear) and unstable if they grow.

Let’s apply this to our pendulum.

2.1.6 Linear stability analysis of the damped pendulum

2.1.7 Linear stability analysis of the frictionless pendulum

If there is no friction, \(\nu =0\), steps – will be the same, except that our linearised solution is \[\begin{equation} \theta_1(t) = A\mathrm{e}^{\lambda_1 t} + B \mathrm{e}^{\lambda_2 t}, \end{equation}\] where \[\begin{equation} \lambda_{1} = \sqrt{\mp \frac{g}{\ell}}, \quad \lambda_{2} = - \sqrt{\mp \frac{g}{\ell}}. \end{equation}\]

2.1.8 Summary

1. To test the feasibility (stability) of equilibrium solutions we linearise the equation about the equilibrium state.

2. We analyse the solutions for exponential decay or growth. If there is only decay, then the solution is stable in that it is resistant to small changes. If there is any growth it is unstable as small changes destroy the solution.

3. Some special systems will have neither growth nor decay. In this case we know, from the nonlinear behaviour of the frictionless \(\theta_0=0\) pendulum equation, that this is because neighbouring solutions are periodic.

We can make these conclusions far more general…

3.1.1 Phase portraits

It’s important to keep in mind that the linear analysis about equilibria only gives a local description of the dynamics. To get a better picture of the global dynamics, we’d like to understand how the equilibria are connected, or if closed orbits exist. In one-dimensional systems, \(\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = f(x)\), this was easy as we only needed to examine the sign of \(f\) between the equilibria and the phase portrait was clear.

In two dimensions, the phase portrait will typically include:

• The equilibria and some indication of their stability,

• The trajectories near equilibria and periodic orbits that show the flow pattern.

While uniqueness of the solution dictates that trajectories cannot cross, compiling the phase portrait in 2D can be challenging. We’re going to construct the example phase portrait in 3.2. A good recipe looks like this:

1. Draw on the equilibria, coloured in according to their stability from a linear stability analysis

2. Draw the nullclines. These are the curves where either \[\begin{equation} \mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = \dot x = 0 \quad\text{or}\quad \mathchoice{\frac{\mathrm{d}y}{\mathrm{d}t}}{\mathrm{d}y/\mathrm{d}t}{\mathrm{d}y/\mathrm{d}t}{\mathrm{d}y/\mathrm{d}t} = \dot y = 0. \end{equation}\] Since the equilibria satisfy both \(\dot x=0\) and \(\dot y = 0\), the nullclines will intersect at the equilibria.

   Particularly helpful is that when trajectories intersect the nullcline, they must be either vertical (\(\dot x = 0\)) or horizontal (\(\dot y = 0\)).

3. The nullclines provide the boundary between regions of phase space where \(\dot x>0\) or \(\dot x<0\), and \(\dot y>0\) or \(\dot y<0\). Thus, the nullclines give us some sense of the direction of the trajectories in different regions of the phase plane. However, there is no way to see which sign \(\dot x,\dot y\) have in each region a priori, so we have to test at least one point in each region.

4. Infer the directions from this information. For example, in a region where \(\dot x<0\) and \(\dot y<0\), we’d expect trajectories heading down and to the left.

See how we’ve done this in 3.2 for a made-up scenario. At this point we can try to plot the full trajectories (although I haven’t done so here).

With this in mind, let’s look at Lotka–Volterra and its variants in more detail.

3.3.1 Reflections on our new model

We saw in the original Lotka–Volterra model, [reducedlotvol], that it was impossible for any population to become extinct. In the competitive Lotka–Volterra model, we introduced the notion of independent growth limitation on the individual populations to our system through the logistic model.

This has reintroduced extinction: we now have a model which allows one species to ‘win’, in addition to cases where both populations settle on fixed values. But in doing so we have found limitations. The main one is that we appear to have lost the periodic/cyclic nature of the original Lotka–Volterra model. (This is in fact true on account of the so-called Dulac criterion, but we won’t prove it in this course.)

4.1.1 Logistic growth with harvesting

We saw in 1.2 that the logistic equation has two equilibria, one which is stable and one which is unstable. Let’s now make a small modification to the logistic equation which accounts for constant-rate harvesting, \[\begin{equation} \label{logistic-harvest} \mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = ax\left(1-\frac{x}{K}\right) - H. \end{equation}\] Let’s look at the phase space in 4.1(a).

As \(H\) increases, we can see graphically we go from two equilibria, to one, to none! This corresponds to the roots of \(\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = 0\) in [logistic-harvest] becoming equal and then imaginary: the quadratic equation tells us \[\begin{equation} \label{bifurcation-equilibrium-harvesting} x_0 = \frac{K \pm \sqrt{K^2 - 4HK/a}}{2}. \end{equation}\] So if we increase \(H\), as soon as we go over \(H=aK/4\), \(\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}\) becomes always negative and we reach extinction in finite time. An interesting result about the dangerous of overharvesting! This discontinuous jump in the long-time solution is known as a catastrophe.

In general, when changes in parameter values in our models lead to qualitative changes in the equilibria and their stability, we say a bifurcation has occurred. The parameter values at which we see the bifurcation is called the bifurcation point.

To emphasise this point: bifurcation points are about parameters. In this harvesting model, we could ask this question more directly: ‘what are the equilibria, \(x_0\), for any given parameter, \(H\)?’ And we’ve actually answered this already, in [bifurcation-equilibrium-harvesting].

So why don’t we plot \(x_0\) as a function of \(H\)? That’s what we’ve done in 4.1(b): you’ll notice that we’ve plotted both of the ‘\(\pm\)’ branches of the solution. Plotting equilibria against parameters is what’s known as a bifurcation diagram and when people draw these, they note which branches refer to stable equilibria, and which refer to unstable equilibria. The stable equilibrium branch gets a solid line and the unstable equilibrium branch gets a dotted line. They can be a bit hard to get your head around to start with, so I will always draw them in purple to remind you that it’s not a phase space.

Can you spot where the bifurcation point is on the diagram? It’s where the unstable and stable lines meet. We’ll see a few more of these diagrams shortly.

4.1.2 Insect outbreak: the spruce budworm

One cannot take a mathematical biology course without encountering the spruce budworm. It shows how a simple, one-dimensional systems can yield rich dynamics with relevant predictions of a biological phenomenon. The spruce budworm, shown in 4.2, is a moth that infects spruce trees (a type of fir tree) in eastern Canada and the US. The moths produce larvae which feed on the needles of the conifers: a serious infestation can lead to complete defoliation of a forest in about four years, and so the budworm is considered one of the most destructive pests in North America.

The budworm is preyed upon by birds that for low budworm population feed only upon the budworms once the budworm population has reached a certain level. Assuming logistic growth for the budworm population in the absence of birds, the budworm population size, \(x\), is governed by the variation on the logistic equation, \[\begin{equation} \mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = rx\left(1-\frac{x}{K}\right) - \frac{Bx^2}{A^2 + x^2}. \end{equation}\] What does the new term tell us? As \(x \to \infty\), the predation rate becomes \(B > 0\), similar to our harvesting model. Meanwhile, \(A > 0\) provides a measure of the threshold population size where predation suddenly increases.

Nondimensionalising with \[\begin{equation} x = A \widehat{x}, \quad t = \frac{A}{B}\widehat{t}, \end{equation}\] and dropping hats, we get \[\begin{equation} \label{budworm-nondim} \mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = Rx\left(1-\frac{x}{k}\right) - \frac{x^2}{1+x^2}, \end{equation}\] where \(R=rA/B\) and \(k=K/A\). Once again nondimensionalisation has reduced the parameter space, from four to two!

Find equilibria

The equilibria are given by \(x=0\) and \[\begin{equation} R\left(1-\frac{x}{k}\right) = \frac{x}{1+x^2}, \label{spruce-rk} \end{equation}\] where the logistic growth is balanced by the predation. Thus, the equilibria occur when the line given by the left-hand side of [spruce-rk] intersects the curve given by the right-hand side. Naturally, we’d like to explore how these equilibria vary with \(R\) and \(k\) and, in doing so, we see the benefit of our choice of nondimensionalisation.

The curve given by \(g(x) = x/(1 + x^2)\) does not depend on either parameter, and as for the line \(h(x) = R(1-x/k)\), \(R\) is the \(y\)-intercept, while \(k\) is the where it crosses the \(x\)-axis. This is plotted in 4.3, and as you can see, this makes graphical evaluation quite easy!

Just like in the harvesting model (where we changed the parameter \(H\)), we see that changing \(R\) and \(k\) can lead to us going from three additional equilibria (we already have \(x=0\)), to two, to one: more bifurcations.

4.1.3 General bifurcations of single-species models

In general, a single-species model (that is, a single ODE), can only do so many different things as parameters are varied. In particular, the only long-time behaviour that can occur for a finite population is that it tends to a steady state value.

But what sorts of different things can happen as parameters vary? Well it turns out that in general you can have a number of different scenarios, depending on how many parameters vary, and how ‘generic’ or commonplace we expect these behaviours to be. Rather than formally and systematically listing all of these bifurcations, let’s just consider the two most common examples, and one example of a ‘complex’ bifurcation diagram formed by combining these basic bifurcations together. Importantly, these are not just examples for single-species models, but they can and do occur commonly in larger models of many species.

Saddle node bifurcation

A saddle node bifurcation (sometimes called a fold bifurcation or annihilation point) occurs when an unstable equilibrium and a stable equilibrium ‘collide’ at the same value, and then both cease to exist afterwards. This is what happened in 4.1.1 as \(H\) increased beyond \(aK/4\), our bifurcation point. Confirm you agree with this by looking at 4.1.

The model there was a little awkward to draw the bifurcation diagram for but actually the prototypical model of the saddle node bifurcation is the much simpler equation \[\begin{equation} \label{prototypical-saddle-node} \mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = p+x^2, \end{equation}\] where \(p\) is a parameter.

Draw the phase space of [prototypical-saddle-node] and check how many equilibria it has, and what their stabilities are, as \(p\) varies. Don’t worry about negative \(x\) not being feasible for representing biological populations: here we are looking at the broader mathematical principle.

By plotting the equilibria against the parameter, we create a bifurcation diagram, and you can see it in 4.5(a).

Hysteresis (again)

We can combine these two simple models to find more complex behaviour, the type of which we have already seen in the spruce budworm model.

Let’s take the pitchfork prototype and add the shift from the saddle node prototype: in other words, let’s consider the two-parameter model \[\begin{equation} \mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = p + x(q-x^2), \end{equation}\] for two parameters \(p\) and \(q\). For this discussion, we’ll fix \(q=1\) but we could do the same analysis as what follows by fixing \(p\). So we have \[\begin{equation} \label{prototypical-hysteresis} \mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = p + x(1-x^2). \end{equation}\] What does the phase space look like and how many equilibria do we have for any \(p\)?

The phase space is drawn in 4.7 and we can see that equilibria appear and disappear as we change \(p\). We can find the equilibria, \(x_0\), in terms of \(p\) by once again solving this equation when \(\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}=0\). Plotting these \(x_0\) when they are real as a function of \(p\) gives us the bifurcation diagram in 4.5(c).

Just like in the spruce budworm, we see we can vary \(p\) such that returning \(p\) to its initial value doesn’t return the system to its original state: hysteresis. If you don’t see why, read again the discussion under the spruce budworm on page and see that it applies exactly to this model.

These three bifurcation phenomena – saddle nodes, pitchforks and hysteresis – will reappear next term: for now, we’re learning about them because with these bifurcation diagrams you can read off essentially the entire system’s behaviour for any given parameter. Pretty cool!

5.1.1 Discrete random walks

As a concrete example, let us consider an individual on a line initially at a position \(x=0\). The individual moves left or right in integer steps with a probability of moving right \(p\) and hence left of \((1-p)\). After \(n\) steps a path can be encoded as LRLLRRR…. If we repeat this process a large number of times we get a set of paths – a set of random walks. In 5.1, we show five random walks of length \(100\) for \(p=0.5\), \(0.6\) and \(0.9\) respectively. Here you see the individual paths; in 5.2 we see the same for 200 walks! For \(p=0.5\), the end positions are reasonably evenly spread about \(x=0\). For \(p=0.6\), the story is similar except that the spread is about the mean positions, \(x=2np-n\), and the spread is not equal either side of this mean. This nonzero mean, indicated by the dotted line, can be interpreted as a natural drift of the set of individuals all starting at the same point due to the probability bias.

5.1.2 Continuum limit of random walks

In our course, we will focus on the continuum limit of such a model. Rather than making the steps of size \(1\), we specify them to be \(\mathrm{d}x\) which will be vanishingly small, and each step is taken in a time \(\mathrm{d}t\). We let \(c(x,t)\) be the continuous probability that, at a time \(t\), a particle (population member) reaches a displacement \(x\) at a time \(t\). This implies that at a time \(t-\mathrm{d}t\), the particle must have been at either \(x-\mathrm{d}x\) or \(x+\mathrm{d}x\). If \(p\) is the probability that the particle moves right, we therefore have \[\begin{equation} c(x,t)= p c(x-\mathrm{d}x,t-\mathrm{d}t) + (1-p) c(x+\mathrm{d}x,t-\mathrm{d}t). \end{equation}\] If we Taylor expand this to \(\mathcal{O}(\varepsilon^2)\), where \(\varepsilon= \max(\mathrm{d}t,\mathrm{d}x)\), we obtain \[\begin{align} c(x,t) &\approx p\left[c(x,t) - \mathrm{d}x \mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} - \mathrm{d}t \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} + \frac{1}{2}\left( \mathrm{d}x^2 \mathchoice{\frac{\partial^2 c}{\partial x^2}}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2} + 2 \, \mathrm{d}x \, \mathrm{d}t \mathchoice{\frac{\partial^2 c}{\partial x \partial t}}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t} + \mathrm{d}t^2 \mathchoice{\frac{\partial^2 c}{\partial t^2}}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2} \right)\right] \\ & \quad + (1-p)\left[c(x,t) + \mathrm{d}x \mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} - \mathrm{d}t \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} + \frac{1}{2}\left( \mathrm{d}x^2 \mathchoice{\frac{\partial^2 c}{\partial x^2}}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2} - 2 \, \mathrm{d}x \, \mathrm{d}t \mathchoice{\frac{\partial^2 c}{\partial x \partial t}}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t} + \mathrm{d}t^2 \mathchoice{\frac{\partial^2 c}{\partial t^2}}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2} \right)\right] \nonumber\\ &= c(x,t) + (1-2p)\mathrm{d}x \mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} - \mathrm{d}t \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} + \frac{1}{2}\mathrm{d}x^2 \mathchoice{\frac{\partial^2 c}{\partial x^2}}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2} - (1-2p)\, \mathrm{d}x \, \mathrm{d}t \mathchoice{\frac{\partial^2 c}{\partial x \partial t}}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t} + \frac12 \mathrm{d}t^2 \mathchoice{\frac{\partial^2 c}{\partial t^2}}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2}, \end{align}\] so that to order \(\mathcal{O}(\varepsilon^2)\), \[\begin{equation} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} = (1-2 p)\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}\mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} + \frac12 \frac{(\mathrm{d}x)^2}{\mathrm{d}{t}} \mathchoice{\frac{\partial^2 c}{\partial x^2}}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2} - (1-2p)\,\mathrm{d}x\mathchoice{\frac{\partial^2 c}{\partial x \partial t}}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t} + \frac12 \mathrm{d}t \mathchoice{\frac{\partial^2 c}{\partial t^2}}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2}. \end{equation}\] Taking the limit \(\mathrm{d}x\to 0\) and \(\mathrm{d}t\to 0\), the final two terms will vanish. The coefficient \((1-2p)\,\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}\) has the dimensions of a velocity, so we set \((2p-1)\,\mathchoice{\frac{\mathrm{d}x}{\mathrm{d}t}}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t}{\mathrm{d}x/\mathrm{d}t} = v\) (with the flipped sign so that it points in the positive direction for high \(p\) biases) and we define a constant \(D\), \[\begin{equation} D = \frac12 \frac{(\mathrm{d}x)^2}{\mathrm{d}t}, \end{equation}\] which we call the diffusion constant for reasons which will follow. One should note that the ratio \(v/D\) will diverge as \(\mathrm{d}x \to 0\) unless \(p\) is very close to \(1/2\), so as the step size gets smaller, the velocity term (the so-called advective term) dominates. Then we have \[\begin{equation} \label{advectiondiffusion} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} = -v\mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} + D \mathchoice{\frac{\partial^2 c}{\partial x^2}}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}, \end{equation}\] a partial differential equation determining the behaviour of \(c(x,t)\). Note that if \(p=1/2\), there is no velocity: we shall come to associate this term with a drift. The physical description of the term, \(D\), can be seen by deriving [advectiondiffusion] in a continuum setting, as we’ll do now.

5.4.1 The Fourier transform

The Fourier transform of a function \(f(x)\), denoted \(\mathcal{F}[f(x)](k)\), is given by \[\begin{equation} \mathcal{F}[f(x)](k) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\mathrm{e}^{-\mathrm{i}k x}f(x)\,\mathrm{d}x. \end{equation}\] Its inverse, \(\mathcal{F}^{-1}[g(k)](x)\), is \[\begin{equation} \mathcal{F}^{-1}[g(k)](x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\mathrm{e}^{\mathrm{i}k x}g(k)\,\mathrm{d}{k}. \end{equation}\] You might have seen a slightly different version of these where only one of the terms has a \(1/2\pi\) factor, but either way, crucially it can be shown that \[\begin{equation} \mathcal{F}^{-1}[\mathcal{F}[f(x))]] = f(x). \end{equation}\] You’re not expected to derive this result, merely to be aware of it. The critical property of the Fourier transform which makes it of use in solving PDEs is its effect on derivatives. The following result can be fairly easily demonstrated: \[\begin{equation} \label{fderiv} \mathcal{F}\left[\mathchoice{\frac{\partial^{n} f}{\partial x^{n}}}{\partial^{n} f/\partial x^{n}}{\partial^{n} f/\partial x^{n}}{\partial^{n} f/\partial x^{n}}\right](k) = (\mathrm{i}k)^{n}\mathcal{F}[f(x)](k). \end{equation}\] Again, you don’t need to derive this result, you just need to know how to apply it.

5.4.2 Solving the diffusion equation

Applying [fderiv] to the first term in [diffprob] we obtain \[\begin{equation} \mathchoice{\frac{\partial\mathcal{F}[u](k)}{\partial t}}{\partial\mathcal{F}[u](k)/\partial t}{\partial\mathcal{F}[u](k)/\partial t}{\partial\mathcal{F}[u](k)/\partial t} = -Dk^2\mathcal{F}[u](k). \end{equation}\] Just like that, we have transformed a second order PDE into a first order ODE. Its solution is clearly an exponential, \[\begin{equation} \mathcal{F}[u(x,t)](k) = \mathcal{F}[u(x,0)](k)\mathrm{e}^{-D k^2 t}. \end{equation}\]

To get the real space solution, we must apply the inverse operator. This is the price we pay for using the transform to simplify the system. This is not always straightforward, but in this case this is not so problematic. First we must find the initial condition \(\mathcal{F}[c(x,0)](k)\), which, from [deltafunc], is the Fourier transform of the Dirac delta function, \[\begin{equation} \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\delta(x)\mathrm{e}^{-\mathrm{i}kx}\,\mathrm{d}x= \frac{1}{\sqrt{2\pi}}. \end{equation}\] Thus our final solution is \[\begin{equation} u(x,t) = \mathcal{F}^{-1}[\mathcal{F}[u(x,t)]] = \frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{e}^{\mathrm{i}kx}\mathrm{e}^{-D k^2 t}\,\mathrm{d}{k}. \end{equation}\] In order to perform this integral we note the following known result: \[\begin{equation} \mathcal{F}[\mathrm{e}^{-\alpha x^2}](k) = \frac{1}{\sqrt{2\alpha}}\mathrm{e}^{-k^2/4\alpha}. \end{equation}\] We therefore have \[\begin{equation} \label{fundamentalsol} u(x,t) = \frac{1}{2 \sqrt{\pi Dt}} \mathrm{e}^{-x^2/4Dt}. \end{equation}\]

This tells us the density spreads out over time following a Gaussian distribution, as depicted in 5.6.

As we remarked at the start of the chapter, the limit of the random walk as the step size reduces to zero is a Gaussian distribution. We then showed, using a slightly hand-wavy argument, that a limit of this distribution was the advection–diffusion equation. We have now squared the circle so to speak, in the sense that we see that the diffusion equation gives Gaussian-like spreading of a point source. In fact it was Einstein (and Marian Smoluchowski) who first made this link precise at the turn of the previous century.

This solution, [fundamentalsol], is known as the fundamental solution of the diffusion equation. Why is it so fundamental? Because it allows us to answer the question ‘what about different initial conditions?’.

5.4.3 Fundamental solutions (or Green’s functions)

Why is a delta function which spreads out over time called the fundamental solution? In general we might have some linear differential operator \(\mathcal{L}\). When acting on a function \(u(\mathbfit{x},t)\) we want to solve the following initial value problem on a domain \(V\) with boundary \(\partial V\): \[\begin{equation} \label{fullprob} \mathcal{L}u(\mathbfit{x},t) = 0,\quad u(\mathbfit{x},0) = g(\mathbfit{x}),\quad u(\partial{V},t) = 0, \end{equation}\] that is to say, \(g(\mathbfit{x})\) is our initial condition.

For example, the diffusion equation operator is \(\mathcal{L}= \mathchoice{\frac{\partial}{\partial t}}{\partial/\partial t}{\partial/\partial t}{\partial/\partial t}- D\nabla^2\), so that \[\begin{equation} \mathcal{L}u(\mathbfit{x},t) = \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} -D\nabla^2 u = 0. \end{equation}\] \(\mathcal{L}\) is linear so for \(a,b\) scalar constants and \(u,v\) scalar densities, \[\begin{equation} \mathcal{L}(au+bv) = a\mathcal{L}u + b\mathcal{L}v. \end{equation}\] The fundamental solution or Green’s function of \(\mathcal{L}\) can be defined as the solution, \(u_f\), to \[\begin{equation} \label{fundprob} \mathcal{L}u_f(\mathbfit{x},t) = 0, \quad u_f(\mathbfit{x},0) = \delta(\mathbfit{x}), \quad u_f(\partial{V},t) = 0, \end{equation}\] where specifically the initial condition is the delta function. See how our setup for the diffusion equation, [diffprob], is in this form.

We can relate the two initial conditions by the delta function identity, \[\begin{equation} \label{deltaid} \int_{V}\delta(\mathbfit{x}- \mathbfit{s})g(\mathbfit{s})\,\mathrm{d}\mathbfit{s} = g(\mathbfit{x}). \end{equation}\] This combination of the two functions \(\delta\) and \(g\) on the left-hand is known as a convolution. We now propose that our solution \(u\) can be written as another convolution: \[\begin{equation} \label{gensol} u(\mathbfit{x},t) = \int_{V}u_f(\mathbfit{x}-\mathbfit{s},t) g(\mathbfit{s})\,\mathrm{d}\mathbfit{s}. \end{equation}\] To see this satisfies \(\mathcal{L}u = 0\), we note that, since the linear operator acts on \((\mathbfit{x},t)\) but the integral is over \(\mathbfit{s}\), we can take it inside, \[\begin{equation} \label{lingreen} \mathcal{L}u(\mathbfit{x},t) = \int_{V}g(\mathbfit{s})\underbrace{\mathcal{L}u_f(\mathbfit{x}-\mathbfit{s},t)}_{=\,0}\,\mathrm{d}\mathbfit{s} = 0. \end{equation}\] In addition, as \(\lim_{t \to 0} u_f = \delta(\mathbfit{x})\), we will obtain the right initial condition through [deltaid]. Thus, solving [fullprob] can be reduced to solving the fundamental problem, [fundprob]: generally an easier task! Integrating [gensol] then leaves us with the full solution. This can often be of benefit for analytic solutions and asymptotic approximations.

5.4.4 Solving the diffusion equation with different initial conditions

Let’s say we wish to solve the following problem: \[\begin{equation} \label{rectheat} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} = D \mathchoice{\frac{\partial^2 c}{\partial x^2}}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2},\quad c(x,0) = \left\{\begin{array}{ll} 0 & x< -a\\ 1 & -a \leq x \leq a\\ 0 & x>a \end{array} \right. , \quad c(\pm \infty ,t) = 0. \end{equation}\] The initial condition is a ‘population’ spread evenly over a domain \(x\in[-a,a]\). Using the fundamental solution in [fundamentalsol], and [gensol], the solution to the problem in [rectheat] is \[\begin{equation} c(x,t) = \frac{1}{ 2\sqrt{\pi Dt}}\int_{-a}^{a} \mathrm{e}^{-(x-s)^2/4Dt}\,\mathrm{d}{s} = \frac{1}{2}\left[\operatorname{erf}\left(\frac{a-x}{\sqrt{4Dt}}\right) + \operatorname{erf}\left(\frac{a+x}{\sqrt{4Dt}}\right)\right], \end{equation}\] where \[\begin{equation} \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}}\int_{0}^{x}\mathrm{e}^{-s^2}\mathrm{d}{s}. \end{equation}\] is the error function, a function with well-known properties. The solution is shown to vary with \(t\) in 5.7.

5.5.1 The porous medium equation

Rather than carrying on with the standard one-dimensional case we now turn to a radially symmetric two-dimensional case. We will again ignore source terms (\(f\equiv 0\)) and we shall assume no velocity \(\mathbfit{v}\), so only diffusive behaviour is present. We seek a fundamental solution with particle density \(Q\), i.e., we are seeking a solution to the problem \[\begin{equation} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} = D_0\boldsymbol{\nabla}\cdot\left[\left(\frac{c}{c_0}\right)^m \boldsymbol{\nabla}c\right],\quad c(\mathbfit{x},0) = Q\delta(\mathbfit{x}), \quad c(\mathbfit{x} \to \infty,t)=0. \label{radial-diff-eqn} \end{equation}\] By \(c(\mathbfit{x} \to \infty,t)\) we mean to imply that the distribution must eventually go to zero in all directions.

In polar coordinates, \((r,\theta)\), you will remember that \[\begin{equation} \boldsymbol{\nabla}c = \mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r}\widehat{\mathbfit{e}}_r + \frac{1}{r}\mathchoice{\frac{\partial c}{\partial\theta}}{\partial c/\partial\theta}{\partial c/\partial\theta}{\partial c/\partial\theta}\widehat{\mathbfit{e}}_\theta, \end{equation}\] and for some vector \(\mathbfit{A} = A_r\widehat{\mathbfit{e}}_r + A_{\theta}\widehat{\mathbfit{e}}_\theta\), \[\begin{equation} \boldsymbol{\nabla}\cdot\mathbfit{A} = \frac{1}{r}\left(\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}(rA_r) + \mathchoice{\frac{\partial A_{\theta}}{\partial\theta}}{\partial A_{\theta}/\partial\theta}{\partial A_{\theta}/\partial\theta}{\partial A_{\theta}/\partial\theta}\right). \end{equation}\] If we assume radial symmetry, \(c(r,\theta,t) \equiv c(r,t)\), then the right-hand side of [radial-diff-eqn] yields \[\begin{equation} \label{porusradial} D_0\boldsymbol{\nabla}\cdot\left[\left(\frac{c}{c_0}\right)^m\mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r} \widehat{\mathbfit{e}}_r\right] = \frac{D_0}{r}\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}\left[r\left(\frac{c}{c_0}\right)^m\mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r} \right], \end{equation}\] and our equation becomes \[\begin{equation} \label{porous-media-eqn} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} = \frac{D_0}{r}\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}\left[r\left(\frac{c}{c_0}\right)^m\mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r} \right]. \end{equation}\] This is an example of the porous media equation. We will shortly show that the solution is \[\begin{equation} \label{porussol} c(r,t) = \left\{ \begin{array}{ll} \displaystyle\frac{c_0}{\lambda^2(t)}\left[1 - \left(\frac{r}{r_0 \lambda(t)}\right)^2\right]^{1/m} & \mbox{} r\leq r_0\lambda(t) , \\[6pt] 0 & \mbox{} r > r_0 \lambda(t) , \end{array} \right. \end{equation}\] where \[\begin{equation} \lambda(t) = \left(\frac{t}{t_0}\right)^{\frac{1}{2(m+1)}}, \quad t_0 = \frac{r_0^2 m}{4D_0(m+1)}, \quad r_0^2=\frac{Q}{\pi c_0}\left(1+\frac{1}{m}\right). \end{equation}\] The solution has compact support – it is zero outside of a finite radius – and this finite radius grows in time. The volume, \(Q\), is conserved and we see in 5.8 the solution spreading out over time.

5.5.2 Similarity solutions

Deriving this solution involves using another technique for solving partial differential equation: a symmetry reduction technique known as looking for similarity solutions. The basic idea is that if we can find some transformation of our solution which changes its shape but still solves the basic equations, then we should be able to use this transformation to eliminate one of the independent variables. This is a bit like with Fourier transforms when this elimination turns a PDE into an ODE. Unfortunately, of course, first you must such a scaling, and this is not always easy! But for a significant class of PDEs with one spatial dimension the generic transformation of \(u(x,t)\) given by \[\begin{equation} \label{dilation} U = \lambda^{-\alpha} u, \quad X = \lambda^\beta x, \quad T = \lambda t %u(x,t) = \lambda^{\alpha} v(\lambda^{\beta}x,\lambda t) \end{equation}\] for some suitable constants \(\alpha\) and \(\beta\) will give the original PDE for \(U(X,T)\) (a dilation).

For example, consider the 1D diffusion equation from the previous chapter, \[\begin{equation} \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} = D\mathchoice{\frac{\partial^2 u}{\partial x^2}}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}. \end{equation}\] Applying [dilation] we get \[\begin{equation} \lambda^{\alpha+1}\mathchoice{\frac{\partial U}{\partial T}}{\partial U/\partial T}{\partial U/\partial T}{\partial U/\partial T} = D\lambda^{\alpha+2\beta}\mathchoice{\frac{\partial^2 U}{\partial X^2}}{\partial^2 U/\partial X^2}{\partial^2 U/\partial X^2}{\partial^2 U/\partial X^2}. \end{equation}\] We can spot that if \(\beta=1/2\) then \[\begin{equation} \mathchoice{\frac{\partial U}{\partial T}}{\partial U/\partial T}{\partial U/\partial T}{\partial U/\partial T} = D \mathchoice{\frac{\partial^2 U}{\partial X^2}}{\partial^2 U/\partial X^2}{\partial^2 U/\partial X^2}{\partial^2 U/\partial X^2}, \end{equation}\] i.e., \(U\) satisfies the diffusion equation! So we see a transformation in the form [dilation] can lead to the same differential equation, just in \(U\) instead of \(u\).

Note that \(\alpha\) turned out to be arbitrary here: as we will see, we can pick it to suit the initial conditions.

But why do we care that under this transformation, \(U\) satisfies the same PDE? Well, note that the groupings \[\begin{equation} \eta = \frac{X}{T^\beta} = \frac{x}{t^\beta} \quad \text{and} \quad \xi = UT^\alpha = ut^{\alpha} \label{sim-sol-sub} \end{equation}\] are invariant (i.e. unchanged) under the transformation. By a theorem sometimes known as Morgan’s theorem,³ we can write one invariant as a function of the other, \[\begin{equation} \xi = v(\eta), \end{equation}\] for some function \(v\). Substituting in the definitions of \(\eta\) and \(\xi\), we get \[\begin{equation} u(x,t) = \frac{1}{t^\alpha} v \left( \frac{x}{t^\beta} \right), \end{equation}\] and this is the key result! This implies that the solution \(u(x,t)\) has the same shape as the ‘master shape’ \(v(\eta)\), just stretched by a factor \(t^{\beta}\) along the \(x\)-axis and squashed by a factor \(t^{\alpha}\) along the \(y\)-axis. This shape can then be appropriately scaled to give the time-varying behaviour of the system (determined by the parameters \(\alpha\) and \(\beta\), depending on the actual equation).

We have seen this behaviour already. Recall the fundamental solution of the diffusion equation, [fundamentalsol] (with \(u\) taking the \(c\) role), \[\begin{equation} u(x,t) =\frac{1}{2 \sqrt{\pi Dt}} \mathrm{e}^{-x^2/4Dt}. \end{equation}\] We have seen above that for the diffusion equation, \(\beta=1/2\), so if we write \(\eta = x/t^{1/2}\) and \(v = ut^{\alpha}\) (from [sim-sol-sub]) then \[\begin{equation} v(\eta,t) =\frac{t^{\alpha}}{2\sqrt{\pi D t}} \mathrm{e}^{-\eta^2/4D}. \end{equation}\] We said earlier that \(\alpha\) can be set to fit the initial conditions. Well, the fundamental solution has a very specific initial condition which gives form to this solution, and here if we set \(\alpha=1/2\) then \[\begin{equation} v(\eta) =\frac{1}{2 \sqrt{\pi D}}\mathrm{e}^{-\eta^2/4D}. \end{equation}\]

So we see our solution \(u(x,t)\) can be written in the form \[\begin{equation} u(x,t) = \frac{1}{t^\alpha}v\left(\frac{x}{t^{\beta}}\right) = \frac{1}{t^{1/2}}v\left(\frac{x}{t^{1/2}}\right), \end{equation}\] and the solution will still satisfy the diffusion equation.

As we saw in 5.6, this confirms that the solution \(u(x,t)\) has the same shape as the master shape \(v(\eta)\), stretched by a factor \(t^{1/2}\) along the \(x\)-axis and squashed by a factor \(t^{1/2}\) along the \(y\)-axis.

5.5.3 Back to the porous media equation

Let’s go back to [porous-media-eqn], \[\begin{equation} \label{radialpde} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} = \frac{D_0}{r}\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}\left[r \left(\frac{c}{c_0}\right)^m \mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r}\right], \end{equation}\] and, using the observed similarity solution, make the substitution \[\begin{equation} c(r,t) = \frac{1}{t^{\alpha}}v\left(\frac{r}{t^{\beta}}\right). \end{equation}\] Denoting \(\eta=r/t^{\beta}\), so that \(v=v(\eta)\), we first change the left-hand side: \[\begin{align} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} &= -\alpha\frac{1}{t^{\alpha+1}}v + \frac{1}{t^{\alpha}}\mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}\mathchoice{\frac{\partial\eta}{\partial t}}{\partial\eta/\partial t}{\partial\eta/\partial t}{\partial\eta/\partial t} \\ &= -\alpha\frac{1}{t^{\alpha+1}}v - \frac{r \beta}{t^{\alpha}t^{\beta+1}}\mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}\\ \label{fincaldcdt}& = -\alpha\frac{1}{t^{\alpha+1}}v - \frac{\eta \beta}{t^{\alpha+1}}\mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}. \end{align}\] On the right-hand side, \[\begin{equation} \mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r} = \frac{1}{t^\beta}\mathchoice{\frac{\partial}{\partial\eta}}{\partial/\partial\eta}{\partial/\partial\eta}{\partial/\partial\eta} \implies \mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r}= \frac{1}{t^{\alpha+\beta}}\mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}, \end{equation}\] and so we have \[\begin{align} \frac{D_0}{r}\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}\left[r \left(\frac{c}{c_0}\right)^m \mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r}\right]&=\frac{D_0}{\eta t^{2\beta}}\mathchoice{\frac{\partial}{\partial\eta}}{\partial/\partial\eta}{\partial/\partial\eta}{\partial/\partial\eta}\left[\frac{t^\beta \eta}{t^{\alpha m}}\left(\frac{v}{c_0}\right)^m \frac{1}{t^{\alpha+\beta}}\mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}\right]\\ \label{finalrhs}&= \frac{1}{t^{\alpha(m+1)+2\beta}}\frac{D_0}{\eta}\mathchoice{\frac{\mathrm{d}}{\mathrm{d}\eta}}{\mathrm{d}/\mathrm{d}\eta}{\mathrm{d}/\mathrm{d}\eta}{\mathrm{d}/\mathrm{d}\eta}\left[\eta\left(\frac{v}{c_0}\right)^m \mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}\right]. \end{align}\] Comparing [fincaldcdt] and [finalrhs] we see we want to set the \(t\)-exponents equal to each other, \[\begin{equation} \alpha + 1 = \alpha(m+1) + 2\beta \implies 1 = \alpha m + 2\beta. \end{equation}\] We can multiply through by \(t^{\alpha+1}\) to obtain \[\begin{equation} \alpha v + \beta \mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}\eta +\frac{D_0}{\eta}\mathchoice{\frac{\mathrm{d}}{\mathrm{d}\eta}}{\mathrm{d}/\mathrm{d}\eta}{\mathrm{d}/\mathrm{d}\eta}{\mathrm{d}/\mathrm{d}\eta}\left[\eta\left(\frac{v}{c_0}\right)^m\mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}\right] =0. \label{rel-between-alpha-and-beta} \end{equation}\] So now we have a relationship between \(\alpha\) and \(\beta\), but still some freedom in choosing one of them.

We now manipulate [rel-between-alpha-and-beta] to simplify it with an appropriate choice of \(\alpha\). In particular, we want to write our equation in the form \(\mathchoice{\frac{\mathrm{d}[\cdots]}{\mathrm{d}\eta}}{\mathrm{d}[\cdots]/\mathrm{d}\eta}{\mathrm{d}[\cdots]/\mathrm{d}\eta}{\mathrm{d}[\cdots]/\mathrm{d}\eta}\) so that we can then integrate and create a first-order equation.

Multiply our equation through by \(\eta\), \[\begin{equation} \label{newode} \alpha v \eta + \beta \mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta} \eta^2 + D_0\mathchoice{\frac{\mathrm{d}}{\mathrm{d}\eta}}{\mathrm{d}/\mathrm{d}\eta}{\mathrm{d}/\mathrm{d}\eta}{\mathrm{d}/\mathrm{d}\eta}\left[\eta\left(\frac{v}{c_0}\right)^m \mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}\right]=0. \end{equation}\] The idea is to turn the first two terms into a derivative so we can integrate. You might notice that for some \(\gamma\), \[\begin{equation} \gamma\mathchoice{\frac{\mathrm{d}}{\mathrm{d}\eta}}{\mathrm{d}/\mathrm{d}\eta}{\mathrm{d}/\mathrm{d}\eta}{\mathrm{d}/\mathrm{d}\eta}[v\eta^2] = \gamma \mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta} \eta^2 + 2\gamma v\eta, \end{equation}\] which is almost the form we want. If we set \(\gamma=\beta\), then we have two simultaneous equations, \[\begin{align} 1 & = \alpha m + 2 \beta,\\ \alpha & = 2\beta. \end{align}\] This gives us \(\alpha = 1/(m+1)\) and \(\beta=1/[2(m+1)]\), so [newode] becomes \[\begin{equation} \beta\mathchoice{\frac{\mathrm{d}}{\mathrm{d}\eta}}{\mathrm{d}/\mathrm{d}\eta}{\mathrm{d}/\mathrm{d}\eta}{\mathrm{d}/\mathrm{d}\eta}[v\eta^2] + D_0\mathchoice{\frac{\mathrm{d}}{\mathrm{d}\eta}}{\mathrm{d}/\mathrm{d}\eta}{\mathrm{d}/\mathrm{d}\eta}{\mathrm{d}/\mathrm{d}\eta}\left[\eta\left(\frac{v}{c_0}\right)^m \mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}\right] = 0. \end{equation}\] Integrating with respect to \(\eta\) gives \[\begin{equation} \beta v\eta^2 + D_0\eta\left(\frac{v}{c_0}\right)^m \mathchoice{\frac{\mathrm{d}v}{\mathrm{d}\eta}}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta}{\mathrm{d}v/\mathrm{d}\eta} = C, \end{equation}\] with \(C\) a constant of integration.

We choose the imposed boundary condition that \(v(\eta) \to 0\) sufficiently fast as \(\eta(r) \to \infty\) which implies the constant of integration is \(C=0\). For any finite \(t\), this is equivalent to saying the population vanishes as \(r\to \infty\).

We solve the differential equation by noting either \(v=0\) or, using separation of variables, \[\begin{align} \beta \int \eta\,\mathrm{d}{\eta} &= - \int \frac{D_0}{c_0^m}v^{m-1}\mathrm{d}{v},\\ \implies \frac{\beta \eta^2}{2} &=- \frac{D_0}{m c_0^m}v^{m} + A,\\ \implies v^m &= \frac{m c_0^m}{D_0}\left(A - \frac{\beta\eta^2}{2}\right). \end{align}\] Finally, using \(c= v/t^\alpha\) and \(\eta= r/t^{\beta}\), we get \[\begin{equation} c(r,t)= \frac{m^{1/m} c_0}{t^{\alpha}D_0^{1/m}}\left(A - \frac{\beta r^2}{2t^{2\beta}}\right)^{1/m}, \end{equation}\] where \(\alpha = 1/(m+1)\) and \(\beta=1/[2(m+1)]\). The concentration, \(c\), is only non-negative (and real) if the term is the bracket is non-negative, i.e. \[\begin{equation} r^2 \leq \frac{2A}{\beta}t^{2\beta}, \end{equation}\] and so outside this region, the solution must be the \(c=0\) solution we ignored earlier. In other words, this solution has compact support, and we can see the size of the support grows with time. Thinking back, this is what we might expect from an equation that represents a type of diffusion.

To get \(A\) and hence the final solution, [porussol], we enforce that the total density (mass) is conserved, \[\begin{equation} \int c \, \mathrm{d}S = 2 \pi \int_0^{\sqrt{\frac{2A}{\beta}}t^\beta} c r \, \mathrm{d}r = Q. \end{equation}\] This is a little tedious and I don’t really want you to spend your time doing this. We’ve shown that this solution has compact support, which is the important result. And the main method we’re interested in here is the method of similarity solutions.