API Reference

Public interface

odecast.var(name: str, order: int | None = None, shape: tuple | None = None) → Variable

Create a dependent variable for use in differential equations.

The variable represents an unknown function y(t) where t is the independent variable. You can create derivatives using the .d(n) method.

Parameters:

• name – Name of the variable (e.g., “y”, “x”, “theta”)

• order – Maximum order of derivatives expected (optional, for validation). If specified and equations use higher-order derivatives, an OrderMismatchError will be raised.

• shape – Shape tuple for vector/matrix variables (e.g., (3,) for 3D vector, (2, 2) for 2x2 matrix). None for scalar variables.

Returns:

• Arithmetic operations: y + z, 2*y, y/3, etc.

• Derivative notation: y.d() for y’, y.d(2) for y’’, etc.

• Use in equations: Eq(y.d(2) + y, 0)

• Component access for vectors: u[0], u[1] for vector components

• Vectorized operations: u.d() + 2*u for vector equations

Return type:

Variable object that supports

Examples

>>> y = var("y")           # Create scalar variable y(t)
>>> y.d()                  # First derivative y'(t)
>>> y.d(2)                 # Second derivative y''(t)
>>> Eq(y.d(2) + y, 0)     # Simple harmonic oscillator

>>> u = var("u", shape=(2,))  # Create 2D vector variable u(t)
>>> u[0]                      # First component u[0](t)
>>> u.d()                     # Vector derivative u'(t)
>>> Eq(u.d() + u, 0)         # Vector equation u' + u = 0

>>> z = var("z", order=1)  # Declare z as first-order only
>>> Eq(z.d() - z, 0)       # Valid: uses only first derivative
>>> Eq(z.d(2) + z, 0)      # Would raise OrderMismatchError

class odecast.Eq(lhs, rhs=0)

Bases: object

Represents an ordinary differential equation.

The Eq class is used to create equations of the form lhs = rhs, where the left and right sides can contain variables, derivatives, and expressions. This mirrors SymPy’s Eq class but works with odecast’s Variable and Derivative objects.

Examples

Simple harmonic oscillator: >>> y = var(“y”) >>> eq = Eq(y.d(2) + y, 0)

Damped oscillator with forcing: >>> eq = Eq(y.d(2) + 0.3*y.d() + y, sp.sin(t))

First-order ODE: >>> eq = Eq(y.d() - 2*y, 0)

Coupled system: >>> y, z = var(“y”), var(“z”) >>> eq1 = Eq(y.d() - z, 0) >>> eq2 = Eq(z.d() + y, 0)

sympy() → Equality

Convert this equation to a SymPy equation.

Returns:

SymPy Eq object

odecast.t

The distinguished independent variable, typically time.

It can be used directly in expressions and in SymPy functions, e.g. sp.sin(2*t).

odecast.solve(equation, *, ivp=None, bvp=None, tspan=None, backend=None, params=None, events=None, sensitivities=None, **kwargs)

Solve an ordinary differential equation or system of ODEs.

This is the main entry point for solving differential equations. It supports both symbolic and numeric solutions through different backends.

Parameters:

• equation – Single Eq object or list of Eq objects representing the ODE system

• ivp – Dictionary of initial conditions for initial value problems (IVP). Keys can be Variable objects (for y(t0)) or Derivative objects (for y’(t0)). Example: {y: 1.0, y.d(): 0.0} for y(0)=1, y’(0)=0

• bvp – List of BC objects for two-point boundary value problems (BVP). Example: [BC(y, t=0, value=0), BC(y, t=1, value=1)]

• tspan – Tuple of (t_start, t_end). Required for numeric solutions; initial conditions apply at t_start. For the sympy backend it is optional and only sets the time of the initial conditions (default 0).

• backend –

  Solution backend to use: - “scipy”: Numeric solution using SciPy (default) - “sympy”: Symbolic solution using SymPy - “auto”: Try SciPy first, fall back to SymPy if needed - “scipy_bvp”: Same as “scipy” with bvp

  SciPy uses solve_ivp for an IVP and solve_bvp for a BVP.

• **kwargs – Additional backend-specific options: For an IVP: method, rtol, atol, max_step, etc. For a BVP: guess, mesh, tol, max_nodes, bc_tol

Returns:

• SolutionIVP: For numeric IVPs, supports sol[y], sol[y.d()], etc.

• SolutionBVP: For numeric BVPs, with the same interface

• SolutionExpr: For symbolic solutions, supports sol.as_expr(y)

Return type:

Solution object that depends on the backend

Raises:

• MissingInitialConditionError – When required initial conditions are missing

• OverdeterminedConditionsError – When too many initial conditions provided

• OrderMismatchError – When declared variable order conflicts with equation usage

• BackendError – When the chosen backend cannot solve the equation

• ValueError – When required arguments are missing (e.g., tspan for scipy)

Examples

Numeric solution of second-order ODE: >>> y = var(“y”) >>> eq = Eq(y.d(2) + 0.3*y.d() + y, 0) >>> sol = solve(eq, ivp={y: 1.0, y.d(): 0.0}, tspan=(0, 10)) >>> y_values = sol[y] # Get y(t) values >>> yprime_values = sol[y.d()] # Get y’(t) values

Symbolic solution: >>> eq = Eq(y.d(2) + y, 0) >>> sol = solve(eq, backend=”sympy”, ivp={y: 1, y.d(): 0}) >>> expr = sol.as_expr(y) # cos(t)

Boundary value problem: >>> sol = solve(eq, bvp=[BC(y, t=0, value=0), BC(y, t=1, value=1)], tspan=(0, 1))

Coupled system: >>> y, z = var(“y”), var(“z”) >>> eq1 = Eq(y.d(2) + z, 0) >>> eq2 = Eq(z.d() - y, 0) >>> sol = solve([eq1, eq2], ivp={y: 1.0, y.d(): 0.0, z: 0.0}, tspan=(0, 5))

