viva_math/ode

Numerical ODE / SDE solvers.

Single-step integrators for scalar systems dx/dt = f(t, x) and stochastic systems dx = f(t, x) dt + g(t, x) dW. Step functions are pure: they take the current (t, x) and produce the next state.

Methods

Function	Order	Use case
`euler`	1	Quick, stiff systems with small dt
`rk2_midpoint`	2	Mid-accuracy, half the cost of RK4
`rk2_heun`	2	Trapezoidal/improved Euler
`rk4`	4	Workhorse for non-stiff systems
`euler_maruyama`	0.5*	SDE with additive noise (* strong order)
`milstein`	1.0*	SDE with multiplicative noise

For full trajectories use integrate / integrate_sde.

Types

Diffusion

</>

Diffusion function g(t, x) for SDEs.

pub type Diffusion =
  fn(Float, Float) -> Float

Drift

</>

Drift function f(t, x).

pub type Drift =
  fn(Float, Float) -> Float

ForceLaw

</>

Force law F(q): acceleration as a function of position (unit mass).

pub type ForceLaw =
  fn(Float) -> Float

Values

dop54

</>

pub fn dop54(
  f: fn(Float, Float) -> Float,
  t: Float,
  x: Float,
  dt: Float,
) -> #(Float, Float)

One step of Dormand-Prince 5(4) with embedded error estimate.

Returns #(x_new, err) where x_new is the 5th-order solution and err is the absolute difference from the embedded 4th-order estimate, suitable for step-size control. Per-step truncation error is O(dt⁶).

euler

</>

pub fn euler(
  f: fn(Float, Float) -> Float,
  t: Float,
  x: Float,
  dt: Float,
) -> Float

Euler step: xₙ₊₁ = xₙ + dt · f(tₙ, xₙ).

euler_maruyama

</>

pub fn euler_maruyama(
  f: fn(Float, Float) -> Float,
  g: fn(Float, Float) -> Float,
  t: Float,
  x: Float,
  dt: Float,
  seed: random.Seed,
) -> #(Float, random.Seed)

Euler-Maruyama step for dx = f dt + g dW with normal increment dW = √dt·Z.

Returns the new state and the advanced PRNG seed.

integrate

</>

pub fn integrate(
  method: fn(fn(Float, Float) -> Float, Float, Float, Float) -> Float,
  f: fn(Float, Float) -> Float,
  t0: Float,
  x0: Float,
  dt: Float,
  steps: Int,
) -> List(#(Float, Float))

Integrate a deterministic ODE with a fixed-step method.

integrate_sde

</>

pub fn integrate_sde(
  f: fn(Float, Float) -> Float,
  g: fn(Float, Float) -> Float,
  t0: Float,
  x0: Float,
  dt: Float,
  steps: Int,
  seed: random.Seed,
) -> #(List(#(Float, Float)), random.Seed)

Integrate a stochastic system with Euler-Maruyama.

integrate_symplectic

</>

pub fn integrate_symplectic(
  method: fn(fn(Float) -> Float, Float, Float, Float) -> #(
    Float,
    Float,
  ),
  force: fn(Float) -> Float,
  t0: Float,
  q0: Float,
  v0: Float,
  dt: Float,
  steps: Int,
) -> List(#(Float, Float, Float))

Integrate a Hamiltonian system with a symplectic method.

leapfrog

</>

pub fn leapfrog(
  force: fn(Float) -> Float,
  q: Float,
  v: Float,
  dt: Float,
) -> #(Float, Float)

Leapfrog (kick-drift-kick) — order 2, symplectic.

milstein

</>

pub fn milstein(
  f: fn(Float, Float) -> Float,
  g: fn(Float, Float) -> Float,
  dg_dx: fn(Float, Float) -> Float,
  t: Float,
  x: Float,
  dt: Float,
  seed: random.Seed,
) -> #(Float, random.Seed)

Milstein step — corrects Euler-Maruyama with the Itô derivative of g.

dx = f dt + g dW + ½g·g’(t,x)·(dW² - dt)

dg_dx is the partial derivative ∂g/∂x. If g doesn’t depend on x (additive noise), use euler_maruyama instead — they coincide.

position_verlet

</>

pub fn position_verlet(
  force: fn(Float) -> Float,
  q_prev: Float,
  q: Float,
  dt: Float,
) -> Float

Position Verlet — q(t+dt) = 2·q(t) - q(t-dt) + a(q(t))·dt².

rk2_heun

</>

pub fn rk2_heun(
  f: fn(Float, Float) -> Float,
  t: Float,
  x: Float,
  dt: Float,
) -> Float

Heun (improved Euler) RK2: trapezoidal predictor-corrector.

rk2_midpoint

</>

pub fn rk2_midpoint(
  f: fn(Float, Float) -> Float,
  t: Float,
  x: Float,
  dt: Float,
) -> Float

Midpoint Runge-Kutta (RK2): evaluates f at the half-step midpoint.

rk4

</>

pub fn rk4(
  f: fn(Float, Float) -> Float,
  t: Float,
  x: Float,
  dt: Float,
) -> Float

Classical 4th-order Runge-Kutta.

rkf45

</>

pub fn rkf45(
  f: fn(Float, Float) -> Float,
  t: Float,
  x: Float,
  dt: Float,
) -> #(Float, Float)

Deprecated alias retained for backwards compatibility. Returns the same pair as dop54. To be removed once external callers migrate. One step of RKF45 with an error estimate.

Returns #(x5, error) where x5 is the 5th-order estimate and error is |x5 - x4|, useful for step-size control.

velocity_verlet

</>

pub fn velocity_verlet(
  force: fn(Float) -> Float,
  q: Float,
  v: Float,
  dt: Float,
) -> #(Float, Float)

Velocity Verlet — order 2, symplectic, time-reversible.

yoshida4

</>

pub fn yoshida4(
  force: fn(Float) -> Float,
  q: Float,
  v: Float,
  dt: Float,
) -> #(Float, Float)

Yoshida 4th-order symplectic integrator.

Composes three leapfrog steps with Yoshida (1990) coefficients. Order 4 global error, ~3× the cost of leapfrog per step but allows ~10× larger dt for the same accuracy.