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](VERN7_BENCHMARK_REPORT.md) 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):

```rust
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:

```rust
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.

## 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:

```rust
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. `f64` implements it by default,
and so do the forward-mode dual numbers of [`num-dual`](https://crates.io/crates/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:

```rust,ignore
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:

```rust
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.