BRUSSELATOR

Overview

The BRUSSELATOR function numerically solves the Brusselator system of ordinary differential equations, a theoretical model for autocatalytic chemical reactions that can exhibit sustained oscillations. The Brusselator (a portmanteau of “Brussels” and “oscillator”) was proposed by Nobel laureate Ilya Prigogine and collaborators at the Université Libre de Bruxelles in 1968.

The model describes the dynamics of two chemical species, X and Y, whose concentrations evolve according to the following system of differential equations:

\frac{dX}{dt} = a - (b + 1)X + X^2 Y

\frac{dY}{dt} = bX - X^2 Y

where a and b are positive rate constants. The system has a fixed point at (X, Y) = (a, b/a). The stability of this fixed point depends on the parameter values: when b > 1 + a^2, the fixed point becomes unstable and the system exhibits self-sustained oscillations approaching a limit cycle. Unlike the Lotka-Volterra equations, these oscillations are independent of initial conditions.

The Brusselator is closely related to real chemical oscillators such as the Belousov-Zhabotinsky reaction, a famous example of a clock reaction that exhibits visually striking color oscillations.

This implementation uses SciPy’s solve_ivp function to integrate the ODE system. The default Runge-Kutta method (RK45) works well for most cases, but alternative methods including RK23, DOP853, Radau, BDF, and LSODA are available for different accuracy requirements or stiff systems.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=BRUSSELATOR(x_initial, y_initial, a, b, t_start, t_end, timesteps, solve_ivp_method)

x_initial (float, required): Initial concentration of species X.
y_initial (float, required): Initial concentration of species Y.
a (float, required): Reaction rate parameter a (positive constant).
b (float, required): Reaction rate parameter b (positive constant).
t_start (float, required): Start time for integration.
t_end (float, required): End time for integration.
timesteps (int, optional, default: 10): Number of timesteps to solve for.
solve_ivp_method (str, optional, default: “RK45”): Integration method for solving the ODE system.

Returns (list[list]): 2D list with header [t, X, Y], or error message string.

Examples

Example 1: Demo case 1

Inputs:

x_initial	y_initial	a	b	t_start	t_end
1	1	1	2	0	1

Excel formula:

=BRUSSELATOR(1, 1, 1, 2, 0, 1)

Expected output:

t	X	Y
0	1	1
0.111	0.9002	1.106
0.222	0.8197	1.201
0.333	0.7544	1.29
0.444	0.7008	1.374
0.556	0.6566	1.454
0.667	0.6205	1.531
0.778	0.5921	1.603
0.889	0.5709	1.671
1	0.5544	1.736

Example 2: Demo case 2

Inputs:

x_initial	y_initial	a	b	t_start	t_end	timesteps
2	1	2	3	0	1	10

Excel formula:

=BRUSSELATOR(2, 1, 2, 3, 0, 1, 10)

Expected output:

t	X	Y
0	2	1
0.111	1.818	1.192
0.222	1.695	1.343
0.333	1.609	1.468
0.444	1.551	1.573
0.556	1.514	1.662
0.667	1.495	1.736
0.778	1.491	1.796
0.889	1.501	1.843
1	1.524	1.874

Example 3: Demo case 3

Inputs:

x_initial	y_initial	a	b	t_start	t_end	timesteps
1	2	1	2	0	2	10

Excel formula:

=BRUSSELATOR(1, 2, 1, 2, 0, 2, 10)

Expected output:

t	X	Y
0	1	2
0.222	1	2
0.444	1	2
0.667	1	2
0.889	1	2
1.111	1	2
1.333	1	2
1.556	1	2
1.778	1	2
2	1	2

Example 4: Demo case 4

Inputs:

x_initial	y_initial	a	b	t_start	t_end	timesteps	solve_ivp_method
1	1	1	2	0	1	10	RK23

Excel formula:

=BRUSSELATOR(1, 1, 1, 2, 0, 1, 10, "RK23")

Expected output:

t	X	Y
0	1	1
0.111	0.9001	1.106
0.222	0.8195	1.202
0.333	0.754	1.291
0.444	0.701	1.374
0.556	0.658	1.453
0.667	0.6234	1.527
0.778	0.5959	1.598
0.889	0.574	1.666
1	0.557	1.731

Python Code

from scipy.integrate import solve_ivp as scipy_solve_ivp
import numpy as np

def brusselator(x_initial, y_initial, a, b, t_start, t_end, timesteps=10, solve_ivp_method='RK45'):
    """
    Numerically solves the Brusselator system of ordinary differential equations for autocatalytic chemical reactions.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        x_initial (float): Initial concentration of species X.
        y_initial (float): Initial concentration of species Y.
        a (float): Reaction rate parameter a (positive constant).
        b (float): Reaction rate parameter b (positive constant).
        t_start (float): Start time for integration.
        t_end (float): End time for integration.
        timesteps (int, optional): Number of timesteps to solve for. Default is 10.
        solve_ivp_method (str, optional): Integration method for solving the ODE system. Valid options: RK45, RK23, DOP853, Radau, BDF, LSODA. Default is 'RK45'.

    Returns:
        list[list]: 2D list with header [t, X, Y], or error message string.
    """
    # Define valid methods
    valid_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']

    # Validate input types
    try:
        x0 = float(x_initial)
        y0 = float(y_initial)
        param_a = float(a)
        param_b = float(b)
        t0 = float(t_start)
        t1 = float(t_end)
        num_timesteps = int(timesteps)
    except Exception:
        return "Invalid input: All initial values, parameters, and timesteps must be numbers."
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    if num_timesteps <= 0:
        return "Invalid input: timesteps must be a positive integer."
    if solve_ivp_method not in valid_methods:
        return f"Invalid input: solve_ivp_method must be one of {valid_methods}."

    # Create time array for evaluation
    t_eval = np.linspace(t0, t1, num_timesteps)

    # Brusselator ODE system (defined inside main function per guidelines)
    def brusselator_ode(t, state):
        x, y = state
        dxdt = param_a - (param_b + 1) * x + x * x * y
        dydt = param_b * x - x * x * y
        return [dxdt, dydt]
    # Integrate
    try:
        sol = scipy_solve_ivp(
            brusselator_ode,
            [t0, t1],
            [x0, y0],
            method=solve_ivp_method,
            dense_output=False,
            t_eval=t_eval
        )
    except Exception as e:
        return f"Integration error: {e}"
    if not sol.success:
        return f"Integration failed: {sol.message}"

    # Helper function to check for valid numeric values
    def is_valid_number(val):
        return np.isfinite(val)

    # Format output: header row then data rows
    result = [["t", "X", "Y"]]
    for i in range(len(sol.t)):
        t_val = float(sol.t[i])
        x_val = float(sol.y[0][i])
        y_val = float(sol.y[1][i])
        # Disallow nan/inf
        if not all(is_valid_number(v) for v in [t_val, x_val, y_val]):
            return "Invalid output: nan or inf detected."
        result.append([t_val, x_val, y_val])
    return result

Overview

Excel Usage

Examples

Python Code

Online Calculator