# Brusselator

System type: polynomial continuous system

State dimension: 2

Application domain: Chemical kinetics

## Model description

A chemical reaction is said to be *autocatalytic* if one of the reaction products is also a catalyst for the same or a coupled reaction, and such a reaction is called an autocatalytic reaction. We refer to the wikipedia article Autocatalysis for details.

The Brusselator is a mathematical model for a class of autocatalytic reactions. The dynamics of the Brusselator is given by the two-dimensional ODE

\[ \left\{ \begin{array}{lcl} \dot{x} & = & A + x^2\cdot y - B\cdot x - x \\ \dot{y} & = & B\cdot x - x^2\cdot y \end{array} \right.\]

The numerical values for the model's constants (in their respective units) are given in the following table.

Quantity | Value |
---|---|

A | 1 |

B | 1.5 |

```
using ReachabilityAnalysis
const A = 1.0
const B = 1.5
const B1 = B + 1
@taylorize function brusselator!(du, u, p, t)
x, y = u
x² = x * x
aux = x² * y
du[1] = A + aux - B1*x
du[2] = B*x - aux
return du
end
```

The auxiliary variables `B1`

, `x²`

and `aux`

have been defined to make better use of `@taylorize`

and help to reduce allocations.

## Reachability settings

The initial set is defined by $x \in [0.8, 1]$, $y \in [0, 0.2]$. These settings are taken from ^{[CAS13]}.

```
U₀ = (0.8 .. 1.0) × (0.0 .. 0.2);
prob = @ivp(u' = brusselator!(U), u(0) ∈ U₀, dim: 2);
```

## Results

We use `TMJets`

algorithm with sixth-order expansion in time and second order expansion in the spatial variables.

`sol = solve(prob, T=18.0, alg=TMJets(orderT=6, orderQ=2));`

```
using Plots
plot(sol, vars=(1, 2), xlab="x", ylab="y", lw=0.2, color=:blue, lab="Flowpipe", legend=:bottomright)
plot!(U₀, color=:orange, lab="Uo")
```

We observe that the system converges to the equilibrium point `(1.0, 1.5)`

.

Below we plot the flowpipes projected into the time domain.

```
plot(sol, vars=(0, 1), xlab="t", lw=0.2, color=:blue, lab="x(t)", legend=:bottomright)
plot!(sol, vars=(0, 2), xlab="t", lw=0.2, color=:red, lab="y(t)")
```

## Changing the initial volume

The model was considered in ^{[GCLASG20]} but using a different set of initial conditions. Let us parametrize the initial states as a ball centered at $x = y = 1$ and radius $r > 0$:

`U0(r) = Singleton([1.0, 1.0]) ⊕ BallInf(zeros(2), r)`

U0 (generic function with 1 method)

The parametric initial-value problem is defined accordingly.

`bruss(r) = @ivp(u' = brusselator!(u), u(0) ∈ U0(r), dim: 2)`

bruss (generic function with 1 method)

First we solve for $r = 0.01$:

```
sol_01 = solve(bruss(0.01), T=30.0, alg=TMJets(orderT=6, orderQ=2))
LazySets.set_ztol(Float64, 1e-15)
plot(sol_01, vars=(1, 2), xlab="x", ylab="y", lw=0.2, color=:blue, lab="Flowpipe (r = 0.01)", legend=:bottomright)
plot!(U0(0.01), color=:orange, lab="Uo", xlims=(0.6, 1.3))
```