Brusselator Model Analysis: From ODEs to Bayesian Inference

This project provides a comprehensive numerical analysis of the Brusselator model, a theoretical framework for autocatalytic chemical reactions. The implementation spans deterministic solvers, stochastic differential equations (SDEs), and Bayesian parameter estimation using Markov Chain Monte Carlo (MCMC).

---

The Report

Report.pdf is the final submitted coursework. It achieved a grade of 82%.

---

Requirements

• Python 3.13.11
• See requirements.txt for full list of dependencies

Install dependencies

pip install -r requirements.txt

---

Usage

Run the full pipeline:

python src/pipeline.py

---

1. Deterministic Analysis (Forward Euler)

The first stage implements the Forward Euler method to solve the system of ODEs. By varying the time step $dt$, we analyze the order of convergence and verify that the numerical approximation approaches the true solution as the grid refines.

Figure 1: Comparison of X(t) approximations across different time steps.

Key Technical Details:

• Convergence: The global truncation error is measured at $T=50$. A log-log plot confirms the expected first-order accuracy.
• Uncertainty Quantification: We perform a Monte Carlo simulation (5,000 samples) to estimate the expected time-averaged variance of $X(t)$ when parameter $B$ is uncertain ($B \sim \text{Uniform}[1.5, 2.5]$).

---

2. Stochastic Modeling (Euler-Maruyama)

To account for molecular fluctuations in finite volumes, we extend the deterministic model into a Stochastic Differential Equation (SDE) using the Euler-Maruyama method.

Figure 2: Euler-Maruyama sample path illustrating stochastic oscillations around the deterministic limit.

Quantities of Interest (QoI):

We track two primary metrics over the interval $[30, 50]$:

1. $Q_1$: The maximum value of the reaction rate term $x^2y$.
2. $Q_2$: The time-averaged variance of $X(t)$.

The relationship between these stochastic variables is visualized using a 2D density histogram to identify correlations in the system's noise profile.

---

3. Bayesian Inference & MCMC

The final component addresses the inverse problem: estimating the parameter $B$ given noisy observation data at $T = {10, 15, 20}$.

The Bayesian Framework:

• Prior: $B \sim \mathcal{N}(2.0, 0.2^2)$.
• Likelihood: Based on Gaussian noise ($\sigma = 0.1$) between the Forward Euler observation operator and the provided data.
• Algorithm: We implement a Random Walk Metropolis-Hastings (RWMH) sampler with 100,000 iterations.

Figure 3: Trace plot of the RWMH chain showing mixing and convergence to the posterior distribution.

Post-Processing:

After removing a burn-in period of 2,000 samples, we calculate the posterior statistics and the 95% confidence interval for $B$.

Metric	Value
Posterior Mean ($B$)	~1.9123
Posterior Variance	~0.0016
95% Confidence Interval	Computed via sample std deviation

---