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{dp}{dt} = -kp\,\,\,\textrm{and}\,\,\, \frac{da}{dt} = 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{da}{dt} = -kab\,\,\,\textrm{and}\,\,\, \frac{db}{dt} = kab }[/math]

We can solve these equations by observing that

[math]\displaystyle{ \frac{d(a+b)}{dt} = 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{db}{dt} = 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 \hat{c}₀, ..., \hat{c}_2N−1 by the DFT according to the formula:

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

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{d^2 \alpha}{dy^2} + v \frac{d \alpha}{dy} = \alpha\beta }[/math]

and

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

If we add these equations we obtain

[math]\displaystyle{ \frac{d^2 (\alpha+\beta)}{dy^2} + v \frac{d (\alpha+\beta)}{dy} = 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] We can then obtain the following equation

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

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

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

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

To find a travelling wave we need to find a hetero clinic connection between