Integrator Architecture and Design

The numerical solution of differential equations is a cornerstone of QuTiP’s time-evolution solvers. Because different quantum systems present varying behaviours, no single numerical algorithm is optimal for all scenarios.

The Integrator framework provides a unified, backend-agnostic interface that wraps external or custom Ordinary Differential Equation (ODE) libraries. This design allows solvers to select routines from multiple distinct libraries through a single, interchangeable class.

Motivation

While Python packages like SciPy offer a robust collection of integration routines, they present several architectural limitations when applied directly to quantum dynamics:

Inconsistent API Interfaces: Across SciPy, interfaces vary drastically (e.g., the older scipy.integrate.ode vs. the modern scipy.integrate.solve_ivp), and some legacy algorithms do not natively support complex numbers without manual splitting into real and imaginary parts.
Hardware Limitations: Standard SciPy integration routines are strictly bound to CPU execution, preventing seamless acceleration via GPUs or distributed architectures.
Specialized Quantum Algorithms: Certain high-performance routines tailored for quantum mechanics —such as Krylov subspace methods or custom unitary propagators— are too domain-specific to be included in general-purpose mathematical libraries.

The Integrator class resolves these challenges by abstracting numerical backends into an interchangeable component, ensuring that switching from a SciPy Fortran solver to a GPU-accelerated backend requires only changing a configuration string.

Core Design

The abstract base class Integrator serves as the common template for all numerical integration backends.

Deterministic Time-Evolution

For standard solvers, the integrator acts sequentially. Given an initial time $t_0$ and an initial state vector or density matrix, it propagates the state forward to requested target times. This behaviour is governed by four core methods:

set_state(t, state): Configures the initial conditions of the integrator, seeding it with the starting time and the state array.
get_state: Returns the current internal numerical state of the integration loop as a tuple of (t, state).
integrate(t): Steps the numerical solver explicitly to the target time t using the internal state configuration.
run(tlist): A generator method that yields successive (t, state) tuples at each timestamp specified in the user-provided array tlist.

Note

The output state returned by get_state, integrate, and run is consistently packed as a tuple: (t, state).

Note

All inputs and outputs for states are handled as qutip.core.data.Data objects. However, an integrator is permitted to change the underlying representation format internally during evolution.

For example, SciPy-based integrators utilize NumPy arrays internally; they will automatically convert an incoming initial state from its initial storage type to a dense format during setup, and will subsequently return only dense states.

Monte Carlo Trajectories

To support Quantum Monte Carlo simulations (mcsolve), integrators must provide hook interfaces capable of handling non-continuous, conditional integration steps:

mcstep(t): Governs evolution inside stochastic loops. Unlike standard integration, mcstep must safely handle non-increasing time targets or backtrack queries (inside the last step range) to assist the parent solver in isolating the exact physical moment a quantum collapse event occurs.

For historical reasons, QuTiP delegates the root-finding search for a collapse to the high-level MCSolver loop rather than relying on native ODE “event” triggers. This design has the advantage of allowing QuTiP to wrap external ODE packages that lack native event-detection features, guaranteeing consistent behavior across all ODE providers.

Configuration and State Resetting

Numerical routines depend on parameters (such as tolerance limits atol/rtol, step limits nsteps, etc.) that must be tuned to achieve optimal efficiency. These parameters are configurable by the developer or end user via the options attribute passed down to the integrator instance. The Integrator.options attribute operates as a live dictionary containing the configuration state.

Because used libraries maintain internal state for their ODE solvers, modifying values within the options dictionary does not always take effect immediately. Activating the newly assigned options requires an explicit call to reset. This tears down the active backend instance and restarts it with the new options while preserving the active state and time.

Derivative Format

While ordinary ODE interfaces strictly expect functions to represent the right-hand side (RHS), quantum solvers in QuTiP can profit from alternative formats. Every integrator exposes an rhs_format class attribute string that dictates the type of RHS object required during initialization:

