odecast

odecast is a Python library for solving ordinary differential equations of arbitrary order. Equations are written in their natural form, e.g. y'' + 0.3 y' + y = 0, and odecast performs the reduction to a first-order system that numerical integrators require. Numerical solutions are computed with SciPy, closed-form solutions with SymPy.

Full documentation: https://maroba.github.io/odecast

Installation

pip install odecast

This installs NumPy, SciPy and SymPy. Python 3.10 or later is required. Plotting and to_pandas() additionally need Matplotlib and pandas (pip install "odecast[examples]"), partial differential equations need findiff (pip install "odecast[pde]").

Example

The damped harmonic oscillator y'' + 0.3 y' + y = 0 with y(0) = 1, y'(0) = 0:

from odecast import var, Eq, solve

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

sol.t          # time points
sol[y]         # y(t)
sol[y.d()]     # y'(t)
sol.eval(y, 5.0)

For comparison, the same problem written directly against scipy.integrate.solve_ivp requires the state vector to be set up by hand:

def rhs(t, z):
    y, dy = z
    return [dy, -0.3 * dy - y]

solve_ivp(rhs, (0, 10), [1.0, 0.0])

This is manageable for a single equation but becomes tedious and error-prone for coupled systems of mixed order, which is the case odecast is designed for.

Features

• Arbitrary order. The order of each variable is inferred from the equations; no manual reduction is needed.
• Coupled systems of mixed order. Pass a list of equations; each variable may appear with a different highest derivative.
• Nonlinear equations. Powers and functions of the unknowns, such as y**2 or sp.sin(y) (or np.sin(y)), can be used directly.
• Vector and matrix variables. var("u", shape=2) or var("M", shape=(2, 2)) declares an array-valued unknown. Equations such as Eq(u.d(2) + u, 0) are expanded component-wise; matrix products are written A @ u.
• Parameters. SymPy symbols in the equations take values from params=. The system is compiled once and reused for further values.
• Initial and boundary value problems. Events are located during integration and can reset the state (impacts, switching, dosing). Two-point boundary value problems may have unknown parameters, e.g. eigenvalues.
• Differential-algebraic equations. Variables without derivatives are algebraic. Constraints such as Eq(x**2 + y**2, L**2) are handled by index reduction with stabilization; algebraic equations without closed-form solution are solved numerically.
• Stiff problems. Implicit solvers receive the exact Jacobian.
• Symbolic solutions. The SymPy backend returns closed-form solutions, with initial or boundary conditions applied.
• Sensitivities. Derivatives of the solution with respect to parameters, computed exactly alongside the solution.
• Analysis. Equilibria and their stability, linearization and transfer functions, phase portraits, continuation of equilibria with detection of bifurcations, periodic orbits with Floquet multipliers, and fitting parameters to measured data.
• Delay and stochastic equations. y.delay(tau) for delayed values; solve_sde() for equations with noise (Euler-Maruyama, Milstein).
• Partial differential equations. Time-dependent PDEs on grids in one to three dimensions, with the method of lines and finite differences from findiff.
• Validation. Missing, superfluous or inconsistent conditions are reported before solving, with a message naming the offending condition.

Coupled systems

x, y = var("x"), var("y")
eqs = [
    Eq(x.d() - x * (1 - y), 0),
    Eq(y.d() - y * (x - 1), 0),
]
sol = solve(eqs, ivp={x: 1.0, y: 0.5}, tspan=(0, 15))

Nonlinear and time-dependent terms

The independent variable is available as odecast.t:

import sympy as sp
from odecast import t

theta = var("theta")
eq = Eq(theta.d(2) + 0.1 * theta.d() + sp.sin(theta), 0.5 * sp.cos(2 * t))

Vector and matrix variables

u = var("u", shape=2)
eq = Eq(u.d(2) + 0.1 * u.d() + u, 0)
sol = solve(eq, ivp={u: [1.0, 0.0], u.d(): [0.0, 1.0]}, tspan=(0, 10))

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

Boundary value problems

from odecast import BC

y = var("y")
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)

For nonlinear problems an initial guess can be passed, e.g. guess={y: lambda t: t * (1 - t)}.

Symbolic solutions

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

Without conditions the general solution is returned, e.g. C1*sin(t) + C2*cos(t).

Parameters and matrix products

SymPy symbols in the equations are parameters. Solving again with other values reuses the compiled system, so parameter studies are fast:

import numpy as np

u = var("u", shape=2)
k = sp.Symbol("k")
A = np.array([[0, 1], [-k, -0.3]], dtype=object)

for kv in [1.0, 2.0, 4.0]:
    sol = solve(Eq(u.d(), A @ u), ivp={u: [1, 0]}, tspan=(0, 10), params={k: kv})

