ordinary-diffeq/readme.md

ordinary-diffeq

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, whose effects can change the state, change the parameters, or end the solution at the event
• 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
• Runtime integrator selection, the method can be a value, named in a configuration file rather than fixed as a type

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 mut 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(&mut value_too_high);
let solution = problem.solve();

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

Events that change the state

A callback watches a function of the state for a zero crossing. When one is found the step is shortened to land on the crossing, the state there is recorded, and the effect runs with the ODE at that point — so the effect can move the state, and the integration continues from where it leaves it.

That is what it takes to integrate through a discontinuity. A rotating body parameterized by Modified Rodrigues Parameters, for instance, has a singularity at a full turn, avoided by swapping the parameter vector for -sigma / |sigma|^2 — the same attitude, different parameters — whenever |sigma| passes 1:

use nalgebra::{Matrix3, Vector3};
use ordinary_diffeq::prelude::*;

/// MRP kinematics for a constant angular velocity.
fn attitude_rate(_t: f64, s: Vector3<f64>, omega: &Vector3<f64>) -> Vector3<f64> {
    let skew = Matrix3::new(0.0, -s[2], s[1], s[2], 0.0, -s[0], -s[1], s[0], 0.0);
    let b = Matrix3::identity() * (1.0 - s.dot(&s)) + skew * 2.0 + s * s.transpose() * 2.0;
    b * omega * 0.25
}

let mut shadow_set = FnCallback {
    event: &|_t: f64, s: &Vector3<f64>, _p: &Vector3<f64>| 1.0 - s.norm(),
    effect: &|ode: &mut ODE<'_, f64, 3, Vector3<f64>>| {
        ode.y = -ode.y / ode.y.norm_squared();
    },
};

let ode = ODE::new(&attitude_rate, 0.0, 5.0, Vector3::zeros(), Vector3::new(0.0, 0.0, 1.0));
let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
let mut problem = Problem::new(ode, vern7, PIController::default())
    .with_callback(&mut shadow_set);
let solution = problem.solve();

// The parameters never approach the singularity, having swapped at half a turn.
assert!(solution.states.iter().all(|s| s.norm() <= 1.0 + 1e-9));

A jump like that one puts the same time in the solution twice, holding the state on each side of it, and makes the step size controller forget a history that described the solution before the jump. Interpolating at the jump reports the state the effect left, the one the trajectory goes on with. effect takes &mut self, so a callback that implements the trait itself rather than going through FnCallback can record that it has fired. stop ends the solution at the crossing, with the crossing as its last point and the problem's own t_end left alone.

If two events fall inside one step the earlier one fires and the step ends there; the later one is found again on the step after it.

Borrowed or owned

with_callback borrows, which is what lets you read a stateful callback's own fields once the solve is done:

# use nalgebra::Vector1;
# use ordinary_diffeq::prelude::*;
struct Count {
    fired: usize,
}

impl Callback<f64, 1, ()> for Count {
    fn event(&self, t: f64, _y: &Vector1<f64>, _p: &()) -> f64 {
        1.0 - t
    }

    fn effect(&mut self, _ode: &mut ODE<'_, f64, 1, ()>) {
        self.fired += 1;
    }
}

fn derivative(_t: f64, _y: Vector1<f64>, _p: &()) -> Vector1<f64> {
    Vector1::new(1.0)
}

let ode = ODE::new(&derivative, 0.0, 2.0, Vector1::new(0.0), ());
let mut counter = Count { fired: 0 };
let mut problem = Problem::new(ode, DormandPrince45::new(), PIController::default())
    .with_callback(&mut counter);
problem.solve();

assert_eq!(counter.fired, 1);

with_owned_callback takes a Box instead, for a caller that has nothing to lend. A type that wraps this one has to keep the Problem in a field, since solve hands back a solution borrowing from it, and it then cannot also hold a callback in a field and lend it to the problem — that is a self-referential struct, and moving the wrapper would dangle the reference. Moving the callback into the problem sidesteps it:

let problem = Problem::new(ode, integrator, controller)
    .with_owned_callback(Box::new(my_callback));
self.problem = Some(problem);
self.problem.as_mut().unwrap().solve()

An owned callback may hold borrows of its own, so a boxed FnCallback works. Nothing can read its state afterwards; a callback with something to report either stays borrowed, or shares an Rc<Cell<_>> with its owner, which can be read while the solution is still alive.

Choosing the integrator at runtime

The integrator is normally a type parameter, fixed when the code is compiled. Method is an enum over the integrators here that itself implements Integrator, so the choice can be a value instead:

use nalgebra::Vector1;
use ordinary_diffeq::prelude::*;

fn derivative(_t: f64, y: Vector1<f64>, _p: &()) -> Vector1<f64> {
    -y
}

// `MethodName` derives `Serialize`/`Deserialize`, so this can come out of a
// configuration file instead of being written here.
let name: MethodName = serde_json::from_str(r#""Vern7""#).unwrap();
let method: Method<1> = name.instantiate().a_tol(1e-12).r_tol(1e-12);

let ode = ODE::new(&derivative, 0.0, 2.0, Vector1::new(1.0), ());
let mut problem = Problem::new(ode, method, PIController::default());
let solution = problem.solve();

// And the method that ran can be recorded in a form that selects it again
assert_eq!(method.name(), MethodName::Vern7);

Dispatch is a match on the variant, not a trait object, so nothing is erased: the same Method works at a dual scalar and differentiates exactly as the integrator it wraps does. A Vec<Method<D>> holds several methods for comparison over the same problem.

Automatic differentiation

State, time and step size are generic over the Scalar trait, which is a blanket extension of DualNum and DualStruct from num-dual: those two supply the arithmetic, the transcendentals and re, which discards derivative information, and Scalar adds only lift. Both are implemented for f64 as well as for every one of that crate's forward-mode dual numbers, so Scalar covers all of them, and a crate depending on this one alongside another built on num-dual can satisfy both with one bound. 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.