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_{-\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_{-\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{ 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{ 2N }[/math] points [math]\displaystyle{ x_m = -L + Lm/N }[/math]. We then use this to approximate the integral above and obtain

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

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

If we compare this to the formula for the discrete Fourier transform

The sequence of N complex numbers x₀, ..., x_N−1 is transformed into the sequence of N complex numbers X₀, ..., X_N−1 by the DFT according to the formula:

[math]\displaystyle{ X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}{N} k n} \quad \quad k = 0, \dots, N-1 }[/math]

where i is the imaginary unit and [math]\displaystyle{ e^{\frac{2 \pi i}{N}} }[/math] is a primitive N'th root of unity. (This expression can also be written in terms of a DFT matrix; when scaled appropriately it becomes a unitary matrix and the X_k can thus be viewed as coefficients of x in an orthonormal basis.)

The transform is sometimes denoted 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{ x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k e^{\frac{2\pi i}{N} k n} \quad \quad n = 0,\dots,N-1. }[/math]