Reaction-Diffusion Systems

This is an incomplete page

We present here a brief theory of reaction diffusion waves.

Law of Mass Action

The law of mass action states that equation rates are proportional to the concentration of reacting species and the ratio in which they combined. It is discussed in detail in Billingham and King 2000. We will present here a few simple examples.

Example 1: Simple Decay

Suppose we have of chemical [math]\displaystyle{ P }[/math] which decays to [math]\displaystyle{ A }[/math], i.e.

[math]\displaystyle{ P \to A }[/math]

with rate [math]\displaystyle{ k[P] }[/math] where [math]\displaystyle{ [P] }[/math] denotes concentration. Then if we set [math]\displaystyle{ p=[P] }[/math] and [math]\displaystyle{ a = [A] }[/math] we obtain the equations

[math]\displaystyle{ \frac{\mathrm{d}p}{\mathrm{d}t} = -kp\,\,\,\textrm{and}\,\,\, \frac{\mathrm{d}a}{\mathrm{d}t} = kp }[/math]

which has solution

[math]\displaystyle{ p = p_0 e^{-kt}\,\,\,\textrm{and}\,\,\, a = a_0 + p_0(1-e^{-kt}) }[/math]

Example 2: Quadratic Autocatalysis

This example will be important when we consider reaction diffusion problems. We consider the reaction

[math]\displaystyle{ A + B \to 2B }[/math]

with rate proportional to [math]\displaystyle{ k[A][B] }[/math]. If we define [math]\displaystyle{ a = [A] }[/math] and [math]\displaystyle{ b = [B] }[/math] we obtain the following equations

[math]\displaystyle{ \frac{\mathrm{d}a}{\mathrm{d}t} = -kab\,\,\,\textrm{and}\,\,\, \frac{\mathrm{d}b}{\mathrm{d}t} = kab }[/math]

We can solve these equations by observing that

[math]\displaystyle{ \frac{\mathrm{d}(a+b)}{\mathrm{d}t} = 0 }[/math]

so that [math]\displaystyle{ a = b = a_0 + b_0 }[/math]. We can then eliminate [math]\displaystyle{ a }[/math] to obtain

[math]\displaystyle{ \frac{\mathrm{d}b}{\mathrm{d}t} = k(a_0 + b_0 - b)b }[/math]

[math]\displaystyle{ \lt math\gt Insert formula here }[/math]</math>

which is separable with solution

[math]\displaystyle{ b = \frac{b_0(a_0 + b_0)e^{k(a_0 + b_0)t}}{a_0 + b_0e^{k(a_0 + b_0)t}} }[/math]

and

[math]\displaystyle{ a = \frac{a_0(a_0 + b_0)}{a_0 + b_0e^{k(a_0 + b_0)t}} }[/math]

Note that [math]\displaystyle{ a\to 0 }[/math] and [math]\displaystyle{ b\to a_0 + b_0 }[/math] as [math]\displaystyle{ t\to \infty. }[/math]

Diffusion

The equation for spatially homogeneous diffusion of a chemical with concentration [math]\displaystyle{ c }[/math] is

[math]\displaystyle{ \partial_t c = D\nabla^2 c }[/math]

which is the heat equation. We will consider this in only one spatial dimension. Consider it on the boundary [math]\displaystyle{ -\infty \lt x \lt \infty }[/math]. In this case we can solve by the Fourier transform and obtain

[math]\displaystyle{ \partial_t \hat{c} = -D k^2 \hat{c} }[/math]

where [math]\displaystyle{ \hat{c} }[/math] is the Fourier transform of [math]\displaystyle{ c }[/math]. This has solution

[math]\displaystyle{ \hat{c} = \hat{c}_0 e^{-D k^2 t} }[/math]

We can find the inverse transform using convolution and obtain

[math]\displaystyle{ c(x,t) = \frac{1}{\sqrt{4\pi D t}} \int_{-\infty}^{\infty} c_0(x) e^{(x-s)^2/4Dt}ds }[/math]

Solution of the dispersion equation using FFT

We can solve the dispersion equation using the discrete Fourier transform and its closely related numerical implementation the FFT. We consider the concentration on the finite domain [math]\displaystyle{ -L \leq x \leq L }[/math] and use a Fourier series expansion

