User Guide

Systems of equations

A system is passed to solve() as a list of equations. The variables need not have the same order; odecast determines the highest derivative of each variable separately and builds the first-order system accordingly.

from odecast import var, Eq, solve

x, z = var("x"), var("z")
eqs = [
    Eq(x.d(2) + z, 0),   # second order in x
    Eq(z.d() - x, 0),    # first order in z
]
sol = solve(eqs, ivp={x: 1.0, x.d(): 0.0, z: 0.0}, tspan=(0, 5))

Each equation must be solvable for the highest derivative of one variable.

Declaring the order

The order of a variable may be declared explicitly. This is optional and serves as a consistency check: if an equation uses a higher derivative, an OrderMismatchError is raised.

w = var("w", order=1)

Vector and matrix variables

A variable with shape represents a vector- or matrix-valued function. Its components are scalar variables and can be used anywhere a scalar variable can:

u = var("u", shape=2)
u[0], u[1]          # components
u.d(), u.d(2)       # vector derivatives

M = var("M", shape=(2, 3))
M[1, 2]             # component in row 1, column 2

An equation involving a vector variable as a whole is expanded into one equation per component. The following are equivalent:

Eq(u.d(2) + u, 0)

[Eq(u[0].d(2) + u[0], 0),
 Eq(u[1].d(2) + u[1], 0)]

Arithmetic operations are element-wise, and all array-valued terms in one equation must have the same shape. Scalars and scalar variables are broadcast. Functions such as sin must be applied to components, sp.sin(u[0]), not to the variable as a whole.

Matrix products are written with @. The matrix may be a NumPy array, a nested list or a SymPy matrix, and may contain parameters (see Parameters):

A = np.array([[0.0, 1.0], [-2.0, -0.3]])
eq = Eq(u.d(), A @ u)                        # u' = A u

K = np.array([[2, -1], [-1, 2]])
eq = Eq(u.d(2) + K @ u, 0)                   # coupled oscillators

M = var("M", shape=(2, 2))
eq = Eq(M.d(), A @ M)                        # matrix equation M' = A M

Initial conditions for array-valued variables are given as (nested) sequences or NumPy arrays of matching shape:

sol = solve(
    Eq(u.d(2) + u, 0),
    ivp={u: [1.0, 0.5], u.d(): [0.0, -0.2]},
    tspan=(0, 10),
)

sol[u]      # shape (2, n)
sol[u[0]]   # shape (n,)

For a matrix variable, sol[M] has shape (2, 3, n).

Scalar and vector variables can be mixed in one system:

x = var("x")
u = var("u", shape=2)

eqs = [
    Eq(x.d(2) + x - u[0], 0),
    Eq(u[0].d() + u[1], x),
    Eq(u[1].d() + u[0], 0),
]
sol = solve(eqs, ivp={x: 1.0, x.d(): 0.0, u: [0.5, 0.0]}, tspan=(0, 5))

Nonlinear and time-dependent terms

Equations may contain powers and functions of the unknowns. Both SymPy and NumPy functions can be used:

import numpy as np
import sympy as sp

# Van der Pol oscillator: y'' - mu (1 - y^2) y' + y = 0
y = var("y")
mu = 1.0
eq = Eq(y.d(2) - mu*(1 - y**2)*y.d() + y, 0)

# Pendulum
theta = var("theta")
eq = Eq(theta.d(2) + sp.sin(theta), 0)   # or np.sin(theta)

The independent variable is exported as odecast.t and can be used in the same way:

from odecast import t

eq = Eq(y.d(2) + y, sp.cos(2*t))   # driven oscillator
sol = solve(eq, ivp={y: 0.0, y.d(): 1.0}, tspan=(0, 10))

The SymPy symbol for \(t\) is t.symbol, which is useful when working with the expressions returned by the SymPy backend.

Each equation must be solvable for the highest derivative of one of the variables; implicit equations such as sp.sin(y.d()) = y are solved with sympy.solve(), and the first solution is used.

Parameters

SymPy symbols in the equations are parameters. Their values are passed to solve() with params, keyed by symbol or by name:

k, c = sp.symbols("k c")
eq = Eq(y.d(2) + c*y.d() + k*y, 0)

sol = solve(eq, ivp={y: 1, y.d(): 0}, tspan=(0, 10), params={k: 4.0, c: 0.1})

Initial and boundary values may contain parameters as well, e.g. ivp={y: y0, y.d(): 0} with params={y0: 2.0, ...}.