class odecast.BC(variable, *, t=None, value=None)

Bases: object

Boundary condition for boundary value problems (BVP).

Boundary conditions specify the value of a variable or its derivative at specific points in the domain. This is used for BVP problems where conditions are given at multiple points rather than all at the initial point.

Parameters:

• variable – The variable or derivative to constrain (Variable or Derivative object)

• t – The time/position where the condition applies

• value – The required value at that position

Examples

>>> y = var("y")
>>> bc1 = BC(y, t=0, value=0)      # y(0) = 0
>>> bc2 = BC(y, t=1, value=1)      # y(1) = 1
>>> bc3 = BC(y.d(), t=0, value=0)  # y'(0) = 0

BVP problem: >>> eq = Eq(y.d(2) + y, 0) >>> sol = solve(eq, bvp=[bc1, bc2], tspan=(0, 1))

Vector variables take a sequence of values: >>> u = var(“u”, shape=2) >>> BC(u, t=0, value=[0, 1])

class odecast.Event(condition, *, terminal: bool = False, direction: int = 0, reset: dict | None = None)

Bases: object

A condition whose occurrences are located during integration.

An event occurs when lhs - rhs of the condition changes sign. Its times are available as sol.t_events and the values of the solution there through sol.event_values.

Parameters:

• condition – An Eq, e.g. Eq(y, 0), or an expression that is zero at the event

• terminal – Stop the integration at the first occurrence

• direction – Only count zero crossings from negative to positive (1), from positive to negative (-1), or both (0)

• reset – New values for states at the event, e.g. {y.d(): -0.8 * y.d()}. Right-hand sides are evaluated with the values just before the event; integration then continues from the new state. Only variables and derivatives below the highest order can be reset.

Examples

>>> 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]

A bouncing ball that loses 20% of its speed at each bounce:

>>> bounce = Event(Eq(y, 0), direction=-1, reset={y.d(): -0.8 * y.d()})

expr() → Expr

The expression whose zeros are the event.

odecast.latex(expr) → str

LaTeX code for an equation or expression.

Parameters:

expr – Eq, odecast expression, or SymPy expression

Returns:

LaTeX string without surrounding $ signs

Analysis

odecast.equilibria(equation, params: dict | None = None, *, tol: float = 1e-10)

Find the equilibria of an autonomous system and classify their stability.

Equilibria are points where all derivatives vanish. They are found with sympy.solve(), so for periodic nonlinearities such as sin(y) only a finite selection is returned. Complex solutions are discarded.

Parameters:

• equation – Eq or list of Eq objects, as for odecast.solve()

• params – Values for parameters; parameters without a value stay symbolic in the results

• tol – Tolerance for treating an eigenvalue’s real part as zero

Returns:

List of Equilibrium objects

Raises:

ValueError – If the system depends explicitly on t, or SymPy cannot solve for the equilibria

Examples

>>> x, y = var("x"), var("y")
>>> equilibria([Eq(x.d(), x - x*y), Eq(y.d(), x*y - y)])
[Equilibrium(x=0, y=0, saddle), Equilibrium(x=1, y=1, center)]

class odecast.stability.Equilibrium(point, jacobian, eigenvalues, stability)

Bases: object

An equilibrium of an autonomous system.

point

Values of the variables; all their derivatives are zero

jacobian

Jacobian of the first-order system at the equilibrium

eigenvalues

Eigenvalues of the Jacobian (complex numbers if the Jacobian is numeric, SymPy expressions otherwise)

stability

Classification from the eigenvalues, e.g. “stable”, “unstable”, “stable focus”, “saddle”, “center”. “marginal” means the linearization is inconclusive; “unknown” is used if the Jacobian contains symbols.

odecast.linearize(equation, around, inputs=(), outputs=None, params: dict | None = None) → LinearSystem

Linearize a system around an operating point.

Parameters:

• equation – Eq or list of Eq objects. Inputs appear in them as SymPy symbols, e.g. a force F.

• around – Operating point: values of the variables and inputs, e.g. {x: 0, F: 0}, or an Equilibrium from odecast.equilibria(). Derivatives default to 0, inputs to 0.

• inputs – Input symbols

• outputs – Output expressions in the variables and their derivatives below the highest order; all variables if omitted

• params – Values of the remaining parameters; parameters without a value stay symbolic

Returns:

LinearSystem

Raises:

ValueError – If the operating point is not an equilibrium, or the system depends explicitly on t

Examples

>>> x = var("x")
>>> F = sp.Symbol("F")
>>> lin = linearize(Eq(2*x.d(2) + 0.5*x.d() + 8*x, F), around={x: 0},
...                 inputs=[F], outputs=[x])
>>> lin.transfer_function()   # 1/(2*s**2 + 0.5*s + 8)

class odecast.linearize.LinearSystem(A, B, C, D, states, inputs, outputs)

Bases: object

Linear state-space model x’ = A x + B u, y = C x + D u.

The states x, inputs u and outputs y are deviations from the operating point. Matrices are SymPy matrices; they are numeric unless parameters were left without a value.

A, B, C, D

State-space matrices

states

Names of the states, e.g. [“x”, “x’”]

inputs

Input symbols

outputs

Output expressions

eigenvalues() → list

Eigenvalues of A (the poles of the system).

to_scipy()

The model as scipy.signal.StateSpace (requires numeric matrices).

transfer_function(s: Symbol | None = None) → MutableDenseMatrix

Transfer function matrix G(s) = C (sI - A)^-1 B + D.

Parameters:

s – Laplace variable; sp.Symbol("s") if omitted

Returns:

SymPy matrix with one row per output and one column per input

