BRUSSELATOR

Overview

The BRUSSELATOR function numerically solves the Brusselator system of ordinary differential equations, which models autocatalytic chemical reactions. The Brusselator is a classic example in nonlinear dynamics and chemical kinetics, describing the interaction of two species, X and Y, with parameters $a$ and $b$ controlling the reaction rates. The system is defined by:

\begin{align*} \frac{dx}{dt} &= a - (b + 1)x + x^2y \\ \frac{dy}{dt} &= b x - x^2y \end{align*}

where $x$ and $y$ are the concentrations of the two species, and $a$ , $b$ are positive constants. The function uses SciPy’s ODE solver (scipy.integrate.solve_ivp ) to integrate the system over a specified time interval. Only the basic two-variable model is exposed; spatial extensions and complex networks are not supported. The integration method can be selected from several common solvers, but only the most used options are exposed. This example function is provided as-is without any representation of accuracy.

Usage

To use the function in Excel:


=BRUSSELATOR(x_initial, y_initial, a, b, t_start, t_end, [method])

x_initial (float, required): Initial value of X.
y_initial (float, required): Initial value of Y.
a (float, required): Parameter a (reaction rate).
b (float, required): Parameter b (reaction rate).
t_start (float, required): Start time for integration.
t_end (float, required): End time for integration.
method (str, optional, default=‘RK45’): Integration method. One of RK45, RK23, DOP853, Radau, BDF, LSODA.

The function returns a 2D array (table) with columns: t, X, Y, representing time and the concentrations of X and Y at each step. If the input is invalid or integration fails, a string error message is returned.

Examples

Example 1: Basic Case

Inputs:

x_initial	y_initial	a	b	t_start	t_end	method
1.0	1.0	1.0	2.0	0.0	1.0	RK45

Excel formula:


=BRUSSELATOR(1.0, 1.0, 1.0, 2.0, 0.0, 1.0)

Expected output:

t	X	Y
0.000	1.000	1.000
0.091	0.916	1.088
0.804	0.587	1.620
1.000	0.554	1.736

Example 2: RK23 Method

Inputs:

x_initial	y_initial	a	b	t_start	t_end	method
1.0	1.0	1.0	2.0	0.0	1.0	RK23

Excel formula:


=BRUSSELATOR(1.0, 1.0, 1.0, 2.0, 0.0, 1.0, "RK23")

Expected output:

t	X	Y
0.000	1.000	1.000
0.019	0.982	1.018
0.204	0.832	1.187
0.447	0.700	1.376
0.736	0.605	1.572
1.000	0.557	1.731

Example 3: Different Parameters

Inputs:

x_initial	y_initial	a	b	t_start	t_end	method
2.0	1.0	2.0	3.0	0.0	1.0	RK45

Excel formula:


=BRUSSELATOR(2.0, 1.0, 2.0, 3.0, 0.0, 1.0)

Expected output:

t	X	Y
0.000	2.000	1.000
0.072	1.874	1.131
0.493	1.532	1.614
0.989	1.521	1.872
1.000	1.524	1.874

Example 4: Longer Time Interval

Inputs:

x_initial	y_initial	a	b	t_start	t_end	method
1.0	2.0	1.0	2.0	0.0	2.0	RK45

Excel formula:


=BRUSSELATOR(1.0, 2.0, 1.0, 2.0, 0.0, 2.0)

Expected output:

t	X	Y
0.000	1.000	2.000
0.000	1.000	2.000
0.000	1.000	2.000
0.000	1.000	2.000
0.001	1.000	2.000
0.011	1.000	2.000
0.111	1.000	2.000
1.111	1.000	2.000
2.000	1.000	2.000

Python Code


from typing import List, Union
from scipy.integrate import solve_ivp
 
def brusselator(x_initial: float, y_initial: float, a: float, b: float, t_start: float, t_end: float, method: str = 'RK45') -> Union[List[List[float]], str]:
    """
    Numerically solves the Brusselator system of ordinary differential equations for autocatalytic chemical reactions.
 
    Args:
        x_initial: Initial x value.
        y_initial: Initial y value.
        a: Parameter a.
        b: Parameter b.
        t_start: Start time of integration.
        t_end: End time of integration.
        method: Integration method ('RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'). Default is 'RK45'.
 
    Returns:
        2D list with header row: t, X, Y. Each row contains time and variable values, or an error message (str) if input is invalid.
 
    This example function is provided as-is without any representation of accuracy.
    """
    # Validate input types
    try:
        x0 = float(x_initial)
        y0 = float(y_initial)
        a = float(a)
        b = float(b)
        t0 = float(t_start)
        t1 = float(t_end)
    except Exception:
        return "Invalid input: All initial values and parameters must be numbers."
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    if method not in ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']:
        return "Invalid input: method must be one of 'RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'."
    # Brusselator ODE system
    def brusselator_ode(t, y):
        x, yv = y
        dxdt = a - (b + 1) * x + x * x * yv
        dydt = b * x - x * x * yv
        return [dxdt, dydt]
    # Integrate
    try:
        sol = solve_ivp(
            brusselator_ode,
            [t0, t1],
            [x0, y0],
            method=method,
            dense_output=False,
            t_eval=None
        )
    except Exception as e:
        return f"Integration error: {e}"
    if not sol.success:
        return f"Integration failed: {sol.message}"
    # 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 any([
            not isinstance(t_val, float) or not isinstance(x_val, float) or not isinstance(y_val, float),
            t_val != t_val or x_val != x_val or y_val != y_val,  # NaN check
            abs(t_val) == float('inf') or abs(x_val) == float('inf') or abs(y_val) == float('inf')
        ]):
            return "Invalid output: nan or inf detected."
        result.append([t_val, x_val, y_val])
    return result

Example Workbook

Link to Workbook