[math]\displaystyle{ c(x,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(t) e^{\mathrm{i} k_n x} }[/math]

where [math]\displaystyle{ k_n = \pi n /L }[/math]. If we substitute this into the diffusion equation we obtain

[math]\displaystyle{ c(x,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(0)e^{-k_n^2 D t} e^{\mathrm{i} k_n x} }[/math]

Note that this is not the same solution as we obtained on the infinite domain because of the boundary conditions on the finite domain. The coefficients [math]\displaystyle{ \hat{c}_n(0) }[/math] are found using the initial conditions so that

[math]\displaystyle{ \hat{c}_n(0) = \frac{1}{2L} \int_{-L}^{L} e^{-\mathrm{i} k_n x} c_0(x) dx }[/math]

The key to the numerical solution of this equation is the use of the FFT. We begin by discretising the domain into a series of [math]\displaystyle{ N }[/math] points [math]\displaystyle{ x_m = -L + 2Lm/N }[/math]. We then use this to approximate the integral above and obtain

[math]\displaystyle{ \hat{c}_n(0) = \frac{1}{N} \sum_{m=0}^{N-1} e^{-\mathrm{i} k_n x_m} c_0(x_m) }[/math]

[math]\displaystyle{ = \frac{1}{N} \sum_{m=0}^{N-1} e^{-2\mathrm{i} \pi nm/N} e^{-\mathrm{i} \pi n} c_0(x_m) }[/math]

We also get

[math]\displaystyle{ c(x_m,t) = \sum_{n=-\infty}^{\infty} \hat{c}_n(0) e^{-k_n^2 D t} e^{\mathrm{i} k_n x_m} }[/math]

[math]\displaystyle{ = \sum_{n=-\infty}^{\infty} \hat{c}_n(0) e^{-k_n^2 D t} e^{2\mathrm{i} \pi nm/N} e^{\mathrm{i} \pi n} }[/math]

but we know that

[math]\displaystyle{ \hat{c}_n(0) e^{\mathrm{i} \pi n} = \frac{1}{N} \sum_{m=0}^{N-1} e^{-2\mathrm{i} \pi nm/N} c_0(x_m) }[/math]

The discrete Fourier transform

The discrete Fourier transform of a sequence of 2N complex numbers c₀, ..., c_2N−1 is transformed into the sequence of N complex numbers [math]\displaystyle{ \hat{c} }[/math]₀, ..., [math]\displaystyle{ \hat{c} }[/math]_N−1 by the DFT according to the formula:

[math]\displaystyle{ \hat{c}_m = \sum_{n=0}^{N-1} c_n e^{-2\pi \mathrm{i}mn/N} \quad \quad m = 0, \dots, N-1 }[/math]

We denote the transform by the symbol [math]\displaystyle{ \mathcal{F} }[/math], as in [math]\displaystyle{ \mathbf{X} = \mathcal{F} \left \{ \mathbf{x} \right \} }[/math] or [math]\displaystyle{ \mathcal{F} \left ( \mathbf{x} \right ) }[/math] or [math]\displaystyle{ \mathcal{F} \mathbf{x} }[/math].

The inverse discrete Fourier transform (IDFT) is given by

[math]\displaystyle{ c_n = \frac{1}{2N} \sum_{m=0}^{2N-1} \hat{c}_m e^{2\pi \mathrm{i}mn/2N} \quad \quad n = 0,\dots,N-1. }[/math]

Therefore we can write

[math]\displaystyle{ c(x_m,t) = \mathcal{F} \left\{ e^{-k_n^2 D t} \mathcal{F}^{-1} \left\{ c_0(x_m) \right\} \right\} }[/math]

The only difficulty is that we need to define carefully the values of [math]\displaystyle{ k_n }[/math]

The real power of this method lies with the Fast Fourier Transform or FFT algorithm. A naive implementation of the discrete Fourier transforms above will involve order [math]\displaystyle{ N^2 }[/math] operations. If the value of N is chosen to be a multiple of 2 this can be reduced to order [math]\displaystyle{ N \log(N) }[/math]. This is an incredible speed up, for example if N = 1024 the speed up is 147 (and this speed up just keeps increasing). This is reason this algorithm is used so extensively).

Reaction Diffusion Equations

We consider an auto catalytic reaction where the chemical species also diffuse. In this case the equations are

[math]\displaystyle{ \partial_t a = D\partial_x^2 a - kab }[/math]

[math]\displaystyle{ \partial_t b = D\partial_x^2 b + kab }[/math]

We can non-dimensionalise these equations scaling the variables as

[math]\displaystyle{ z = x/x^*\,\,\,\tau = t/t^*\,\,\,\alpha = a/a_0\,\,\,\beta = b/a^* }[/math]

So that the equations become

[math]\displaystyle{ \frac{1}{t^*}\partial_\tau \alpha = \frac{a_0}{(x^*)^2}D\partial_z^2 \alpha - k a_0^2 \alpha\beta }[/math]

[math]\displaystyle{ \frac{1}{t^*}\partial_\tau \beta = \frac{a_0}{(x^*)^2}D\partial_z^2 \beta + k a_0^2 \alpha\beta }[/math]

If we choose

[math]\displaystyle{ x^* = \sqrt{\frac{D}{ka_0}}\,\,\,t^*=\frac{1}{ka_0} }[/math]

then we obtain the system

[math]\displaystyle{ \partial_\tau \alpha = \partial_z^2 \alpha -\alpha\beta }[/math]

[math]\displaystyle{ \partial_\tau \beta =\partial_z^2 \beta + \alpha\beta }[/math]

Solution via split step method

Solutions for [math]\displaystyle{ \alpha(z,0) =1 \, \beta(z,0) = \exp(-10z^2) }[/math]

We can solve this equations numerically using a split step method. We assume that at time [math]\displaystyle{ \tau }[/math] we know [math]\displaystyle{ \alpha(z,\tau) }[/math] and [math]\displaystyle{ \beta(z,\tau) }[/math]. We then solve first the following equation

[math]\displaystyle{ \partial_\tau \alpha = \partial_z^2 \alpha }[/math]

from [math]\displaystyle{ \tau }[/math] to [math]\displaystyle{ \tau + \Delta\tau }[/math] (which we can do exactly using the spectral methods just discussed for the dispersion equation). We write this solution as [math]\displaystyle{ \tilde{\alpha}(z,\tau + \Delta\tau) }[/math] Then we solve

[math]\displaystyle{ \partial_\tau \alpha = -\alpha\beta }[/math]

by assuming that [math]\displaystyle{ \beta }[/math] is constant and subject to the boundary condition that [math]\displaystyle{ \alpha(z,\tau) = \tilde{\alpha}(z,\tau + \Delta\tau) }[/math]. This gives

[math]\displaystyle{ \alpha(z,\tau + \Delta\tau) = e^{-\beta(z,\tau) \Delta\tau} \tilde{\alpha}(z,\tau+ \Delta\tau)\, }[/math]

and we do likewise for the equation for [math]\displaystyle{ \beta }[/math]. Note that while both steps are exact the result from the split step method is an approximation with error which becomes smaller as the step size becomes smaller.

We can easily implement this split step method in matlab and we obtain a pair of travelling waves.

Travelling Waves solution

When we solve the equations we found the solution formed travelling waves and we now consider this phenomena in detail.

We define a new coordinate [math]\displaystyle{ y = z - v\tau }[/math] (so we will consider only waves travelling to the right, although we could analyse waves travelling to the left in a similar fashion). We seek stationary solutions in [math]\displaystyle{ \alpha(y) }[/math] and [math]\displaystyle{ \beta(y) }[/math] which satisfy

[math]\displaystyle{ \frac{\mathrm{d}^2 \alpha}{\mathrm{d}y^2} + v \frac{\mathrm{d} \alpha}{\mathrm{d}y} = \alpha\beta }[/math]

and

[math]\displaystyle{ \frac{\mathrm{d}^2 \beta}{\mathrm{d}y^2} + v \frac{\mathrm{d} \beta}{\mathrm{d}y} = -\alpha\beta }[/math]

If we add these equations we obtain

[math]\displaystyle{ \frac{\mathrm{d}^2 (\alpha+\beta)}{\mathrm{d}y^2} + v \frac{\mathrm{d} (\alpha+\beta)}{\mathrm{d}y} = 0 }[/math]

so that [math]\displaystyle{ \alpha + \beta = c_0 + c_1 e^{-vy} }[/math]. Boundary conditions are that as [math]\displaystyle{ y\to\infty }[/math] [math]\displaystyle{ \alpha = 1 }[/math] and [math]\displaystyle{ \beta = 0 }[/math]. We also require the solution to be bounded as [math]\displaystyle{ y\to-\infty }[/math] so that [math]\displaystyle{ \alpha + \beta = 1 }[/math]. This means that, since [math]\displaystyle{ \alpha \geq 0 }[/math], we must have [math]\displaystyle{ 0\leq \beta \leq 1 }[/math]. We can then obtain the following equation

[math]\displaystyle{ \frac{\mathrm{d}^2 \beta}{\mathrm{d}y^2} + v \frac{\mathrm{d} \beta}{\mathrm{d}y} + \beta(1-\beta)= 0 }[/math]

which we can write as the system of first order equations.

Phase portrait for out system showing the equilibrium points and the heteroclinic connection

We define the variable [math]\displaystyle{ \gamma = \frac{\mathrm{d}\beta}{\mathrm{d}y} }[/math] and we obtain

[math]\displaystyle{ \begin{align} \frac{\mathrm{d}\beta}{\mathrm{d} y} &= \gamma&\\ \frac{\mathrm{d}\gamma}{\mathrm{d} y} &= -v\gamma + \beta(\beta -1)& \\ \end{align} }[/math]

This dynamical system has equilibrium points at [math]\displaystyle{ (0,0) }[/math] and [math]\displaystyle{ (1,0) }[/math]. We can analyse these equilibrium points by linearization. The Jacobian matrix is

[math]\displaystyle{ J =\begin{pmatrix} 0 & 1 \\ -1 + 2\beta & -v \end{pmatrix} }[/math]

We can easily see that the Jacobian evaluated at our first equilibrium point is

[math]\displaystyle{ J_{(0,0)} =\begin{pmatrix} 0 & 1 \\ -1 & -v \end{pmatrix} }[/math]

which has eigenvalues [math]\displaystyle{ \mu_{\pm} = -1/2 (v \mp \sqrt{v^2-4}) }[/math]. Therefore this point is a nodal sink (possibly a spiral)

Additionally,

[math]\displaystyle{ J_{(1,0)} =\begin{pmatrix} 0 & 1 \\ 1 & -v \end{pmatrix} }[/math]

which has eigenvalues [math]\displaystyle{ \lambda_{\pm} = -1/2 (v \mp \sqrt{v^2+4}) }[/math]. This is a a saddle point. The unstable and stable separatrices leave the equilibrium point at [math]\displaystyle{ (1,0) }[/math] in the directions [math]\displaystyle{ \begin{pmatrix}\lambda_{\pm} \\ 1\end{pmatrix} }[/math]. The only path on which [math]\displaystyle{ \beta }[/math] is bounded as [math]\displaystyle{ y\to-\infty }[/math] are the unstable separatrices. Also, only the unstable separatrix which enters the region [math]\displaystyle{ \beta\lt 1 }[/math] is physically meaningful.

To find a travelling wave we need to find a heteroclinic connection between the two equilibrium points which also has to satisfy the conditions that [math]\displaystyle{ 0\leq \beta \leq 1 }[/math].

We need to show that the heteroclinic connection does not cross the [math]\displaystyle{ \beta }[/math] axis. Consider the region

[math]\displaystyle{ R = \left\{(\beta,\gamma)\,|\, \beta\lt 1,\,-k\beta\lt \gamma\lt 0\right\} }[/math]

On the line [math]\displaystyle{ \beta = 1,d\beta/dy\lt 0 }[/math] and hence all flow it into [math]\displaystyle{ R }[/math]. On the line [math]\displaystyle{ \gamma = 0, d\gamma/dy \lt 0 }[/math] for [math]\displaystyle{ 0\lt \beta\lt 1 }[/math]. On the line [math]\displaystyle{ \gamma = -k \beta }[/math] we know that [math]\displaystyle{ d\beta/dy \lt 0 }[/math] so that integral paths enter the region if and only if [math]\displaystyle{ d\gamma/d\beta \lt \gamma/\beta }[/math]. We know that

[math]\displaystyle{ d\gamma/d\beta - \gamma/\beta = -v - \frac{\beta(1-\beta)}{\gamma} -\frac{\gamma}{\beta} = \frac{1}{k} (k^2 - vk +1 -\beta) }[/math]

when [math]\displaystyle{ \gamma = -k\beta }[/math]