odecast.phase_portrait(equation, x=None, y=None, *, xlim: tuple[float, float], ylim: tuple[float, float], params: dict | None = None, trajectories=(), tspan: tuple[float, float] = (0.0, 20.0), density: int = 25, nullclines: bool = True, equilibria: bool = True, ax=None)

Draw the phase portrait of a system with two states.

The system must be autonomous and have exactly two first-order states, e.g. two first-order equations, or one second-order equation (states y and y’).

Parameters:

• equation – Eq or list of Eq objects

• x – The states on the horizontal and vertical axis, e.g. theta and theta.d(); the two states in order if omitted

• y – The states on the horizontal and vertical axis, e.g. theta and theta.d(); the two states in order if omitted

• xlim – Plot ranges

• ylim – Plot ranges

• params – Parameter values

• trajectories – Initial values (dicts as for ivp) of trajectories to draw over tspan

• tspan – Time span of the trajectories

• density – Grid points per axis for the direction field and nullclines

• nullclines – Draw the curves on which x’ = 0 and y’ = 0

• equilibria – Mark the equilibria found by odecast.equilibria(), filled for stable, open for unstable, crosses for saddles. For periodic nonlinearities only representatives are found (e.g. 0 and pi for sin).

• ax – Matplotlib axes; the current axes if omitted

Returns:

The Matplotlib axes

Examples

>>> theta = var("theta")
>>> phase_portrait(Eq(theta.d(2) + 0.3*theta.d() + sp.sin(theta), 0),
...                xlim=(-7, 7), ylim=(-3, 3))

odecast.continuation(equation, param, start: dict, p_range: tuple[float, float], *, params: dict | None = None, step: float | None = None, max_points: int = 2000)

Follow a branch of equilibria as a parameter varies.

Uses pseudo-arclength continuation, so the branch is followed around folds where the parameter turns back. Folds, Hopf bifurcations and branch points are detected along the way and located precisely.

Parameters:

• equation – Eq or list of Eq objects of an autonomous system

• param – The parameter to vary (symbol or name)

• start – An equilibrium (or a point close to one) to start from: values of the variables and of param, e.g. {x: -1, r: -1}

• p_range – The branch is followed until the parameter leaves this interval; it starts in the direction of p_range[1]

• params – Values of the other parameters

• step – Initial arclength step; 1/200 of the parameter range by default

• max_points – Maximum number of points on the branch

Returns:

Branch

Examples

>>> x = var("x")
>>> r = sp.Symbol("r")
>>> branch = continuation(Eq(x.d(), r + x**2), r, {x: -1, r: -1}, (-1, 1))
>>> branch.bifurcations   # [Bifurcation(fold at param=0, x=0)]

class odecast.continuation.Branch(param, points, stable, eigenvalues, bifurcations, symbol)

Bases: object

A branch of equilibria, as returned by odecast.continuation().

param

Parameter values along the branch

points

Values of the variables along the branch, keyed by variable

stable

Whether each point is (linearly) stable

eigenvalues

Eigenvalues of the Jacobian at each point

bifurcations

Bifurcation points found on the branch

plot(variable, ax=None, **kwargs)

Bifurcation diagram: stable parts solid, unstable parts dashed.

Parameters:

• variable – Variable on the vertical axis

• ax – Matplotlib axes; the current axes if omitted

• **kwargs – Passed to ax.plot

Returns:

The Matplotlib axes

class odecast.continuation.Bifurcation(kind, param, point, eigenvalues)

Bases: object

A bifurcation point on a branch of equilibria.

kind

“fold” (saddle-node), “hopf”, or “branch point” (e.g. transcritical or pitchfork)

param

Parameter value

point

Values of the variables

eigenvalues

Eigenvalues of the Jacobian at the point

odecast.periodic_orbit(equation, ivp: dict, period: float, *, params: dict | None = None, transient: float = 0.0, mesh: int = 100, tol: float = 1e-06)

Find a periodic orbit (limit cycle) of an autonomous system.

A trajectory is computed from ivp for transient + period time units. Its last stretch of length period is the starting guess for a boundary value problem in which the period is unknown: x(0) = x(T), with the phase fixed by requiring that the derivative of the most strongly varying state vanishes at t = 0.

Parameters:

• equation – Eq or list of Eq objects

• ivp – Initial values for the starting trajectory

• period – Estimate of the period

• params – Parameter values

• transient – Time to integrate before the guess is taken, so that the trajectory approaches the orbit

• mesh – Number of mesh nodes for the boundary value problem

• tol – Tolerance for scipy.integrate.solve_bvp

Returns:

PeriodicOrbit

Raises:

• ValueError – If the system depends explicitly on t

• RuntimeError – If the boundary value problem does not converge

Examples

>>> y = var("y")
>>> orbit = periodic_orbit(Eq(y.d(2) - (1 - y**2)*y.d() + y, 0),
...                        ivp={y: 2, y.d(): 0}, period=6.5, transient=20)
>>> orbit.period      # 6.6633
>>> orbit.stability   # "stable"

class odecast.periodic.PeriodicOrbit(period, solution, multipliers, stability)

Bases: object

A periodic orbit, as returned by odecast.periodic_orbit().

period

The period T

solution

Solution over one period, t in [0, T]; supports sol[y], sol.eval(y, t) and sol.plot() like other solutions

multipliers

Floquet multipliers, the eigenvalues of the monodromy matrix. One of them is 1 (perturbations along the orbit).

stability

“stable” if all other multipliers lie inside the unit circle, “unstable” if one lies outside, “marginal” otherwise

odecast.fit(equation, data: dict, guess: dict, *, ivp: dict, tspan: tuple[float, float] | None = None, params: dict | None = None, bounds: dict | None = None, sigma: dict | None = None, **solve_options) → FitResult

Fit parameters of an initial value problem to measured data.

The parameters are SymPy symbols in the equations and, optionally, in the initial conditions. They are determined by nonlinear least squares (scipy.optimize.least_squares()).