The equations are reduced and compiled once; solving them again with other parameter values reuses the compiled system. This makes parameter studies cheap:

for kv in np.linspace(1, 4, 100):
    sol = solve(eq, ivp={y: 1, y.d(): 0}, tspan=(0, 10), params={k: kv, c: 0.1})

The SymPy backend keeps parameters without a value symbolic. Declaring assumptions helps SymPy simplify the result:

k = sp.Symbol("k", positive=True)
sol = solve(Eq(y.d(2) + k*y, 0), backend="sympy", ivp={y: 1, y.d(): 0})
sol.as_expr(y)   # cos(sqrt(k)*t)

Sensitivities

The derivatives of the solution with respect to parameters show how sensitive it is to them. They are computed together with the solution, from the exact Jacobians:

k = sp.Symbol("k")
sol = solve(Eq(y.d(2) + k*y, 0), ivp={y: 1, y.d(): 0}, tspan=(0, 10),
            params={k: 4.0}, sensitivities=[k])

sol.sensitivity(y, k)            # dy/dk at sol.t
sol.sensitivity(y, k, t=5.0)     # dy/dk at t = 5

Parameters that appear in the initial values are taken into account, e.g. ivp={y: y0, ...} with sensitivities=[y0]. Sensitivities are available for initial value problems without event resets, and not for algebraic variables. fit() uses them for exact derivatives of the residuals.

Boundary value problems

A two-point boundary value problem specifies conditions at both ends of the interval tspan. Conditions are given as a list of BC objects:

from odecast import BC

# y'' + y = 0, y(0) = 0, y(1) = 1
sol = solve(
    Eq(y.d(2) + y, 0),
    bvp=[BC(y, t=0, value=0), BC(y, t=1, value=1)],
    tspan=(0, 1),
)
sol.eval(y, 0.5)

Conditions may involve derivatives below the highest order, e.g. BC(y.d(), t=0, value=0), and array-valued variables, e.g. BC(u, t=0, value=[0, 1]). The total number of conditions must equal the total order of the system.

The problem is solved with scipy.integrate.solve_bvp(). Further keyword arguments are passed on (tol, max_nodes, bc_tol), and two more are available:

mesh

Number of nodes of the initial mesh (default 50), or the mesh itself.

guess

Initial guess, keyed like an ivp dictionary. Values may be numbers, arrays over the mesh, or functions of \(t\). States without a guess start at zero.

Nonlinear boundary value problems may have several solutions, and the guess determines which one is found. The Bratu problem \(y'' + e^y = 0\), \(y(0) = y(1) = 0\) has two:

eq = Eq(y.d(2) + sp.exp(y), 0)
bcs = [BC(y, t=0, value=0), BC(y, t=1, value=0)]

lower = solve(eq, bvp=bcs, tspan=(0, 1))   # y(0.5) = 0.14
upper = solve(eq, bvp=bcs, tspan=(0, 1),
              guess={y: lambda t: 12*t*(1 - t)})   # y(0.5) = 4.09

Unknown parameters

Parameters can be determined together with the solution by passing them as unknowns with a starting value. Each unknown requires one boundary condition in addition to one per first-order state. This covers eigenvalue problems:

# y'' + lam y = 0, y(0) = y(pi) = 0; y'(0) = 1 fixes the amplitude
lam = sp.Symbol("lam")
sol = solve(
    Eq(y.d(2) + lam*y, 0),
    bvp=[BC(y, t=0, value=0), BC(y, t=np.pi, value=0), BC(y.d(), t=0, value=1)],
    tspan=(0, np.pi),
    unknowns={lam: 3.8},
    guess={y: lambda t: np.sin(2*t) / 2},
)
sol.unknowns   # {lam: 4.000...}

The starting values of the unknowns and the guess for the solution select which eigenvalue is found.

Differential-algebraic equations

A variable that occurs without derivatives is algebraic. It is determined by the equations that contain no highest derivative, the algebraic equations. No special syntax is needed:

# Fast equilibrium A <=> B (B = K A) followed by a slow reaction B -> C.
# T = A + B is the total amount of A and B.
T, C, A, B = var("T"), var("C"), var("A"), var("B")
k, K = 0.5, 3.0
eqs = [
    Eq(T.d(), -k*B),
    Eq(C.d(), k*B),
    Eq(B, K*A),        # algebraic
    Eq(A + B, T),      # algebraic
]
sol = solve(eqs, ivp={T: 1.0, C: 0.0}, tspan=(0, 10))
sol[B]                 # algebraic variables are part of the solution

