# Ordinary Differential Equations (ODE)
The `ode` module provides tools for solving ordinary differential equations (ODEs),
including initial value problems (ODE IVPs) and boundary value problems (ODE BVPs).
## Table of Contents
- [Defining a ODE](#defining-an-ode)
- [Solving an Initial Value Problem (ODE IVP builder)](#solving-an-initial-value-problem)
- [Numerical Quadrature](#numerical-quadrature)
- [Examples](#examples)
- [Benchmarks](#benchmarks)
- [Notation](#notation)
## Defining an ODE
The `ODE` trait defines the differential equation `dydt = f(t, y)` for the solver. The differential equation is used to solve the ordinary differential equation.
### ODE Trait
* `diff` - Differential Equation `dydt = f(t, y)` in the form `f(t, &y, &mut dydt)`.
### Event Trait
For event detection with precise zero-crossing detection, implement the separate `Event` trait:
* `config` - Configure event detection parameters (direction filtering, termination count)
* `event` - Event function `g(t,y)` whose zero crossings are detected using Brent-Dekker root finding
### Solout Trait
* `solout` - function to choose which points to save in the solution. This is useful when you want to save only points that meet certain criteria. Common implementations are included in the `solout` module. The `ODE IVP builder` trait implements methods to use them easily without direct interaction as well e.g. `even`, `dense`, and `t_eval`.
### Implementation
```rust
// Includes required elements and common methods.
// Less common methods are in the `methods` module
use differential_equations::prelude::*;
struct LogisticGrowth {
k: f64,
m: f64,
}
impl ODE for LogisticGrowth {
fn diff(&self, _t: f64, y: &f64, dydt: &mut f64) {
*dydt = self.k * *y * (1.0 - *y / self.m);
}
}
// Separate event detection implementation
impl Event for LogisticGrowth {
fn config(&self) -> EventConfig {
EventConfig::new(CrossingDirection::Positive, Some(1)) // Terminate after first event
}
fn event(&self, _t: f64, y: &f64) -> f64 {
// Event function g(t,y) = y - 0.9*m
// Zero crossing occurs when y = 0.9*m
*y - 0.9 * self.m
}
}
```
Note that the `ODE` is defined with generics `<T, Y>` where `T` is the float type
(e.g. `f64` or `f32`) and `Y` is the state vector. `Y` can be a scalar (`f64`, `f32`),
a fixed-size array (`[T; N]`), a `Vec<T>`, or various matrix types from `nalgebra`,
`ndarray`, or `faer`.
## Solving an Initial Value Problem
The `IVP::ode` builder is used to solve the system using the solver. The builder includes methods to set the initial conditions, output strategy, event handling, and numerical method. The `solve` method returns a `Result<Solution, Error>` where `Solution` contains output times, states, status, statistics, and solve time.
```rust
fn main() {
let method = ExplicitRungeKutta::dop853().rtol(1e-12).atol(1e-12);
let y0 = 1.0;
let t0 = 0.0;
let tf = 10.0;
let ode = LogisticGrowth { k: 1.0, m: 10.0 };
let logistic_growth_problem = IVP::ode(&ode, t0, tf, y0);
match logistic_growth_problem
.even(1.0) // uses EvenSolout to save with dt of 1.0
.event(&ode) // Add event detection using the Event trait implementation
.method(method)
.solve() // Solve the system and return the solution
{
Ok(solution) => {
// Check if the solver stopped due to event detection
if let Status::Interrupted = solution.status {
println!("Solver stopped due to event detection");
}
// Print the solution
println!("Solution:");
for (t, y) in solution.iter() {
println!("({:.4}, {:.4})", t, y);
}
// Print the statistics
println!("Function evaluations: {}", solution.evals.function);
println!("Steps: {}", solution.steps.total());
println!("Rejected Steps: {}", solution.steps.rejected);
println!("Accepted Steps: {}", solution.steps.accepted);
}
Err(e) => panic!("Error: {:?}", e),
};
}
```
### Output
```sh
Solver stopped: Reached 90% of carrying capacity
Solution:
(0.0000, 1.0000)
(1.0000, 2.3197)
(2.0000, 4.5085)
(3.0000, 6.9057)
(4.0000, 8.5849)
(4.3944, 9.0000)
Function evaluations: 325
Steps: 22
Rejected Steps: 2
Accepted Steps: 20
```
## Numerical Quadrature
To compute a numerical integral alongside an ODE, encode the quadrature as an additional equation in your system. If your ODE is `dy/dt = f(t, y)` and you want to compute `Q = integral of g(t, y) dt`, augment the state vector to `[y, Q]` and add the equation `dQ/dt = g(t, y)`:
```rust
use nalgebra::{SVector, vector};
struct MyODE;
impl ODE<f64, SVector<f64, 2>> for MyODE {
fn diff(&self, t: f64, y: &SVector<f64, 2>, dydt: &mut SVector<f64, 2>) {
dydt[0] = y[0]; // dy/dt = f(t, y)
dydt[1] = y[0]; // dQ/dt = g(t, y) = y
}
}
```
The initial condition becomes `y0 = vector![y0, Q0]` where `Q0` is typically `0`.
### Controlling Quadrature Influence on Step Size
Per-component tolerances control whether the quadrature affects adaptive step-size selection. The error norm is computed per-component as `|err_i| / (atol_i + rtol_i * max(|y_i|, |y_new_i|))`, and the maximum across all components determines step acceptance.
- **Tight tolerances on Q** (e.g. `1e-12`): the quadrature participates in error control and can force smaller steps to maintain quadrature accuracy.
- **Loose tolerances on Q** (e.g. `1e30`): the quadrature error contribution is negligible, so step size is driven entirely by the primary state `y`. The quadrature is still computed at the solver's full accuracy — it just doesn't force additional step rejections.
```rust
.method(
ExplicitRungeKutta::dop853()
.rtol([1e-12, 1e-12]) // [y tolerance, Q tolerance]
.atol([1e-12, 1e-12])
)
```
This approach requires no changes to the solver — the existing adaptive error estimation, dense output interpolation, and event detection all work on the combined state vector.
## Examples
For more examples, see the `examples` directory. The examples demonstrate different systems, methods, and output methods for different use cases.
| [Exponential Growth](../../examples/ode/01_exponential_growth/main.rs) | Solves a simple exponential growth equation using the `dop853` method. Demonstrates basic usage of `ODE IVP builder` and `ODE` traits. Manually prints results from `Solution` struct fields. |
| [Harmonic Oscillator](../../examples/ode/02_harmonic_oscillator/main.rs) | Simulates a harmonic oscillator system using `rk4` method. Uses a condensed setup to demonstrate chaining to solve without intermediate variables. Uses `last` method on solution to conveniently get results and print. |
| [Logistic Growth](../../examples/ode/03_logistic_growth/main.rs) | Models logistic growth with a carrying capacity. Demonstrates the use of the `event` function to stop the solver based on a condition. In addition shows the use of `even` output for `ODE IVP builder` setup and `iter` method on the solution for output. |
| [SIR Model](../../examples/ode/04_sir_model/main.rs) | Simulates the SIR model for infectious diseases. Uses the `AdamsPredictorCorrector::v4()` method to solve the system. Uses struct as the `State` via `derive(State)` and a custom event termination enum. |
| [Damped Pendulum](../../examples/ode/05_damped_pendulum/main.rs) | Simulates a simple pendulum using the `rkf45` solver. Shows the use of `problem.t_eval()` to define specific points to be saved e.g. `t_eval(vec![1.0, 3.0, 7.5, 10.0])` |
| [Integration](../../examples/ode/06_integration/main.rs) | Demonstrates the differences between `even`, `dense`, `t_eval`, and the default solout methods for a simple differential equation with an easily found analytical solution. |
| [Cr3bp](../../examples/ode/07_cr3bp/main.rs) | Simulates the Circular Restricted Three-Body Problem (CR3BP) using the `dop853` method. Uses the `hyperplane_crossing` method to log when the spacecraft crosses a 3D plane. |
| [Damped Oscillator](../../examples/ode/08_damped_oscillator/main.rs) | Demonstrates the use of the `crossing` method to use the CrossingSolout to log instances where a crossing occurs. In this case, the example saves points where the position is at zero. |
| [Matrix ODE](../../examples/ode/09_matrix_ode/main.rs) | Solves a system of ODEs using a matrix system. Demonstrates how to define a system of equations using matrices. |
| [Custom Solout](../../examples/ode/10_custom_solout/main.rs) | Demonstrates how to create a custom `Solout` implementation to save points based on a custom condition. In addition inside the Solout struct additional calculations are stored each step and accessible after solving is complete. |
| [Schrodinger](../../examples/ode/11_schrodinger/main.rs) | Solves the time-dependent Schrödinger equation using the `dop853` method. Demonstrates the use of complex numbers in the ODE system. |
| [Brusselator](../../examples/ode/12_brusselator/main.rs) | Demonstrates solving a stiff system using implicit Runge-Kutta methods (`gauss_legendre_6`) with analytical jacobian provided to accelerate Newton iterations. |
| [Reduced Two-Body Problem with STM](../../examples/ode/14_r2bp_stm/main.rs) | Propagates a Kepler two-body orbit together with its state transition matrix. Demonstrates dual-number sensitivities, variational equations, and the solver `filter` hook for step-size post-processing. |
| [Adjoint Sensitivities](../../examples/ode/16_adjoint_sensitivities/main.rs) | Demonstrates Adjoint Sensitivity Analysis (ASA) using backward-integration to compute gradients of a cost function with respect to parameters. |
| [Quadrature via State Augmentation](../../examples/ode/17_quadrature/main.rs) | Computes cumulative energy dissipated by a damped oscillator as a quadrature. Demonstrates per-component tolerances — tight on the ODE state, loose on the quadrature. |
| [Pipe Heat Transfer BVP](../../examples/ode/18_pipe_heat_transfer/main.rs) | Solves a steady pipe heat-transfer boundary value problem using `BVP::ode`, the regular `ODE` trait, the `Boundary` residual trait, and `Shooting`. |
## Benchmarks
Test results from [differential-equations-comparison](https://github.com/Ryan-D-Gast/differential-equations-comparison) github repository which compares the performance of the `differential-equations` library with Fortran implementations of the DOP853 method.
| Van der Pol Osc. | Rust | 8.8 ± 0.3 | 8.2 | 10.3 | 1.00 |
| Van der Pol Osc. | Fortran | 12.5 ± 2.9 | 11.5 | 75.1 | 1.41 ± 0.33 |
| CR3BP | Rust | 7.0 ± 0.4 | 6.5 | 8.6 | 1.00 |
| CR3BP | Fortran | 9.8 ± 1.8 | 9.1 | 46.0 | 1.40 ± 0.26 |
Testing has shown that the `differential-equations` Rust implementation is at about 10% to 40% faster than the Fortran implementations of the DOP853 method.
## Notation
Typical ODE libraries either use `x` or `t` for the independent variable and `y` for the dependent variable. This library uses the following notation:
- `t` - The independent variable, typically time often `x` in other ode libraries.
- `y` - The dependent variable, instead of `x` to avoid confusion with an independent variable in other notations.
- `dydt` - The derivative of `y` with respect to `t`.
- `k` - The coefficients of the solver, typically a derivative such as in the Runge-Kutta methods.
Any future methods added to the library should follow this notation to maintain consistency.