Parameters:

• equation – Eq or list of Eq objects

• data – Measurements as {target: (times, values)}, where target is a variable or derivative. Several targets may be given, each with its own time points.

• guess – Parameters to fit with their starting values, e.g. {k: 1.0}

• ivp – Initial conditions, as for odecast.solve(). Values may be parameters, e.g. {y: y0} with y0 in guess.

• tspan – Integration interval; by default from the earliest to the latest measurement. Initial conditions apply at tspan[0].

• params – Values of parameters that are not fitted

• bounds – Optional {symbol: (lower, upper)}

• sigma – Optional measurement uncertainties {target: sigma} (scalar or array); residuals are divided by them

• **solve_options – Passed to odecast.solve(). The default tolerances are tightened to rtol=1e-8, atol=1e-10, since solver noise disturbs the optimization.

The Jacobian of the residuals is computed from the sensitivities of the solution (see sensitivities in odecast.solve()). Where these are not available, e.g. for algebraic variables computed numerically, finite differences are used.

Returns:

FitResult

Examples

>>> k = sp.Symbol("k")
>>> y = var("y")
>>> result = fit(Eq(y.d(), -k*y), {y: (t_data, y_data)},
...              guess={k: 1.0}, ivp={y: 1.0})
>>> result.params[k]

class odecast.fitting.FitResult(params, stderr, covariance, solution, result)

Bases: object

Result of odecast.fit().

params

Fitted parameter values, keyed by symbol

stderr

Standard errors of the fitted parameters

covariance

Covariance matrix, ordered like params

solution

Solution of the system with the fitted parameters

cost

Half the sum of squared (weighted) residuals

success

Whether the optimizer converged

message

Message from the optimizer

result

The full scipy.optimize.OptimizeResult

Delay and stochastic equations

odecast.solve_sde(equation, noise: dict, ivp: dict, tspan: tuple[float, float], *, dt: float, paths: int = 1, method: str = 'euler', params: dict | None = None, seed: int | None = None) → SDESolution

Simulate a stochastic differential equation (Ito interpretation).

The deterministic part is given by the equations, as for odecast.solve(). noise adds a term g dW to the equation of a state, with an independent Wiener process W for each entry:

• noise={y: sigma * y} for a first-order variable y gives dy = f dt + sigma y dW;

• for a variable of higher order the noise enters the equation of its highest derivative: noise={x.d(): sigma} for x’’ = … gives the Langevin equation dx’ = (…) dt + sigma dW.

Parameters:

• equation – Eq or list of Eq objects (the drift)

• noise – Noise intensities g, keyed by state (variables and derivatives below the highest order, like ivp)

• ivp – Initial values

• tspan – Time interval

• dt – Time step

• paths – Number of sample paths, simulated together

• method – “euler” (Euler-Maruyama, strong order 1/2) or “milstein” (strong order 1; uses dg/dx of each noise term with respect to its own state)

• params – Parameter values

• seed – Seed for the random number generator

Returns:

SDESolution

Examples

>>> y = var("y")
>>> sol = solve_sde(Eq(y.d(), 0.1*y), noise={y: 0.2*y}, ivp={y: 1.0},
...                 tspan=(0, 1), dt=1e-3, paths=1000, seed=1)
>>> sol.mean(y)[-1]   # close to exp(0.1)

class odecast.stochastic.SDESolution(t, states, mapping, W)

Bases: object

Ensemble of sample paths of a stochastic differential equation.

t

Time points (the integration grid)

W

Wiener processes that drove the paths, shape (number of noise terms, paths, len(t))

mean(target) → ndarray

Ensemble mean over the paths, as a function of t.

property paths: int

plot(target, paths: int = 20, ax=None, **kwargs)

Plot sample paths, the ensemble mean, and mean +- one standard deviation.

Parameters:

• target – Scalar variable or derivative

• paths – Number of sample paths to draw

• ax – Matplotlib axes; the current axes if omitted

• **kwargs – Passed to ax.plot for the sample paths

Returns:

The Matplotlib axes

quantile(target, q) → ndarray

Ensemble quantile(s) over the paths, as a function of t.

std(target) → ndarray

Ensemble standard deviation over the paths, as a function of t.

Partial differential equations

odecast.pde.solve_pde(equation, ic: dict, bc: dict, tspan: tuple[float, float], *, params: dict | None = None, acc: int = 2, **options) → PDESolution

Solve time-dependent partial differential equations with the method of lines.

Parameters:

• equation – Eq or list of Eq objects in fields, their time derivatives u.d()/u.d(2) and spatial derivatives u.dx(n), u.diff(x=1, y=1), u.laplace(). Coordinates (grid.x), t and parameters may appear as well.

• ic – Initial values per field: arrays of the grid shape, numbers, or functions of the coordinates. For equations of second order in time, also the initial time derivative: {u: ..., u.d(): ...}.

• bc – Boundary conditions per field and axis, e.g. {u: {"x": (Dirichlet(0), Neumann(0))}} (lower and upper end) or {u: {"x": Periodic()}}. Values may be functions of t.

• tspan – Time interval

• params – Parameter values

• acc – Accuracy order of the finite differences. Higher orders use one-sided stencils near non-periodic boundaries, which can make undamped wave equations unstable; acc=2 is safe there.

• **options – Passed to scipy.integrate.solve_ivp. The default method is “BDF”, with the sparsity pattern of the Jacobian.

Returns:

PDESolution

Examples

>>> grid = Grid(x=np.linspace(0, 1, 101))
>>> u = field("u", grid)
>>> sol = solve_pde(Eq(u.d(), 0.1 * u.dx(2)),
...                 ic={u: lambda x: np.sin(np.pi * x)},
...                 bc={u: {"x": (Dirichlet(0), Dirichlet(0))}},
...                 tspan=(0, 1))

