Numra supports seven types of differential equations, each designed for a different
class of dynamical system. This chapter provides a unified comparison to help you
choose the right equation type and solver for your problem.
Type Crate Equation Form Memory Randomness ODE numra-odey ′ ( t ) = f ( t , y ) y'(t) = f(t, y) y ′ ( t ) = f ( t , y ) None None SDE numra-sded X = f d t + g d W dX = f\,dt + g\,dW d X = f d t + g d W None Wiener process DDE numra-ddey ′ ( t ) = f ( t , y ( t ) , y ( t − τ ) ) y'(t) = f(t, y(t), y(t-\tau)) y ′ ( t ) = f ( t , y ( t ) , y ( t − τ )) Finite (discrete delays) None FDE numra-fdeC D α y = f ( t , y ) {}^C D^\alpha y = f(t, y) C D α y = f ( t , y ) Infinite (power-law kernel) None IDE numra-idey ′ = f ( t , y ) + ∫ 0 t K ( t , s , y ) d s y' = f(t,y) + \int_0^t K(t,s,y)\,ds y ′ = f ( t , y ) + ∫ 0 t K ( t , s , y ) d s Infinite (general kernel) None PDE numra-pde∂ u / ∂ t = L [ u ] \partial u/\partial t = \mathcal{L}[u] ∂ u / ∂ t = L [ u ] None None SPDE numra-spde∂ u / ∂ t = L [ u ] + σ ξ \partial u/\partial t = \mathcal{L}[u] + \sigma\,\xi ∂ u / ∂ t = L [ u ] + σ ξ None Space-time noise
Type Time Space Stochastic ODE t t t — — SDE t t t — W ( t ) W(t) W ( t ) DDE t t t — — FDE t t t — — IDE t t t — — PDE t t t x x x (1D, 2D, 3D)— SPDE t t t x x x (1D)ξ ( x , t ) \xi(x,t) ξ ( x , t )
Type What You Must Provide ODE Point value y ( t 0 ) y(t_0) y ( t 0 ) SDE Point value X ( t 0 ) X(t_0) X ( t 0 ) + random seed DDE History function φ ( t ) \varphi(t) φ ( t ) for t ≤ t 0 t \leq t_0 t ≤ t 0 FDE Point value y ( t 0 ) y(t_0) y ( t 0 ) IDE Point value y ( t 0 ) y(t_0) y ( t 0 ) PDE Initial field u ( x , t 0 ) u(x, t_0) u ( x , t 0 ) + boundary conditions SPDE Initial field u ( x , t 0 ) u(x, t_0) u ( x , t 0 ) + boundary conditions + random seed
Type State Space Dimension ODE R n \mathbb{R}^n R n Finite SDE R n \mathbb{R}^n R n Finite DDE Function space C ( [ − τ , 0 ] , R n ) C([-\tau, 0], \mathbb{R}^n) C ([ − τ , 0 ] , R n ) Infinite-dimensional FDE R n \mathbb{R}^n R n (with full history stored internally)Finite (but O(N) storage) IDE R n \mathbb{R}^n R n (with full history stored internally)Finite (but O(N) storage) PDE Function space (discretized to grid) N grid N_{\text{grid}} N grid SPDE Function space (discretized to grid) N grid N_{\text{grid}} N grid
Type Cost Per Step Total Cost (N N N steps) Memory ODE O ( n ) O(n) O ( n ) O ( N n ) O(Nn) O ( N n ) O ( n ) O(n) O ( n ) SDE O ( n ) O(n) O ( n ) O ( N n ) O(Nn) O ( N n ) O ( n ) O(n) O ( n ) DDE O ( n ) O(n) O ( n ) per stage + interpolationO ( N n ) O(Nn) O ( N n ) O ( N n ) O(Nn) O ( N n ) for historyFDE O ( n ⋅ N ) O(n \cdot N) O ( n ⋅ N ) (history sum)O ( N 2 n ) O(N^2 n) O ( N 2 n ) O ( N n ) O(Nn) O ( N n ) IDE (general) O ( n ⋅ N ) O(n \cdot N) O ( n ⋅ N ) (quadrature)O ( N 2 n ) O(N^2 n) O ( N 2 n ) O ( N n ) O(Nn) O ( N n ) IDE (Prony) O ( n m ) O(nm) O ( nm ) (m m m exponential terms)O ( N n m ) O(Nnm) O ( N nm ) O ( n m ) O(nm) O ( nm ) PDE O ( N grid ) O(N_{\text{grid}}) O ( N grid ) per ODE stepO ( N ⋅ N grid ) O(N \cdot N_{\text{grid}}) O ( N ⋅ N grid ) O ( N grid ) O(N_{\text{grid}}) O ( N grid ) SPDE O ( N grid ) O(N_{\text{grid}}) O ( N grid ) (white) to O ( N grid 2 ) O(N_{\text{grid}}^2) O ( N grid 2 ) (colored)O ( N ⋅ N grid ) O(N \cdot N_{\text{grid}}) O ( N ⋅ N grid ) O ( N grid ) O(N_{\text{grid}}) O ( N grid )
The O ( N 2 ) O(N^2) O ( N 2 ) total cost of FDE and IDE solvers is the main computational bottleneck
for long integrations of memory-dependent equations. For IDEs with sum-of-exponential
kernels, the Prony solver reduces this to O ( N ) O(N) O ( N ) .
Problem Size Recommended Approach Small (n ≤ 20 n \leq 20 n ≤ 20 , N ≤ 10 4 N \leq 10^4 N ≤ 1 0 4 ) Any solver works well Medium (n ≤ 100 n \leq 100 n ≤ 100 , N ≤ 10 5 N \leq 10^5 N ≤ 1 0 5 ) Consider implicit ODE solvers for stiff problems Large (n ≤ 1000 n \leq 1000 n ≤ 1000 , N ≤ 10 6 N \leq 10^6 N ≤ 1 0 6 ) Avoid general IDE/FDE solvers; use Prony or ODE reduction PDE (N grid ≤ 200 N_{\text{grid}} \leq 200 N grid ≤ 200 ) MOL with explicit or implicit ODE solver PDE (N grid > 200 N_{\text{grid}} > 200 N grid > 200 ) MOL with implicit solver (Radau5, Esdirk, Bdf)
Solver Type Order Step Size Best For DoPri5Explicit RK 5(4) Adaptive General purpose Tsit5Explicit RK 5(4) Adaptive Efficient (FSAL) Vern6Explicit RK 6(5) Adaptive High accuracy Vern7Explicit RK 7(6) Adaptive High accuracy Vern8Explicit RK 8(7) Adaptive Very high accuracy Radau5Implicit RK 5 Adaptive Stiff problems Esdirk32ESDIRK 2(1) Adaptive Mildly stiff Esdirk43ESDIRK 3(2) Adaptive Moderately stiff Esdirk54ESDIRK 4(3) Adaptive Stiff, high accuracy BdfMultistep 1—5 Adaptive Stiff, variable order AutoAutomatic Varies Adaptive Stiffness detection
Solver Strong Order Weak Order Step Size Best For EulerMaruyama0.5 1.0 Fixed Quick prototyping Milstein1.0 1.0 Fixed Multiplicative noise Sra11.0—1.5 — Adaptive Additive noise, accuracy Sra2— 2.0 Adaptive Monte Carlo statistics
Solver Type Order Step Size Features MethodOfStepsExplicit RK (DoPri5) 5(4) Adaptive Hermite interpolation, discontinuity tracking, state-dependent delays
Solver Type Convergence Step Size Features L1SolverL1 scheme O ( Δ t 2 − α ) O(\Delta t^{2-\alpha}) O ( Δ t 2 − α ) Fixed Caputo derivative, fixed-point iteration
Solver Time Stepping Quadrature Cost Best For VolterraSolverEuler Trapezoidal O ( N 2 ) O(N^2) O ( N 2 ) Simple problems VolterraRK4SolverRK4 Trapezoidal O ( N 2 ) O(N^2) O ( N 2 ) General kernels PronySolverAugmented ODE Recursive O ( N ) O(N) O ( N ) Sum-of-exponential kernels
Approach Spatial Discretization Time Integration Best For MOLSystemFinite differences (2nd order) Any ODE solver Parabolic, reaction-diffusion MOLSystem2DFinite differences (2nd order) Any ODE solver 2D problems
Solver Method Noise Types Step Size Features MolSdeSolverMOL + EM/Milstein White, Colored, TraceClass Fixed or Adaptive Ensemble support
All equation types follow a common pattern: define a system struct implementing a
trait, configure options with a builder, call a solver, and inspect the result.
fn rhs ( & self , t : S , y : & [ S ], dydt : &mut [ S ]);
fn drift ( & self , t : S , x : & [ S ], f : &mut [ S ]);
fn diffusion ( & self , t : S , x : & [ S ], g : &mut [ S ]);
fn noise_type ( & self ) -> NoiseType ;
// DDE: y'(t) = f(t, y(t), y(t-tau))
fn delays ( & self ) -> Vec < S >;
fn rhs ( & self , t : S , y : & [ S ], y_delayed : & [ & [ S ]], dydt : &mut [ S ]);
// FDE: D^alpha y = f(t, y)
fn rhs ( & self , t : S , y : & [ S ], f : &mut [ S ]);
// IDE: y' = f(t,y) + integral K(t,s,y(s)) ds
fn rhs ( & self , t : S , y : & [ S ], f : &mut [ S ]);
fn kernel ( & self , t : S , s : S , y_s : & [ S ], k : &mut [ S ]);
// PDE: du/dt = F(t, x, u, du/dx, d2u/dx2)
fn rhs ( & self , t : S , x : S , u : S , du_dx : S , d2u_dx2 : S ) -> S ;
// SPDE: du/dt = L[u] + sigma(u) * noise
fn drift ( & self , t : S , u : & [ S ], du : &mut [ S ], grid : & Grid1D < S >);
fn diffusion ( & self , t : S , u : & [ S ], sigma : &mut [ S ], grid : & Grid1D < S >);
fn noise_correlation ( & self ) -> NoiseCorrelation < S >;
Type Options Struct Key Parameters ODE SolverOptionsrtol, atol, h0, h_max, denseSDE SdeOptionsdt, rtol, atol, seed, save_trajectoryDDE DdeOptionsrtol, atol, h0, track_discontinuities, discontinuity_orderFDE FdeOptionsdt, max_steps, tol, max_iterIDE IdeOptionsdt, max_steps, tol, quad_pointsPDE SolverOptions (via ODE)Same as ODE (MOL system is an ODE system) SPDE SpdeOptionsdt, n_output, method, seed, adaptive
All result types provide a consistent interface:
// Common across all result types:
result . y_final () // Final state vector
result . y_at (i) // State at index i
result . t // Time points (Vec)
result . stats // Solver statistics
Type Result Struct Storage Format Extra Fields ODE SolverResultRow-major flat Vec<S> n_accept, n_reject, n_evalSDE SdeResultRow-major flat Vec<S> n_drift, n_diffusion, n_accept, n_rejectDDE DdeResultRow-major flat Vec<S> n_eval, n_accept, n_reject, n_discontinuitiesFDE FdeResultRow-major flat Vec<S> n_rhs, n_stepsIDE IdeResultRow-major flat Vec<S> n_rhs, n_kernel, n_stepsPDE SolverResult (via ODE)Row-major flat Vec<S> Same as ODE SPDE SpdeResultVec<Vec<S>> per time stepn_steps, n_drift, n_diffusion, n_reject
Start here: What does the derivative depend on?
Only the current state y ( t ) y(t) y ( t ) : ODE
Current state + random noise : SDE
Current state + state at past time(s) y ( t − τ ) y(t - \tau) y ( t − τ ) : DDE
Fractional-order derivative (power-law memory): FDE
Current state + integral over history (general kernel): IDE
Current state + spatial derivatives ∂ 2 u / ∂ x 2 \partial^2 u / \partial x^2 ∂ 2 u / ∂ x 2 : PDE
Spatial derivatives + random forcing : SPDE
Type Non-Stiff Solver Stiff Solver ODE DoPri5, Tsit5, Vern6/7/8Radau5, Esdirk, BdfSDE EulerMaruyama— (no implicit SDE solver currently) DDE MethodOfSteps— (explicit only currently) FDE L1Solver (implicit fixed-point)L1Solver handles mild stiffnessIDE VolterraRK4Solver— (explicit only currently) PDE (via MOL) DoPri5 (coarse grid)Radau5, Bdf (fine grid)SPDE EulerMaruyama— (explicit only currently)
If your system has memory, the type of memory determines the equation:
Memory Type DE Type Example No memory ODE y ′ = − k y y' = -ky y ′ = − k y Discrete past states DDE y ′ ( t ) = f ( y ( t ) , y ( t − τ ) ) y'(t) = f(y(t), y(t-\tau)) y ′ ( t ) = f ( y ( t ) , y ( t − τ )) Power-law (Caputo) kernel FDE C D 0.5 y = f ( y ) {}^C D^{0.5} y = f(y) C D 0.5 y = f ( y ) Exponential kernel IDE (Prony) y ′ = f ( y ) + ∫ e − b ( t − s ) y ( s ) d s y' = f(y) + \int e^{-b(t-s)} y(s)\,ds y ′ = f ( y ) + ∫ e − b ( t − s ) y ( s ) d s General kernel IDE (Volterra) y ′ = f ( y ) + ∫ K ( t , s , y ( s ) ) d s y' = f(y) + \int K(t,s,y(s))\,ds y ′ = f ( y ) + ∫ K ( t , s , y ( s )) d s
If your system has spatial structure:
Spatial Structure DE Type Approach No spatial dependence ODE, SDE, DDE, FDE, or IDE Direct Spatial derivatives (deterministic) PDE MOL -> ODE solver Spatial derivatives + noise SPDE MOL -> SDE solver
Several equation types are closely related and can sometimes be converted:
A fractional differential equation C D α y = f ( t , y ) {}^C D^\alpha y = f(t,y) C D α y = f ( t , y ) is equivalent to a
Volterra IDE with a power-law kernel:
y ′ ( t ) = f ( t , y ) + 1 Γ ( 1 − α ) ∫ 0 t ( t − s ) − α y ′ ( s ) d s y'(t) = f(t,y) + \frac{1}{\Gamma(1-\alpha)} \int_0^t (t-s)^{-\alpha} y'(s)\,ds y ′ ( t ) = f ( t , y ) + Γ ( 1 − α ) 1 ∫ 0 t ( t − s ) − α y ′ ( s ) d s
The dedicated FDE solver (L1Solver) uses optimized Caputo weights that are more
efficient than the general Volterra quadrature for this specific kernel.
Through the Method of Lines, any PDE becomes a large ODE system:
d u i d t = F ( u i − 1 , u i , u i + 1 ) , i = 1 , … , N − 1 \frac{du_i}{dt} = F(u_{i-1}, u_i, u_{i+1}), \quad i = 1, \ldots, N-1 d t d u i = F ( u i − 1 , u i , u i + 1 ) , i = 1 , … , N − 1
Numra’s MOLSystem implements OdeSystem, so all ODE solvers work directly.
Through the Method of Lines, an SPDE becomes a system of SDEs (one per grid point).
Numra’s MolSdeWrapper implements SdeSystem, bridging the two crates.
When the IDE kernel is a sum of exponentials K ( τ ) = ∑ a i e − b i τ K(\tau) = \sum a_i e^{-b_i \tau} K ( τ ) = ∑ a i e − b i τ ,
the integral satisfies an auxiliary ODE, so the IDE becomes an augmented ODE system
of dimension n + n m n + nm n + nm where m m m is the number of exponential terms.
When the diffusion coefficient g = 0 g = 0 g = 0 , an SDE reduces to a deterministic ODE.
Use numra-ode for better accuracy and efficiency in this case.
When all delays τ i = 0 \tau_i = 0 τ i = 0 , a DDE reduces to an ODE. Use numra-ode for
better performance.
Consider modeling a system with exponential memory decay. This can be formulated as
three different equation types:
use numra_ide :: { IdeSystem , VolterraSolver , IdeSolver , IdeOptions };
impl IdeSystem < f64 > for MemorySystem {
fn dim ( & self ) -> usize { 1 }
fn rhs ( & self , _t : f64 , y : & [ f64 ], f : &mut [ f64 ]) {
fn kernel ( & self , t : f64 , s : f64 , y_s : & [ f64 ], k : &mut [ f64 ]) {
k[ 0 ] = ( - (t - s)) . exp () * y_s[ 0 ];
let opts = IdeOptions :: default () . dt ( 0.01 );
let result = VolterraSolver :: solve ( & MemorySystem , 0.0 , 5.0 , & [ 1.0 ], & opts) ? ;
use numra_ide :: { IdeSystem , PronySolver , IdeSolver , IdeOptions };
// Same system, but recognized as sum-of-exponentials
// PronySolver reduces the integral to an auxiliary ODE
let opts = IdeOptions :: default () . dt ( 0.01 );
let result = PronySolver :: solve ( & MemorySystem , 0.0 , 5.0 , & [ 1.0 ], & opts) ? ;
use numra_ode :: { DoPri5 , Solver , SolverOptions };
// Manually augment: y' = -y + I, I' = y - I
// where I(t) = integral_0^t exp(-(t-s)) y(s) ds
// ... (implement OdeSystem with dim = 2)
// Cost: O(N), and benefits from adaptive high-order ODE solvers
Feature ODE SDE DDE FDE IDE PDE SPDE Crate numra-odenumra-sdenumra-ddenumra-fdenumra-idenumra-pdenumra-spdeSpatial No No No No No Yes Yes Stochastic No Yes No No No No Yes Memory No No Discrete Power-law General No No Adaptive dt Yes Yes (SRA) Yes No No Yes (via ODE) Yes (step-doubling) Implicit Yes No No Yes (FP) No Yes (via ODE) No Ensemble No Yes No No No No Yes Total cost O ( N ) O(N) O ( N ) O ( N ) O(N) O ( N ) O ( N ) O(N) O ( N ) O ( N 2 ) O(N^2) O ( N 2 ) O ( N 2 ) O(N^2) O ( N 2 ) O ( N ) O(N) O ( N ) O ( N ) O(N) O ( N )