Drivers (Time-dependent functions)

In some models, rates of change can be functions of time. This makes the ODE system “non-autonomous” - an “autonomous” system is one where time isn’t explicitly included in the equations. The trick that mathematicians use to handle this is to add an extra state variable that changes at a constant rate - effectively, a clock. In Cubie, we skip the middle-state and just allow you to use the symbol t in your equations, and it will be automatically evaluated with time as the integration proceeds. If you can neatly express your variable as a function of time, this is the cleanest way to get a forcing term into your problem. As an example, an experiment I was working on recently involved us shaking a tiny MEMS cantilever using a piezoelectric actuator under its base, measuring how much the cantilever bent as a result, then feeding that measured signal back into a heater on the cantilever that caused it to bend. That system could be modelled in Cubie like this:

import cubie as qb
import numpy as np

# Parameters
k = 0.1      # cantilever spring constant
c = 0.01     # cantilever damping coefficient (largely due to air)
alpha = 0.5  # heater coupling coefficient (how much bend the heat caused)
beta = 0.1   # heater dissipation coefficient (how quickly the heat dissipates)


def cantilever(t, y, p):
    f = 1e4 * t + 1e5
    base_wiggle = np.sin(2 * p.pi * f * t)
    dx = y.v
    dv = -p.k * y.x - p.c * y.v + p.alpha * y.T + base_wiggle
    dT = (p.alpha * (p.feedback_strength * y.x + p.feedback_offset)
          - p.beta * y.T)
    di = p.feedback_strength * y.x + p.feedback_offset
    return {"x": dx, "v": dv, "T": dT, "i": di}

parameters = {
    "k": 0.1,
    "c": 0.01,
    "alpha": 0.5,
    "beta": 0.1,
    "pi": np.pi,
    "feedback_strength": 0.5,
    "feedback_offset": 0.1,
}
initial_conditions = {
    "x": 0,    # initial position
    "v": 0,    # initial velocity
    "T": 0,    # initial temperature
    "i": 0,    # initial current
}

sys = qb.create_ODE_system(
    cantilever,
    parameters=parameters,
    states=initial_conditions,
    name="MEMSCantilever",
)

result = qb.solve_ivp(
    sys,
    y0=initial_conditions,
    parameters={
        "feedback_strength": np.linspace(-1, 1, 200),
        "feedback_offset": np.linspace(-1, 1, 200),
    },
    method="crank_nicolson",
    duration=1e-3,
)

This is an easy way to set up and parametrize a forcing term, and see how your system responds to different versions of it.

Arbitrary Values

Sometimes, your forcing function might not be easily expressed as a function of time. You might have some measured data that you want to test, you might want to use a random signal, or you might just want to throw some numbers in and see what comes out. If your integrator is fixed-step and you can supply the value of the forcing term at each time step, this is straightforward in theory - just pick a value for each time step, and use that. It’s a bit limiting to only use fixed-step algorithms, and if you’re running a long-duration simulation at a high time resolution, you might need to store a lot of time points. To give a bit more flexibility, Cubie can accept an array of forcing values and the time points to which they correspond, and will interpolate between them to get the value at any time point. This allows you to use adaptive-step algorithms, and also means you don’t need to store the value at every time point if your forcing function is smooth. Here’s an example of how to do this:

import cubie as qb
import numpy as np

# Create a measured signal as a driver
t_driver = np.linspace(0, 1.0, 1000)
signal = np.sin(2 * np.pi * 5 * t_driver) * np.exp(-t_driver)
drive = qb.DriverSamples({"drive_signal": signal}, time=t_driver)

def driven(t, y, p, d):
    dx = -p.k * y.x + p.amplitude * d.drive_signal
    return [dx]

sys = qb.create_ODE_system(
    driven,
    parameters={"k": 1.0, "amplitude": 1.0},
    states={"x": 0.0},
    drivers=["drive_signal"],
    name="DrivenSystem",
)

result = qb.solve_ivp(
    sys,
    y0={"x": np.array([0.0])},
    parameters={"amplitude": np.linspace(0.1, 2.0, 100)},
    drivers=drive,
    method="dormand-prince-54",
    duration=1.0,
)

Note how the driver appears: you declare it in drivers when you create the system, and read it by name inside the function body through the argument after the parameters, here d. Keep the parameter argument p even when the system has no parameters.

Sampled forcing data is defined using a DriverSamples object. It takes a dict of "[driver name]": [sampled data] pairs, with either the sample time vector given as time, or a starting time and sample spacing given as t0 and driver_sample_period. It is a solver setting, so Solver and Solver.update accept it too, and a solve without drivers reuses the last ones given.

Interpolation Options

Three further solver settings control how CuBIE interpolates between your samples:

order (default 3)

Polynomial degree of the spline fitted over each sample segment. The default cubic is smooth enough for most signals.

wrap (default True)

What happens outside the sampled time range. With wrap=True the signal repeats periodically — useful for periodic forcing. With wrap=False the signal is zero outside the sampled range; CuBIE adds a smooth spline segment from the last data point down to zero (and up from zero to the first point) so there is no discontinuous jump.

boundary_condition

End-point condition for the spline fit: one of "natural", "periodic", "clamped", or "not-a-knot" (the scipy CubicSpline conventions). If you don’t set it, CuBIE picks "periodic" when wrap=True and "clamped" when wrap=False.

result = qb.solve_ivp(
    sys,
    y0={"x": np.array([0.0])},
    parameters={"amplitude": np.linspace(0.1, 2.0, 100)},
    drivers=drive,
    order=3,
    wrap=True,
    method="dormand-prince-54",
    duration=1.0,
)

A worked example, including a sanity-check plot of the interpolated driver, lives in docs/source/examples/array_interpolation_example.py.