class odecast.pde.Grid(**axes)

Bases: object

Rectangular grid of one to three dimensions.

Parameters:

**axes – Coordinates per axis, e.g. Grid(x=np.linspace(0, 1, 101)) or Grid(x=..., y=...). Non-uniform spacing is allowed except on periodic axes, whose grid must not contain the end point.

The coordinates are available as SymPy symbols for use in equations, e.g. grid.x.

axis(name) → int

odecast.pde.field(name: str, grid: Grid) → Field

Create an unknown field u(t, x, …) on a grid.

Examples

>>> grid = Grid(x=np.linspace(0, 1, 101))
>>> u = field("u", grid)
>>> Eq(u.d(), 0.1 * u.dx(2))   # heat equation

class odecast.pde.Dirichlet(value=0.0)

Bases: object

Prescribed value: a number, an array, or a function of t.

class odecast.pde.Neumann(value=0.0)

Bases: object

Prescribed derivative along the axis (not the outward normal).

class odecast.pde.Periodic

Bases: object

Periodic boundary; the grid must not contain the end point.

class odecast.pde.PDESolution(t, grid, fields, interpolant, unpack)

Bases: object

Solution of solve_pde().

t

Time points returned by the solver

grid

The grid

eval(f: Field, times) → ndarray

Values on the grid at arbitrary times.

Returns:

Array of shape grid.shape for a single time, or (*grid.shape, len(times))

plot(f: Field, times=None, ax=None, **kwargs)

Plot a field: profiles at several times on a 1D grid, or an image at one time on a 2D grid.

Parameters:

• f – The field

• times – Times to show; five equally spaced times (1D) or the final time (2D) if omitted

• ax – Matplotlib axes; the current axes if omitted

Returns:

The Matplotlib axes

Solutions

class odecast.solution.SolutionIVP(t: ndarray, y: ndarray, mapping: Any, t0: float, f_compiled: Callable | None = None, jac_compiled: Callable | None = None, x0: ndarray | None = None, interpolant: Callable | None = None, t_events: list[ndarray] | None = None, y_events: list[ndarray] | None = None)

Bases: object

Container for IVP solution data with convenient access methods.

as_first_order() → tuple[Callable, Callable | None, ndarray, float, dict]

Return first-order system representation.

Returns:

• f: Compiled RHS function f(t, x)

• jac: Compiled Jacobian function (or None)

• x0: Initial state vector

• t0: Initial time

• mapping: State mapping dictionary

Return type:

Tuple of (f, jac, x0, t0, mapping) where

constraint_residuals(t: ndarray | None = None) → ndarray

Residuals of the algebraic equations of a differential-algebraic system.

For constraints that odecast had to differentiate (higher index), the numerical solution drifts slowly away from the constraint; the residuals show how far.

Parameters:

t – Time points; the solver’s time points if omitted

Returns:

Array of shape (number of algebraic equations, number of time points)

eval(target, tpoints: float | ndarray) → float | ndarray

Evaluate solution at specific time points using interpolation.

Parameters:

• target – Variable, Derivative, or VectorDerivative to evaluate

• tpoints – Time point(s) to evaluate at

Returns:

Interpolated solution value(s). For vector and matrix targets the component axes come first.

event_values(index: int, target) → ndarray

Values of a variable or derivative at the occurrences of an event.

Parameters:

• index – Position of the event in the events list passed to solve

• target – Variable, Derivative, or vector/matrix variable or derivative

Returns:

Array with one value per occurrence (last axis)

plot(*targets, ax=None, n: int = 500, **kwargs)

Plot variables or derivatives against t using Matplotlib.

Parameters:

• *targets – Variables or derivatives to plot; all variables if omitted

• ax – Matplotlib axes to draw on; the current axes if omitted

• n – Number of evaluation points

• **kwargs – Passed to ax.plot

Returns:

The Matplotlib axes

sensitivity(target, param, t: float | ndarray | None = None)

Derivative of a variable or derivative with respect to a parameter.

Available for parameters passed as sensitivities to solve.

Parameters:

• target – Variable, Derivative, or vector/matrix variable or derivative

• param – Parameter symbol or name

• t – Time point(s); the solver’s time points if omitted

Returns:

Values of d(target)/d(param), shaped like eval or sol[target]

to_pandas(*targets, t: ndarray | None = None)

Tabulate the solution as a pandas DataFrame indexed by t.

Parameters:

• *targets – Variables or derivatives; all variables if omitted

• t – Time points; the solver’s time points if omitted

Returns:

pandas.DataFrame with one column per scalar component

variables() → list[Variable]

The scalar variables (and vector components) of the system.

class odecast.solution.SolutionBVP(t: ndarray, y: ndarray, mapping: Any, interpolant: Callable, f_compiled: Callable | None = None)

Bases: SolutionIVP

Solution of a boundary value problem.

Supports the same access methods as SolutionIVP. eval uses the cubic interpolant returned by scipy.integrate.solve_bvp.

class odecast.backends.sympy_backend.SolutionExpr(solutions: dict[Variable, Expr], t_symbol: Symbol)

Bases: object

Symbolic solution object that wraps SymPy expressions.

Without initial or boundary conditions the expressions contain integration constants C1, C2, ....

as_expr(target) → Expr

Get the SymPy expression for a variable or one of its derivatives.

Parameters:

target – Variable, Derivative, or vector/matrix variable or derivative

Returns:

SymPy expression in terms of t; a SymPy Array for vector and matrix targets

Raises:

KeyError – If variable not found in solution

eval(target, tpoints: float | ndarray) → float | ndarray

Evaluate the solution numerically.

Parameters:

• target – Variable or derivative, as for as_expr

• tpoints – Time point(s)

Returns:

Value(s) at tpoints; for vector targets the component axis comes first

