ode.hpp

Go to the documentation of this file.

/// @file ode/ode.hpp
/// @brief ODE time integrators: Euler, RK4, adaptive RK45 (Dormand-Prince),
///        and symplectic Velocity Verlet / Yoshida-4 for Hamiltonian systems.
///
/// General ODE solvers (for dy/dt = f(t, y)):
///   ode_euler   -- forward Euler, O(h)
///   ode_rk4     -- classic 4th-order Runge-Kutta, O(h⁴)
///   ode_rk45    -- adaptive Dormand-Prince with PI step control, O(h⁵)
///
/// Symplectic integrators (for separable H(q,p) = T(p) + V(q)):
///   ode_verlet   -- velocity Verlet, 2nd-order symplectic, 1 force eval/step
///   ode_yoshida4 -- Yoshida (1990) 4th-order symplectic, 3 force evals/step
///
/// All integrators accept an optional on_step callback invoked after each
/// accepted step, enabling trajectory recording and event detection without
/// storing the full history in the integrator itself.
#pragma once
 
#include "core/types.hpp"
#include "core/vector.hpp"
#include <functional>
 
namespace num {
 
// ── Callable types ────────────────────────────────────────────────────────────
 
/// Right-hand side callable: fills dydt = f(t, y) in-place.
using ODERhsFn = std::function<void(real t, const Vector& y, Vector& dydt)>;
 
/// Called after each accepted step. Use to record trajectories or detect events.
/// @param t  Time at end of step
/// @param y  State at end of step
using StepCallback = std::function<void(real t, const Vector& y)>;
 
/// Acceleration function for symplectic integrators.
/// Fills acc = −∇V(q) / m (generalised force per unit mass) from positions q.
using AccelFn = std::function<void(const Vector& q, Vector& acc)>;
 
/// Per-step callback for symplectic integrators.
/// @param t  Time at end of step
/// @param q  Generalised positions
/// @param v  Generalised velocities
using SymplecticCallback = std::function<void(real t, const Vector& q, const Vector& v)>;
 
// ── Result types ──────────────────────────────────────────────────────────────
 
/// Result returned by general ODE integrators.
struct ODEResult {
    Vector y;          ///< Final state vector
    real   t;          ///< Final time
    idx    steps;      ///< Number of (accepted) steps taken
    bool   converged;  ///< Always true for fixed-step; false if rk45 hit max_steps
};
 
/// Result returned by symplectic integrators.
struct SymplecticResult {
    Vector q;     ///< Final generalised positions
    Vector v;     ///< Final generalised velocities
    real   t;     ///< Final time
    idx    steps; ///< Number of steps taken
};
 
// ── General ODE solvers ───────────────────────────────────────────────────────
 
/// Forward Euler (1st order, fixed step):
///   y_{n+1} = y_n + h · f(t_n, y_n)
ODEResult ode_euler(ODERhsFn f, Vector y0, real t0, real t1, real h,
                    StepCallback on_step = nullptr);
 
/// Classic 4th-order Runge-Kutta (fixed step).
ODEResult ode_rk4(ODERhsFn f, Vector y0, real t0, real t1, real h,
                  StepCallback on_step = nullptr);
 
/// Adaptive Dormand-Prince RK45 with FSAL and PI step-size control.
///
/// Advances from t0 to t1, adjusting h so that the mixed
/// absolute/relative error estimate stays below tolerance:
///   err / (atol + |y|·rtol) ≤ 1
///
/// @param rtol      Relative tolerance (default 1e-6)
/// @param atol      Absolute tolerance (default 1e-9)
/// @param h0        Initial step-size hint (default 1e-3)
/// @param max_steps Hard limit on accepted steps (default 1e6)
/// @param on_step   Optional callback after each accepted step
ODEResult ode_rk45(ODERhsFn f, Vector y0, real t0, real t1,
                   real rtol = 1e-6, real atol = 1e-9,
                   real h0 = 1e-3, idx max_steps = 1000000,
                   StepCallback on_step = nullptr);
 
// ── Symplectic integrators ────────────────────────────────────────────────────
//
// For separable Hamiltonians H(q,p) = T(p) + V(q) only.
// Symplectic (volume-preserving in phase space) and time-reversible —
// energy error is bounded, not growing, over exponentially long runs.
// Compare with RK4 which has secular O(h^4) energy drift per step.
//
// Both methods take q and v separately (not concatenated) and an AccelFn.
 
/// Velocity Verlet — 2nd-order symplectic, 1 force evaluation per step.
///
/// Algorithm (kick-drift-kick):
///   a_n     = accel(q_n)
///   q_{n+1} = q_n + h·v_n + h²/2 · a_n
///   a_{n+1} = accel(q_{n+1})
///   v_{n+1} = v_n + h/2 · (a_n + a_{n+1})
///
/// Acceleration is cached between steps (FSAL — first same as last).
SymplecticResult ode_verlet(AccelFn accel, Vector q0, Vector v0,
                             real t0, real t1, real h,
                             SymplecticCallback on_step = nullptr);
 
/// Yoshida 4th-order symplectic — 3 force evaluations per step.
///
/// Composes three leapfrog sub-steps with Yoshida (1990) coefficients:
///   w₁ = 1 / (2 − ∛2),  w₀ = 1 − 2w₁
///   c₁=c₄=w₁/2,  c₂=c₃=(w₀+w₁)/2,  d₁=d₃=w₁,  d₂=w₀
///
/// Drift-kick sequence per step:
///   q += c₁h·v  →  v += d₁h·a(q)  →  q += c₂h·v  →  v += d₂h·a(q)
///   →  q += c₃h·v  →  v += d₃h·a(q)  →  q += c₄h·v
///
/// Preferred over Verlet when accuracy matters and force is cheap.
SymplecticResult ode_yoshida4(AccelFn accel, Vector q0, Vector v0,
                               real t0, real t1, real h,
                               SymplecticCallback on_step = nullptr);
 
} // namespace num