Solving ODEs

Once you have defined an ODE system (see Creating a “System” of Differential Equations), CuBIE provides two ways to solve it: the convenience function solve_ivp() and the reusable Solver class.

Quick Start with solve_ivp

The fastest way to get results is solve_ivp(). Here we solve the Lotka–Volterra system for a single initial condition:

import cubie as qb
import numpy as np

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

result = qb.solve_ivp(
    LV,
    y0={"x": np.array([0.5]), "y": np.array([0.3])},
    parameters={"b": np.array([0.02]), "d": np.array([0.01])},
    method="dormand_prince_54",
    duration=100.0,
)

result is a SolveResult that holds the time-domain output, summaries, and metadata. See Working with Results for how to inspect it.

Note

CuBIE is a batch solver—it is designed to solve many IVPs at once. A single-system call works, but the real performance gains come from passing arrays of initial values or parameters. See Batching and Parameter Sweeps.

Changing Initial States

Pass arrays for any states you want to vary across the batch:

result = qb.solve_ivp(
    LV,
    y0={
        "x": np.linspace(0.1, 2.0, 500),
        "y": np.array([0.3]),  # same for all runs
    },
    parameters={"b": np.array([0.02]), "d": np.array([0.01])},
    method="dormand_prince_54",
    duration=100.0,
)

States provided as length-1 arrays are broadcast to all runs.

The Solver Class

If you plan to solve the same system multiple times (e.g. exploring different parameter sets interactively), create a Solver once and call solve() repeatedly. This avoids recompiling the CUDA kernel on every call:

solver = qb.Solver(LV, algorithm="dormand_prince_54")

result_a = solver.solve(
    initial_values={"x": np.array([0.5]), "y": np.array([0.3])},
    parameters={"b": np.linspace(0.01, 0.05, 1000),
                 "d": np.array([0.01])},
    duration=100.0,
)

result_b = solver.solve(
    initial_values={"x": np.array([1.0]), "y": np.array([0.5])},
    parameters={"b": np.linspace(0.01, 0.05, 1000),
                 "d": np.array([0.01])},
    duration=200.0,
)

The first call compiles the kernel; subsequent calls reuse it.

The duration Parameter

duration sets the total integration time. You can also specify t0 (default 0) and settling_time. When settling_time is non-zero, the solver integrates for that period before it starts recording output, useful for discarding transients.

Next Steps