Raises:

ValueError – If the solution still contains integration constants

Exceptions

Custom error types for odecast

exception odecast.errors.BackendError

Bases: OdecastError

Raised when a backend-specific error occurs.

exception odecast.errors.CompilationError

Bases: OdecastError

Raised when symbolic to numeric compilation fails.

exception odecast.errors.MissingInitialConditionError

Bases: OdecastError

Raised when required initial conditions are missing.

exception odecast.errors.NonSolvableFormError

Bases: OdecastError

Raised when the equation cannot be solved in the required form.

exception odecast.errors.OdecastError

Bases: Exception

Base exception for odecast errors.

exception odecast.errors.OrderMismatchError

Bases: OdecastError

Raised when variable order doesn’t match equation requirements.

exception odecast.errors.OverdeterminedConditionsError

Bases: OdecastError

Raised when too many initial/boundary conditions are provided.

Internals

The modules below implement the reduction to a first-order system. They are not part of the stable interface but are documented for contributors.

Equation analysis

Analysis functions for ODEs - order inference, validation, etc.

odecast.analyze.collect_variables(eqs: list[Eq]) → set[Variable]

Collect all variables used in the given equations.

Parameters:

eqs – List of equations to analyze

Returns:

Set of Variable objects found in the equations

odecast.analyze.infer_orders(eqs: list[Eq]) → dict[Variable, int]

Infer the maximum derivative order for each variable from the equations.

The equations are inspected in their SymPy form, so variables that only occur inside functions such as sp.sin(y) are found as well.

Parameters:

eqs – List of equations to analyze

Returns:

Dictionary mapping each variable to its maximum derivative order, in order of first appearance

odecast.analyze.resolve_orders(variables: list[Variable], inferred: dict[Variable, int]) → dict[Variable, int]

Resolve final orders by checking intent_order against inferred orders.

Parameters:

• variables – List of variables to resolve orders for

• inferred – Dictionary of inferred orders from equations

Returns:

Dictionary mapping each variable to its resolved order

Raises:

OrderMismatchError – If intent_order is set but is less than inferred order

Compilation

Compilation of symbolic expressions to numerical functions

class odecast.compile.CompiledSystem(f_sym_vec: MutableDenseMatrix, t_symbol: Symbol, extra_exprs=(), exclude=())

Bases: object

First-order system compiled to numerical functions.

Free symbols other than t and the states are parameters; their values are bound when the numerical functions are requested. The Jacobian is computed on first use.

For differential-algebraic systems (see odecast.dae), the algebraic variables have been eliminated from f_sym; algebraic_exprs gives them in terms of t and the states, invariants lists hidden constraints that the initial values must satisfy, and constraints the original algebraic equations.

STATE_PREFIX = '_odecast_x'

algebraic_functions(params: dict | None = None, z_guess: dict | None = None) → dict[str, Callable]

Functions z(t, X) for the algebraic variables, keyed by name.

check_consistency(t0: float, x0: ndarray, params: dict | None) → None

Check that initial values satisfy the hidden constraints.

Raises:

ValueError – If a constraint or one of its derivatives does not vanish at the initial point

constraint_function(params: dict | None = None, z_guess: dict | None = None) → Callable | None

Residuals g(t, X) of the original algebraic equations.

property is_dae: bool

jac(params: dict | None = None, z_guess: dict | None = None) → Callable

Jacobian J(t, x) for implicit solvers.

property jac_sym: MutableDenseMatrix

Symbolic Jacobian df/dx.

param_values(params: dict | None) → list[float]

Order parameter values like param_syms, checking completeness.

rhs(params: dict | None = None, z_guess: dict | None = None) → Callable

f(t, x) for scipy.integrate.solve_ivp.

rhs_vectorized(params: dict | None = None, z_guess: dict | None = None) → Callable

f(t, X) for arrays of time points, as scipy.integrate.solve_bvp requires.

rhs_vectorized_free(params: dict | None, names: list[str]) → Callable

f(t, X, p) for solve_bvp, where p holds the values of the parameters names, which are determined together with the solution.

sensitivity_rhs(params: dict | None, names: list[str]) → Callable

Right-hand side of the system augmented by its sensitivities.

The state is (x, S) with S = dx/dp for the parameters names, stored column by column after x. S obeys S’ = J S + df/dp.

odecast.compile.compile_system(eqs: list, orders: dict) → CompiledSystem

Reduce equations to a first-order system and compile it.

Results are cached by the structure of the equations, so solving the same system repeatedly (e.g. for different parameter values) compiles it only once.

Parameters:

• eqs – Scalar equations (vector equations already expanded)

• orders – Variable orders as returned by resolve_orders

Returns:

CompiledSystem whose states are ordered like build_state_map(orders)

odecast.compile.lambdify_jac(jac_sym: MutableDenseMatrix | None, t_symbol: Symbol, state_syms: list[Symbol], param_names: list[str] | None = None) → Callable | None

Compile symbolic Jacobian matrix to a numerical function.

Parameters:

• jac_sym – SymPy Matrix representing the Jacobian df/dx (can be None)

• t_symbol – SymPy symbol for time

• state_syms – List of SymPy symbols for state variables

• param_names – Optional list of parameter names

Returns:

Callable function jac(t, x) -> np.ndarray, or None if jac_sym is None

odecast.compile.lambdify_rhs(f_sym_vec: MutableDenseMatrix, t_symbol: Symbol, state_syms: list[Symbol], param_names: list[str] | None = None) → Callable

Compile symbolic RHS vector to a numerical function.

Parameters:

• f_sym_vec – SymPy Matrix representing the RHS vector f(t, x)

• t_symbol – SymPy symbol for time

• state_syms – List of SymPy symbols for state variables [x0, x1, …]

• param_names – Optional list of parameter names to include in the function

Returns:

