A library, written in Rust, for integrating ordinary differential equations. For now, this is relatively simple, but it does have key features that are needed for orbit propagation, ray tracing, and field line tracing:

Features

Explicit Runge-Kutta Methods (Non-Stiff Problems)

Method	Order	Stages	Dense Output	Best Use Case
BS3 (Bogacki-Shampine)	3(2)	4	3rd order	Moderate accuracy (rtol ~ 1e-4 to 1e-6)
DormandPrince45	5(4)	7	4th order	General purpose (rtol ~ 1e-6 to 1e-8)
Vern7 (Verner)	7(6)	10+6	7th order	High accuracy (rtol ~ 1e-8 to 1e-12)

Performance at 1e-10 tolerance:

Vern7: 1.1-4.4x faster than DP5
Vern7: 53-83x faster than BS3

See benchmark report for detailed performance analysis.

Other Features

Adaptive time stepping with PI controller
Callback events with zero-crossing detection
Dense output interpolation at any time point
Parameters in derivative and callback functions
Lazy computation of extra interpolation stages (Vern7)
PID step size control as an alternative to the PI controller
Automatic Differentiation, the same problem can be solved with floats or also differentiated with dual numbers

Future Improvements

More algorithms
- Rosenbrock methods (for stiff problems)
- Tsit5
- Runge-Kutta Cash-Karp
Composite Algorithms
Automatic Stiffness Detection
Fixed Time Steps
Boolean callback eventing
Improved solution handling like DifferentialEquations.jl

To Use

For now, here is a simple example of using the propagator to solve a simple second-order system (the pendulum problem):

use nalgebra::Vector2;
use ordinary_diffeq::prelude::*;
use std::f64::consts::PI;

// Define the system (parameters, derivative, and initial state)
type Params = (f64, f64); // Gravity and Length of Pendulum
let params = (9.81, 1.0);

fn derivative(_t: f64, y: Vector2<f64>, p: &Params) -> Vector2<f64> {
    let &(g, l) = p;
    let theta = y[0];
    let d_theta = y[1];
    Vector2::new( d_theta, -(g/l) * theta.sin() )
}

let y0 = Vector2::new(0.0, PI/2.0);

// Set up the problem (ODE, Integrator, Controller, and Callbacks)
let ode = ODE::new(&derivative, 0.0, 6.3, y0, params);
let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-6);
let controller = PIController::default();

let value_too_high = FnCallback {
    event: &|t: f64, _y: &Vector2<f64>, _p: &Params| { 5.0 - t },
    effect: &stop,
};

// Solve the problem
let mut problem = Problem::new(ode, dp45, controller).with_callback(&value_too_high);
let solution = problem.solve();

// Can interpolate solutions to whatever you want
let _interpolated_answer = solution.interpolate(4.4);

Automatic differentiation

State, time and step size are generic over the Scalar trait. f64 implements it by default, and so do the forward-mode dual numbers of num-dual (the autodiff feature, on by default). Solving at a dual scalar propagates exact derivatives alongside the state.

The only change needed is to the signature of the right-hand side. Here is the pendulum from above, before and after:

fn derivative(_t: f64, y: Vector2<f64>, p: &Params) -> Vector2<f64> {
    let &(g, l) = p;
    Vector2::new( y[1], -(g/l) * y[0].sin() )
}

fn derivative<T: Scalar>(_t: T, y: Vector2<T>, p: &Params) -> Vector2<T> {
    let &(g, l) = p;
    Vector2::new( y[1], -y[0].sin() * (g/l) )
}

Params stays (f64, f64) — parameters only need to become T themselves if you want to differentiate with respect to them. Real constants mix into T arithmetic directly; Scalar::lift is only needed to scale a vector by a real, or to build a T out of nothing.

With that one change, the same function can be solved or differentiated:

use nalgebra::{Matrix2, Vector2};
use num_dual::DualSVec64;
use ordinary_diffeq::prelude::*;
use std::f64::consts::PI;

type Params = (f64, f64);

fn derivative<T: Scalar>(_t: T, y: Vector2<T>, p: &Params) -> Vector2<T> {
    let &(g, l) = p;
    Vector2::new(y[1], -y[0].sin() * (g / l))
}

// One propagation in `DualSVec64<2>` yields the whole 2x2 Jacobian of the
// final state with respect to the initial one.
let (y_end, jacobian): (Vector2<f64>, Matrix2<f64>) = num_dual::jacobian(
    |y0: Vector2<DualSVec64<2>>| {
        type T = DualSVec64<2>;
        let ode = ODE::new(&derivative, T::lift(0.0), T::lift(1.0), y0, (9.81, 1.0));
        let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
        let mut problem = Problem::new(ode, vern7, PIController::default());
        let solution = problem.solve();
        *solution.states.last().unwrap()
    },
    &Vector2::new(0.0, PI / 2.0),
);

// The pendulum is Hamiltonian, so its flow is area-preserving
assert!((jacobian.determinant() - 1.0).abs() < 1e-9);

These are exact derivatives of the scheme as it ran, at machine precision, with no finite differencing and no step size to tune.