# Coupled reaction–diffusion equations

The Brusselator is a theoretical model for a type of autocatalytic reaction. The Brusselator model was proposed by Ilya Prigogine and his collaborators at the Free University of Brussels.

Two chemical compounds with concentrations ${C}_{1}$ and ${C}_{2}$ interact according to the coupled reaction–diffusion equations: ${\partial }_{t}{C}_{1}={\nabla }^{2}{C}_{1}+k\left(ka-\left(kb+1\right){C}_{1}+{C}_{1}^{2}{C}_{2}\right)$ ${\partial }_{t}{C}_{2}=D{\nabla }^{2}{C}_{2}+k\left(kb{C}_{1}-{C}_{1}^{2}{C}_{2}\right)$

We will use a Cartesian (multi)grid, the generic time loop and the time-implicit diffusion solver.

``````#include "grid/multigrid.h"
#include "run.h"
#include "diffusion.h"``````

We need scalar fields for the concentrations.

``scalar C1[], C2[];``

We use the same parameters as Pena and Perez-Garcia, 2001

``````double k = 1., ka = 4.5, D = 8.;
double μ, kb;``````

The generic time loop needs a timestep. We will store the statistics on the diffusion solvers in `mgd1` and `mgd2`.

``````double dt;
mgstats mgd1, mgd2;``````

## Parameters

We change the size of the domain `L0` and set the tolerance of the implicit diffusion solver.

``````int main()
{
init_grid (128);
size (64);
TOLERANCE = 1e-4;``````

Here $\mu$ is the control parameter. For $\mu >0$ the system is supercritical (Hopf bifurcation). We test several values of $\mu$.

``````  μ = 0.04; run();
μ = 0.1;  run();
μ = 0.98; run();
}``````

## Initial conditions

``````event init (i = 0)
{``````

The marginal stability is obtained for `kb = kbcrit`.

``````  double ν = sqrt(1./D);
double kbcrit = sq(1. + ka*ν);
kb = kbcrit*(1. + μ);``````

The (unstable) stationary solution is ${C}_{1}=ka$ and ${C}_{2}=kb/ka$. It is perturbed by a random noise in [-0.01:0.01].

``````  foreach() {
C1[] = ka ;
C2[] = kb/ka + 0.01*noise();
}
boundary ({C1, C2});
}``````

## Outputs

Here we create an mpeg animation of the ${C}_{1}$ concentration. The `spread` parameter sets the color scale to $±$ twice the standard deviation.

``````event movie (i = 1; i += 10)
{
static FILE * fp1 = popen ("ppm2mpeg > f.mpg", "w");
output_ppm (C1, fp1, 200, linear = true, spread = 2);
fprintf (stderr, "%d %g %g %d %d\n", i, t, dt, mgd1.i, mgd2.i);
}``````

We make a PNG image of the final “pseudo-stationary” solution.

``````event final (t = 3000)
{
char name[80];
sprintf (name, "mu-%g.png", μ);
output_ppm (C1, file = name, n = 200, linear = true, spread = 2);
}``````

## Time integration

``````event integration (i++)
{``````

We first set the timestep according to the timing of upcoming events. We choose a maximum timestep of 1 which ensures the stability of the reactive terms for this example.

``  dt = dtnext (1.);``

We can rewrite the evolution equations as ${\partial }_{t}{C}_{1}={\nabla }^{2}{C}_{1}+k{k}_{a}+k\left({C}_{1}{C}_{2}-{k}_{b}-1\right){C}_{1}$ ${\partial }_{t}{C}_{2}=D{\nabla }^{2}{C}_{2}+k{k}_{b}{C}_{1}-k{C}_{1}^{2}{C}_{2}$ And use the diffusion solver to advance the system from $t$ to $t+dt$.

``````  scalar r[], β[];

foreach() {
r[] = k*ka;
β[] = k*(C1[]*C2[] - kb - 1.);
}
mgd1 = diffusion (C1, dt, r = r, β = β);
foreach() {
r[] = k*kb*C1[];
β[] = - k*sq(C1[]);
}
const face vector c[] = {D, D};
mgd2 = diffusion (C2, dt, c, r, β);
}``````

## Results

We get the following stable Turing patterns.

 $\mu =0.04$ $\mu =0.1$ (stripes) $\mu =0.98$ (hexagons)