Callable function f(t, x) -> np.ndarray

Order reduction

Reduction of higher-order ODEs to first-order systems

class odecast.reduce.StateMapping

Bases: dict

Enhanced state mapping that supports both Variable and (Variable, level) access patterns.

This class ensures that mapping[var][level] == mapping[(var, level)] for consistency across the API.

class odecast.reduce.VariableIndexer(mapping: StateMapping, variable: Variable)

Bases: object

Helper class that makes mapping[var][level] equivalent to mapping[(var, level)].

odecast.reduce.build_state_map(orders: dict[Variable, int]) → StateMapping

Build a state mapping for converting higher-order ODEs to first-order systems.

Parameters:

orders – Dictionary mapping each Variable to its maximum derivative order

Returns:

• mapping[var][level] returns state index for derivative level

• mapping[(var, level)] returns same state index

Return type:

StateMapping with consistent access patterns

odecast.reduce.isolate_highest_derivatives(eqs: list[Eq], orders: dict[Variable, int]) → dict[Variable, Expr]

Isolate the highest derivatives in the equations.

For each variable with order k>0, find an equation containing the k-th derivative and solve for that derivative in terms of t and lower derivatives.

Parameters:

• eqs – List of equations

• orders – Dictionary mapping variables to their orders

Returns:

Dictionary mapping each variable to a SymPy expression for its highest derivative

Raises:

NonSolvableFormError – If a highest derivative cannot be isolated

odecast.reduce.make_rhs(t_symbol: Symbol, mapping: dict[Variable | tuple, list[int] | int], highest_rules: dict[Variable, Expr], params: dict[str, float] | None = None, state_prefix: str = 'x') → tuple[MutableDenseMatrix, MutableDenseMatrix | None]

Build the right-hand side vector field f(t, x) for the first-order system.

Parameters:

• t_symbol – The SymPy symbol for the independent variable

• mapping – State mapping from build_state_map

• highest_rules – Expressions for highest derivatives from isolate_highest_derivatives

• params – Optional parameter values to substitute

• state_prefix – Name prefix of the state symbols (x0, x1, …)

Returns:

Tuple of (f_sympy_vector, jac_sympy_matrix) - f_sympy_vector: SymPy Matrix representing dx/dt = f(t, x) - jac_sympy_matrix: SymPy Matrix representing the Jacobian df/dx (optional for MVP)

Validation

Validation functions for ODEs and initial/boundary conditions

odecast.validate.normalize_conditions(ivp: dict | None = None, bvp: list | None = None, t0: float = 0.0) → list[tuple[Variable, int, float, object]]

Convert initial or boundary conditions to (variable, level, t, value) tuples.

Parameters:

• ivp – Initial conditions as passed to solve; they apply at t0

• bvp – List of BC objects

• t0 – Initial time for ivp

Returns:

List of (Variable, derivative level, time, value) for scalar components

odecast.validate.normalize_ivp(ivp_dict: dict[Variable | Derivative, float]) → dict[tuple[Variable, int], float]

Normalize IVP dictionary to use (Variable, level) keys.

Parameters:

ivp_dict – Dictionary with Variable (level 0), Derivative (level n≥1), or VectorDerivative keys. Values for vector and matrix variables are nested sequences or arrays of matching shape.

Returns:

Dictionary with (Variable, level) keys for all components

odecast.validate.validate_bvp(orders: dict[Variable, int], conditions: list[tuple[Variable, int, float, object]], tspan: tuple[float, float], n_unknowns: int = 0) → None

Validate boundary conditions for a two-point boundary value problem.

Each unknown parameter requires one condition in addition to one per first-order state.

Raises:

• ValueError – If a condition is not at an end point of tspan

• MissingInitialConditionError – If there are too few conditions

• OverdeterminedConditionsError – If there are too many or duplicate conditions

odecast.validate.validate_condition_levels(orders: dict[Variable, int], conditions: list[tuple[Variable, int, float, object]]) → None

Check that conditions refer to variables of the system and only to derivatives below each variable’s order.

Raises:

• ValueError – For a variable that does not occur in the equations

• OverdeterminedConditionsError – For a condition on the highest derivative

odecast.validate.validate_ivp(orders: dict[Variable, int], ivp: dict[Variable | Derivative, float], t0: float) → None

Validate that IVP has correct number of initial conditions.

Parameters:

• orders – Dictionary mapping variables to their orders

• ivp – IVP dictionary with initial conditions

• t0 – Initial time

Raises:

• MissingInitialConditionError – When required initial conditions are missing

• OverdeterminedConditionsError – When too many conditions are provided

odecast.validate.validate_variable_orders(equations: list[Eq], declared_orders: dict[Variable, int]) → None

Validate that declared variable orders match equation requirements.

Parameters:

• equations – List of equations to check

• declared_orders – Dictionary of explicitly declared variable orders

Raises:

OrderMismatchError – When declared order doesn’t match equation usage

Differential-algebraic equations

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). The system is reduced to ordinary differential equations for the remaining variables:

1. Algebraic equations that do not contain an algebraic variable are hidden constraints on the states. They are differentiated with respect to t, using the differential equations, until an algebraic variable appears (index reduction). The initial values must satisfy the constraint and its derivatives up to that point; drift away from the constraint is damped with Baumgarte stabilization, i.e. the k-th derivative c^(k) = 0 is replaced by (d/dt + alpha)^k c = 0.

2. If the resulting equations are linear in the algebraic variables, these are eliminated symbolically.

3. Otherwise they are computed numerically with Newton’s method whenever the right-hand side is evaluated.

class odecast.dae.NumericDAE(xdot, reduced, z_syms, t_symbol, extra_exprs=())

Bases: CompiledSystem

Differential-algebraic system whose algebraic variables are computed numerically.

