Tutorial 1: Your First Parameter Sweep

This tutorial walks from an ODE definition to a heatmap showing how the system’s behaviour varies across a grid of parameter values. The whole workflow is about 30 lines of code.

The example system is the Lotka–Volterra predator–prey model, in which prey x grows and gets eaten while predators y eat and die off.

Step 1: Define the system

create_ODE_system() accepts either a plain Python function or a string of equations. The function form will feel familiar if you have used scipy.integrate.solve_ivp: it takes time, a state object, and a container of named values, and returns the derivatives in state order. States and named values can be accessed by attribute (y.x), by name (y["x"]), or by index (y[0]):

import numpy as np
import cubie as qb

def lotka_volterra(t, y, p):
    dx = p.a * y.x - p.b * y.x * y.y
    dy = -p.c * y.y + p.d * y.x * y.y
    return [dx, dy]

LV = qb.create_ODE_system(
    lotka_volterra,
    parameters={"a": 0.1, "b": 0.02, "c": 0.3, "d": 0.01},  # default values
    states={"x": 0.5, "y": 0.3},  # initial values
    name="LotkaVolterra",
)

CuBIE reads the function’s source code to build the GPU kernel, so the function form needs to live somewhere Python can read it back: a script, a module, or a notebook cell all work, but a string passed to exec or python -c does not.

The same system can be written as equation strings, which is handy when the model already exists as equations on paper. Anything of the form dx = ... defines a state variable x:

LV = qb.create_ODE_system(
    """
    dx = a*x - b*x*y
    dy = -c*y + d*x*y
    """,
    parameters={"a": 0.1, "b": 0.02, "c": 0.3, "d": 0.01},
    states={"x": 0.5, "y": 0.3},
    name="LotkaVolterra",
)

The keyword arguments sort your symbols into roles:

  • states are the variables the solver integrates. The dict values are default initial conditions; the y0 argument of solve_ivp() overrides them per run. For the string form the state names can also be inferred from the dx = ... left-hand sides; the function form requires states so CuBIE knows what the returned derivatives refer to.

  • parameters are the non-state constants in the model. Alongside initial values, they form the inputs you can batch over: each batch runs one integration per set of parameters and initial conditions. Every parameter that doesn’t vary over a given batch is compiled into the code as a number instead of loaded from memory, which can provide a large speedup in big systems.

  • Any right-hand-side symbol you never declared is inferred as a parameter with a default value of 0.0, and CuBIE emits a warning naming it. Declaring everything explicitly keeps the warning noise down and catches typos early.

Step 2: Sweep two parameters

Pass an array for each parameter you want to sweep. With the default grid_type="combinatorial", CuBIE solves every combination. Here that is 1000 x 1000 = 1,000,000 initial value problems, a batch size that a GPU handles comfortably:

b_values = np.linspace(0.01, 0.05, 1000)
d_values = np.linspace(0.005, 0.02, 1000)

result = qb.solve_ivp(
    LV,
    y0={"x": 0.5, "y": 0.3},
    parameters={"b": b_values, "d": d_values},
    method="ode45",
    duration=50.0,
    save_every=0.5,
)

Initial values accept plain floats; pass an array only when you want to sweep the initial state as well. method selects the integration algorithm, and "ode45" selects the adaptive 5th-order Dormand–Prince pair, the same method behind MATLAB’s ode45 and a good default for non-stiff problems. save_every=0.5 records a snapshot every half time-unit. If the batch outgrows GPU memory, CuBIE splits it into chunks along the run axis and solves the chunks in sequence, so batch size is limited by host memory rather than device memory.

Step 3: Look at the results

The returned SolveResult holds a 3-D array indexed [time, variable, run]:

trajectories = result.time_domain_array
print(trajectories.shape)        # (101, 2, 1000000)
print(result.time_domain_legend) # which variable is which index

# Trajectory of prey (variable 0) in the first run:
prey_run0 = trajectories[:, 0, 0]

Before trusting the numbers, check that every run integrated cleanly:

assert np.all(result.status_codes == 0), result.status_messages

Step 4: Map runs back to parameters

Combinatorial runs are laid out in C-order: the last parameter you passed (d) varies fastest. A reshape therefore recovers the 2-D parameter grid directly:

final_prey = trajectories[-1, 0, :]              # final x, per run
final_prey_grid = final_prey.reshape(1000, 1000)  # rows: b, cols: d

Step 5: Plot

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
im = ax.imshow(
    final_prey_grid,
    origin="lower",
    extent=[d_values[0], d_values[-1], b_values[0], b_values[-1]],
    aspect="auto",
)
ax.set_xlabel("d (predator growth per prey eaten)")
ax.set_ylabel("b (predation rate)")
ax.set_title("Final prey population")
fig.colorbar(im)
fig.savefig("lv_sweep.png", dpi=150)

That is the whole workflow: define the system once, sweep it in a single call, reshape, and plot.

Where to go next