“callable”: The integrator expects a function matching the signature rhs(t: float, state: Data) -> Data returning the derivative of the state at time t. This is the default for most ODE methods. For specific Cython-backed ODE implementations, if the passed function is the method of QobjEvo.matmul_data, it will bypass standard Python overhead and hook directly into the Cython bindings of QobjEvo.matmul using a RHS object.
“matrix”: The integrator accepts a complex matrix formatted as a qutip.core.data.Data object, solving the linear matrix differential equation $frac{dX}{dt} = text{RHS} times X$.
“solver”: The integrator takes the parent Solver instance that instantiated it and builds the RHS internally. This is reserved for specialized integration methods that deliberately blend the underlying quantum physics with the numerical method.

Class Attributes Reference

Every subclass inheriting from Integrator must define the following class-level attributes:

integrator_options: A class-level dictionary defining the comprehensive set of configuration keys supported by the numerical backend along with their default values. Upon instantiation, this is copied to an instance-level options dictionary.

rhs_format: A string identifier defining the expected structural format of the Right-Hand Side (RHS) derivative function (e.g., "callable", "matrix", or "solver").

name: A descriptive, human-readable string representation of the algorithm (e.g., "SciPy zvode Adams solver") used primarily for populating metadata inside the final Result container.

Creating New Integrators

To implement a custom numerical backend, create a subclass derived from Integrator. At a minimum, your class must override integrate, set_state, get_state, and either __init__ or _prepare:

Let’s build a simple linear step integrator that updates the current state by adding dt * derivative(t, state). This integrator won’t provide very accurate results, but it will nicely illustrate what’s required.

from qutip.solver.integrator import Integrator

class LinearStepODE(Integrator):
    rhs_format = "callable"
    name = "Simple Linear Step"
    method = "linear"
    # An integrator may also define `integrator_options` here. If the attribute is not set,
    # a default empty dict is assumed, indicating that their are not special options for
    # this integrator.

    def _prepare(self):
        # __init__ automatically prepares the local options dict,
        # stores the rhs in self.derivative, and initializes self._is_set to False.
        # Apply any specialized options here!
        # This will be called from __init__ and reset.
        pass

    def set_state(self, t, state):
        self.t = t
        self.state = state
        self._is_set = True

    def integrate(self, t):
        # Basic Euler step example
        self.state = self.state + (t - self.t) * self.derivative(self.t, self.state)
        self.t = t

    def get_state(self):
        return self.t, self.state

The base Integrator class automatically provides default implementations for the run reset and mcstep methods (though the fallback version of mcstep can be very inefficient).

To expose your custom integrator to a specific physical solver, register it using the solver class’s target method:

from qutip.solver import SESolver

SESolver.add_integrator(LinearStepODE, "linear")

Once added, the method is immediately accessible via the high-level functional wrappers by configuring the options map (e.g., options={"method": "linear"}).

Registering the class with the parent Solver base class makes it globally available to all solvers except stochastic solver.

The RHS class

To optimize performance, certain Cython-backed integrators reuse memory when computing derivatives by supporting an in-place functional signature:

rhs(t: float, state: Data, out: Data) -> Data

In this paradigm, the operation should be interpreted as returning the computed derivative added directly into the existing out object. The function is permitted to reuse, modify, or overwrite the out container as needed. This pattern eliminates the considerable overhead of frequent memory allocations and deallocations inside tight numerical integration loops. The returned object does not need to the same as the out input if it’s immutable or otherwise not reusable.

These specialized integrators are marked as rhs_format="callable". Internally, they rely on the RHS helper class. This class acts as an interface layer that safely wraps an ordinary, out-of-place python derivative function and transforms it into one that behaves correctly when subjected to in-place calls.

Within the QuTiP’s solver, physical derivatives are usually computed via QobjEvo.matmul_data. The integrators automatically detect that method and bind directly to the low-level Cython implementation with in-place support.

If you are providing a custom Python callable, you must explicitly wrap your function inside an RHS instance to leverage in-place performance optimizations:

from qutip.solver.integrator._rhs import RHS

def custom_diff(t, state, out):
    # Perform in-place operations directly on the 'out' data object
    out += ...
    return out

derivative = RHS(custom_diff, inplace=True)