Initial conditions are only needed for the differential variables. The number of algebraic equations must equal the number of algebraic variables.

odecast reduces the system to ordinary differential equations in three ways, depending on its structure.

Algebraic variables determined directly (index 1)

If the algebraic equations can be solved for the algebraic variables and are linear in them, as above, the algebraic variables are eliminated symbolically. The result is an ordinary system; everything else (events, parameters, the SymPy backend, equilibria, fitting) works as usual.

Constraints on the states (higher index)

An algebraic equation may not contain any algebraic variable at all. The classic example is a pendulum of length \(L\) in Cartesian coordinates, where the string tension \(\lambda\) (per unit mass and length) keeps the bob on the circle:

x, y, lam = var("x"), var("y"), var("lam")
eqs = [
    Eq(x.d(2), -lam*x),
    Eq(y.d(2), -lam*y - 9.81),
    Eq(x**2 + y**2, 1),    # constraint without lam
]
sol = solve(eqs, ivp={x: 0.6, x.d(): 0, y: -0.8, y.d(): 0}, tspan=(0, 10))

Such a constraint is differentiated with respect to \(t\) until an algebraic variable appears: here twice, after which \(\lambda\) follows from the equations of motion. The index of the system is the number of differentiations plus one; the pendulum has index 3.

The initial values must satisfy the constraint and the derivatives that were taken, i.e. the bob starts on the circle with a tangential velocity. odecast checks this and reports inconsistent initial values.

Numerically, only the differentiated constraint is enforced, so the solution would slowly drift away from the original one. odecast damps this drift with Baumgarte stabilization: the \(k\)-th derivative \(c^{(k)} = 0\) of a constraint \(c\) is replaced by \((d/dt + \alpha)^k c = 0\). The rate \(\alpha\) is chosen as five times the largest absolute eigenvalue of the Jacobian at the initial point, so that violations decay faster than the dynamics evolve. solve(..., baumgarte=alpha) sets it explicitly; baumgarte=0 disables the stabilization. The remaining violation is available as sol.constraint_residuals(); it is of the order of the solver tolerance.

Algebraic equations without closed-form solution

If the (possibly differentiated) algebraic equations are not linear in the algebraic variables, they are solved numerically with Newton’s method each time the right-hand side is evaluated, starting from the previous solution. Values given for algebraic variables in ivp are used as the starting guess at the initial time. Where an equation has several solutions, the guess selects one:

y, z = var("y"), var("z")
eqs = [Eq(y.d(), -z), Eq(z**2, y)]

solve(eqs, ivp={y: 1.0, z: 0.5}, tspan=(0, 1))    # follows z = +sqrt(y)
solve(eqs, ivp={y: 1.0, z: -0.5}, tspan=(0, 1))   # follows z = -sqrt(y)

For implicit solvers, the Jacobian is computed from the implicit function theorem. If the algebraic equations have no solution at some point, the integration stops with an error that names the time.

Limitations:

• Only constraints that can be differentiated one at a time are reduced. Systems whose algebraic equations remain singular after this (which would require the Pantelides algorithm) are rejected with an error.

• Boundary value problems, the SymPy backend and equilibria() support only the first case (index 1, linear).

• Events may involve algebraic variables only if they are eliminated symbolically.

Events

An event is a condition whose occurrences are located during integration, for example the time at which a falling body reaches the ground. Events are passed to solve() as a list of Event objects, or simply as equations:

from odecast import Event

y = var("y")
hits_ground = Event(Eq(y, 0), terminal=True, direction=-1)
sol = solve(Eq(y.d(2), -9.81), ivp={y: 10, y.d(): 0}, tspan=(0, 10),
            events=[hits_ground])

sol.t_events[0]               # [1.428]
sol.event_values(0, y.d())    # velocity at impact: [-14.01]

An event occurs where the difference of the two sides of its equation changes sign. terminal=True stops the integration at the first occurrence; direction restricts events to crossings from negative to positive (1) or from positive to negative (-1). Conditions may involve variables, their derivatives below the highest order, \(t\) and parameters. Plain functions g(t, x) of the state vector, as accepted by scipy.integrate.solve_ivp(), are passed through unchanged. An event is not triggered at the very start of the integration.

Resets

An event with reset changes the state when it occurs and then continues the integration; this describes impacts, switching and dosing. The new values are expressions in the values just before the event:

# A ball that loses 20% of its speed at every bounce
bounce = Event(Eq(y, 0), direction=-1, reset={y.d(): -0.8*y.d()})
sol = solve(Eq(y.d(2), -9.81), ivp={y: 1, y.d(): 0}, tspan=(0, 3),
            events=bounce)