Whenever the right-hand side is evaluated, the algebraic equations g(t, x, z) = 0 are solved for z with Newton’s method, starting from the previous solution. The Jacobian for implicit solvers follows from the implicit function theorem: dz/dx = -g_z^{-1} g_x.

algebraic_functions(params: dict | None = None, z_guess: dict | None = None) → dict[str, Callable]

Functions z(t, X) for the algebraic variables, keyed by name.

jac(params: dict | None = None, z_guess: dict | None = None) → Callable

Jacobian J(t, x) for implicit solvers.

property jac_sym

Symbolic Jacobian df/dx.

rhs(params: dict | None = None, z_guess: dict | None = None) → Callable

f(t, x) for scipy.integrate.solve_ivp.

rhs_vectorized(params: dict | None = None, z_guess: dict | None = None) → Callable

f(t, X) for arrays of time points, as scipy.integrate.solve_bvp requires.

rhs_vectorized_free(params, names)

f(t, X, p) for solve_bvp, where p holds the values of the parameters names, which are determined together with the solution.

sensitivity_rhs(params, names)

Right-hand side of the system augmented by its sensitivities.

The state is (x, S) with S = dx/dp for the parameters names, stored column by column after x. S obeys S’ = J S + df/dp.

odecast.dae.compile_dae(eqs: list, orders: dict, t_symbol: Symbol) → CompiledSystem

Reduce a differential-algebraic system and compile it.

Parameters:

• eqs – Scalar equations

• orders – Variable orders; algebraic variables have order 0

• t_symbol – SymPy symbol of the independent variable

Returns:

CompiledSystem, or NumericDAE if the algebraic variables must be computed numerically

Raises:

NonSolvableFormError – If the system cannot be brought into a form in which all highest derivatives and algebraic variables are determined

odecast.dae.default_baumgarte(system, params, t0, x0, z_guess) → float

Baumgarte rate for a system with hidden constraints.

Constraint violations should decay faster than the dynamics evolve, so the rate is five times the largest absolute eigenvalue of the Jacobian at the initial point (computed without stabilization).

Backends

SciPy-based backend for solving initial value problems

class odecast.backends.scipy_ivp.ScipyIVPBackend

Bases: object

Backend for solving IVPs using SciPy’s solve_ivp.

solve(f_compiled: Callable, jac_compiled: Callable | None, x0: ndarray, t0: float, tspan: tuple[float, float], mapping: dict[Variable | tuple[Variable, int], list | int], options: dict[str, Any] | None = None, after_reset: Callable | None = None) → SolutionIVP

Solve the IVP using SciPy’s solve_ivp.

Parameters:

• f_compiled – Compiled RHS function f(t, x) -> dx/dt

• jac_compiled – Compiled Jacobian function (optional)

• x0 – Initial state vector

• t0 – Initial time

• tspan – Tuple of (t_start, t_end)

• mapping – State mapping from build_state_map

• options – Keyword arguments for scipy.integrate.solve_ivp (method, rtol, atol, max_step, t_eval, first_step, events, …)

• after_reset – Called as after_reset(t, x) after an event has reset the state, e.g. to check the new state

Events with a reset attribute end the integration; the state is reset and the integration restarts until the end of tspan.

Returns:

SolutionIVP object with results

odecast.backends.scipy_ivp.convert_ivp_to_state_vector(ivp_dict: dict[Variable | Derivative, float], mapping: dict[Variable | tuple[Variable, int], list | int]) → ndarray

Convert IVP dictionary to initial state vector.

Parameters:

• ivp_dict – Dictionary mapping variables/derivatives to initial values

• mapping – State mapping from build_state_map

Returns:

Initial state vector x0

SciPy-based backend for solving boundary value problems

class odecast.backends.scipy_bvp.ScipyBVPBackend

Bases: object

Backend for solving two-point boundary value problems using SciPy’s solve_bvp.

solve(fun: Callable, n_states: int, conditions: list[tuple[Variable, int, float, object]], tspan: tuple[float, float], mapping: dict, guess: dict | None = None, mesh: int | ndarray = 50, p0: list[float] | None = None, **options) → SolutionBVP

Solve a boundary value problem.

Parameters:

• fun – Right-hand side f(t, X) of the first-order system, vectorized over mesh points

• n_states – Number of first-order states

• conditions – Validated (variable, level, t, value) boundary conditions

• tspan – Interval (a, b); all conditions lie at a or b

• mapping – State mapping from build_state_map

• guess – Initial guess, keyed like an IVP dictionary. Values may be numbers, arrays matching the mesh, or callables of t. Unspecified states start at zero.

• mesh – Number of initial mesh nodes, or the mesh itself

• p0 – Starting values of unknown parameters; fun then takes them as a third argument

• **options – Passed to scipy.integrate.solve_bvp (tol, max_nodes, bc_tol, verbose)

Returns:

SolutionBVP object with results

Raises:

RuntimeError – If solve_bvp does not converge

SymPy-based backend for symbolic ODE solving

class odecast.backends.sympy_backend.SymPyBackend

Bases: object

Backend that uses SymPy for symbolic ODE solving.

solve(equations: list[Eq], t_symbol: Symbol, conditions: list[tuple[Variable, int, float, object]] | None = None, params: dict | None = None, **options) → SolutionExpr

Solve ODEs symbolically using SymPy.

Independent equations, each involving a single variable, are solved one at a time. Coupled systems are passed to sympy.dsolve as a whole, which supports linear systems.

Parameters:

• equations – List of equations to solve

• t_symbol – SymPy symbol for independent variable

• conditions – Optional (variable, derivative level, t, value) tuples, used to determine the integration constants

• params – Optional values for parameters (SymPy symbols) in the equations; parameters without a value stay symbolic

• **options – Additional options (currently unused)

Returns:

SolutionExpr object with symbolic solutions

Raises:

BackendError – If SymPy cannot solve the system