Ordinary Differential Equations (ODE)

May 13, 2026 · View on GitHub

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

// 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.

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

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

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

ExampleDescription & Demonstrated Features
Exponential GrowthSolves 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 OscillatorSimulates 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 GrowthModels 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 ModelSimulates 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 PendulumSimulates 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])
IntegrationDemonstrates the differences between even, dense, t_eval, and the default solout methods for a simple differential equation with an easily found analytical solution.
Cr3bpSimulates 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 OscillatorDemonstrates 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 ODESolves a system of ODEs using a matrix system. Demonstrates how to define a system of equations using matrices.
Custom SoloutDemonstrates 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.
SchrodingerSolves the time-dependent Schrödinger equation using the dop853 method. Demonstrates the use of complex numbers in the ODE system.
BrusselatorDemonstrates 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 STMPropagates 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 SensitivitiesDemonstrates Adjoint Sensitivity Analysis (ASA) using backward-integration to compute gradients of a cost function with respect to parameters.
Quadrature via State AugmentationComputes 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 BVPSolves 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 github repository which compares the performance of the differential-equations library with Fortran implementations of the DOP853 method.

ProblemImplementationMean [ms]Min [ms]Max [ms]Relative
Van der Pol Osc.Rust8.8 ± 0.38.210.31.00
Van der Pol Osc.Fortran12.5 ± 2.911.575.11.41 ± 0.33
CR3BPRust7.0 ± 0.46.58.61.00
CR3BPFortran9.8 ± 1.89.146.01.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.