sol.t_events[0]   # times of the bounces

# A dose at t = 1, 2, 3, ... (the zeros of sin(pi t))
c = var("c")
dose = Event(Eq(sp.sin(sp.pi*t), 0), reset={c: c + 1})
sol = solve(Eq(c.d(), -0.5*c), ivp={c: 1}, tspan=(0, 10), events=dose)

The solution, including sol.eval, is discontinuous at the events. Only states (variables and derivatives below the highest order) can be reset. If events accumulate, as for a bouncing ball coming to rest, the integration stops with an error after many resets; shorten tspan or add a terminal event.

Backends

The backend argument of solve() selects the solver.

"scipy" (default)

Numerical solution with scipy.integrate.solve_ivp() for initial value problems and scipy.integrate.solve_bvp() for boundary value problems. Requires ivp or bvp, and tspan. Further keyword arguments to solve are passed on to the SciPy solver, for example method, rtol, atol, max_step, first_step or t_eval:

solve(eq, ivp=ic, tspan=(0, 10), method="Radau", rtol=1e-8)

For the implicit methods "Radau", "BDF" and "LSODA", which are suited for stiff problems, odecast supplies the exact Jacobian, computed symbolically.

"sympy"

Symbolic solution with sympy.dsolve(). Initial or boundary conditions, if given, determine the integration constants; initial conditions apply at tspan[0] if tspan is given, otherwise at \(t = 0\). Conditions may also be incomplete, in which case the remaining constants stay in the result.

sol = solve(Eq(y.d(2) + y, 0), backend="sympy")
sol.as_expr(y)   # C1*sin(t) + C2*cos(t)

sol = solve(Eq(y.d(2) + y, 0), backend="sympy", ivp={y: 1, y.d(): 0})
sol.as_expr(y)   # cos(t)

Independent equations are solved one at a time. Coupled systems are passed to dsolve as a whole, which handles linear systems but generally not nonlinear ones.

"auto"

Solves numerically if ivp or bvp and tspan are given, and falls back to "sympy" if they are not or if the numerical solution fails. Invalid conditions are reported as errors rather than triggering the fallback.

Working with solutions

A numerical solution provides:

sol.t

The time points returned by the solver; set them with t_eval.

sol[y], sol[y.d()]

Values of a variable or derivative at sol.t.

sol.eval(y, times)

Values at arbitrary times within tspan. times may be a scalar or an array. The values come from the interpolant of the solver (dense_output of solve_ivp, or the cubic interpolant of solve_bvp) and are as accurate as the solution itself.

sol.plot(y, y.d())

Plots variables or derivatives against \(t\) with Matplotlib; all variables if none are given. Returns the axes.

sol.to_pandas()

A pandas.DataFrame with one column per variable, indexed by \(t\). Pass targets and t=... to choose columns and rows.

sol.t_events, sol.event_values(i, y)

Event times and solution values there (see Events).

sol.as_first_order()

The first-order system constructed by odecast: right-hand side function, Jacobian (or None), initial state vector, initial time, and a mapping from (variable, derivative order) to state vector index. Useful for debugging or for passing the system to a different integrator.

Matplotlib and pandas are optional; install them with pip install "odecast[examples]".

A symbolic solution provides sol.as_expr(target), which returns a SymPy expression for a variable or derivative (a SymPy array for vector variables), and sol.eval(target, times) once all integration constants are determined.

In Jupyter, equations and symbolic solutions are displayed as formulas. odecast.latex() returns the LaTeX code:

from odecast import latex

latex(Eq(y.d(2) + 0.3*y.d() + y, sp.sin(t)))
# y + 0.3 y' + y'' = \sin{\left(t \right)}

Errors

Problems with the equations or conditions raise exceptions derived from OdecastError. Conditions are checked before solving:

from odecast.errors import MissingInitialConditionError

y = var("y")
try:
    solve(Eq(y.d(2) + y, 0), ivp={y: 1.0}, tspan=(0, 10))
except MissingInitialConditionError as e:
    print(e)
# Missing initial condition for y^(1). Variable has order 2 and
# requires conditions for levels 0 to 1.

Supplying a condition for the highest derivative (e.g. y.d(2) above) raises OverdeterminedConditionsError. Invalid arguments raise ValueError or TypeError, for example an initial condition whose shape does not match the variable, or a boundary condition that is not at an end point of tspan. If the numerical solver fails, RuntimeError is raised.