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 C1 and C2 interact according to the coupled reaction–diffusion equations: tC1=2C1+k(ka(kb+1)C1+C12C2) tC2=D2C2+k(kbC1C12C2)

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;


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 μ is the control parameter. For μ>0 the system is supercritical (Hopf bifurcation). We test several values of μ.

  μ = 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 C1=ka and C2=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});


Here we create an mpeg animation of the C1 concentration. The spread parameter sets the color scale to ± twice the standard deviation.

event movie (i = 1; i += 10)
  output_ppm (C1, linear = true, spread = 2, file = "f.mp4", n = 200);
  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 tC1=2C1+kka+k(C1C2kb1)C1 tC2=D2C2+kkbC1kC12C2 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, β);


We get the following stable Turing patterns.

Animation of the transitions
μ=0.04 μ=0.1 (stripes) μ=0.98 (hexagons)