Differential-algebraic equations

A pendulum in Cartesian coordinates, with the string tension lam as an algebraic variable:

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),
]
sol = solve(eqs, ivp={x: 0.6, x.d(): 0, y: -0.8, y.d(): 0}, tspan=(0, 10))
sol[lam]                      # tension
sol.constraint_residuals()    # remaining violation of x**2 + y**2 = 1

Events

from odecast import Event

y = var("y")
sol = solve(
    Eq(y.d(2), -9.81),
    ivp={y: 10, y.d(): 0},
    tspan=(0, 10),
    events=[Event(Eq(y, 0), terminal=True)],   # stop when y reaches 0
)
sol.t_events[0]   # [1.428]

Equilibria and fitting

from odecast import equilibria, fit

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

c, k = sp.symbols("c k")
result = fit(Eq(y.d(2) + c * y.d() + k * y, 0), {y: (t_data, y_data)},
             guess={c: 1.0, k: 1.0}, ivp={y: 1, y.d(): 0})
result.params, result.stderr

Analysis

from odecast import equilibria, linearize, continuation, periodic_orbit, phase_portrait

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

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

r = sp.Symbol("r")
continuation(Eq(x.d(), r + x - x**3 / 3), r, start={x: -2.5, r: -3},
             p_range=(-3, 3)).bifurcations
# [Bifurcation(fold at param=0.666667, x=-1), Bifurcation(fold at param=-0.666667, x=1)]

v = var("v")
periodic_orbit(Eq(v.d(2) - (1 - v**2) * v.d() + v, 0), ivp={v: 2, v.d(): 0},
               period=6.5, transient=20)
# PeriodicOrbit(period=6.66329, stable)

phase_portrait() draws direction field, nullclines, equilibria and trajectories of planar systems; fit() estimates parameters from data, with exact derivatives from solve(..., sensitivities=[...]).

Delay and stochastic equations

from odecast import solve_sde

y = var("y")
solve(Eq(y.d(), 1.8 * y * (1 - y.delay(1))), ivp={y: 0.5}, tspan=(0, 50))

sol = solve_sde(Eq(y.d(), 0.1 * y), noise={y: 0.3 * y}, ivp={y: 1.0},
                tspan=(0, 1), dt=1e-3, paths=1000, method="milstein")
sol.mean(y), sol.std(y)

Partial differential equations

With pip install "odecast[pde]", time-dependent PDEs are solved with the method of lines, using finite differences from findiff:

from odecast.pde import Grid, field, solve_pde, Neumann

u = field("u", Grid(x=np.linspace(0, 100, 501)))
sol = solve_pde(
    Eq(u.d(), u.dx(2) + u * (1 - u)),          # Fisher-KPP
    ic={u: lambda x: np.where(x < 5, 1.0, 0.0)},
    bc={u: {"x": (Neumann(0), Neumann(0))}},
    tspan=(0, 40),
)

Output

sol.plot() plots the solution with Matplotlib, sol.to_pandas() returns a DataFrame, and in Jupyter equations and symbolic solutions are displayed as formulas.

Solver options

Additional keyword arguments to solve are passed to scipy.integrate.solve_ivp (or solve_bvp for boundary value problems):

solve(eq, ivp=ic, tspan=(0, 10), method="LSODA", rtol=1e-9, t_eval=t_points)

Limitations

odecast is at an early stage (alpha). In particular:

• Boundary value problems must be two-point problems: all conditions lie at the end points of tspan.
• Functions must be applied to components of vector variables (sp.sin(u[0])), not to the vector as a whole.
• Matrix products require one constant operand: A @ u works, a product of two variables such as M @ u does not.
• The SymPy backend can solve what sympy.dsolve can solve: single equations and linear systems, but generally not nonlinear coupled systems.
• Differential-algebraic systems are reduced by differentiating constraints one at a time; structures that need the Pantelides algorithm are not supported. Boundary value problems, the SymPy backend and equilibria() support only algebraic equations that are linear in the algebraic variables and determine them directly (index 1).
• Delay equations must have constant delays and may not delay the highest derivative (neutral equations).
• The analysis functions find equilibria with sympy.solve; for periodic nonlinearities such as sin only representatives are returned.

Examples

The examples/ directory contains runnable scripts, e.g.

python examples/01_ivp_damped_oscillator.py

More examples with plots are in the documentation.

Development

git clone https://github.com/maroba/odecast.git
cd odecast
pip install -e ".[dev]"
pytest
ruff check .
black --check .

See the contributing guide for details. Bug reports and pull requests are welcome.

License

MIT. See LICENSE.