Adding `rkf78` library as a dependency

Run this command in a terminal, in your project's directory:

cargo add rkf78

To add it manually, edit your project's Cargo.toml file and add to the [dependencies] section:

rkf78 = "0.2.1"

The rkf78 library will be automatically available globally. Read the rkf78 library documentation.

Back to the crate overview.

Readme

RKF78

A high-precision Runge-Kutta-Fehlberg 7(8) ODE integrator in Rust.

Why RKF78?

RKF78 is an 8th-order embedded Runge-Kutta method with 7th-order error estimation — one of the highest-order single-step integrators in common use. It takes large, accurate steps on smooth problems, making it ideal when you need high precision without a huge step count.

Zero runtime dependencies — No nalgebra, no BLAS, no math libraries. Just Rust.
Const-generic [T; N] arrays — Compile-time state dimension, zero heap allocation during integration.
Generic over f32/f64 — One codebase, two precision levels.
NASA-heritage coefficients — From the original 1968 Fehlberg paper (NASA TR R-287).
Event detection — Brent's method + Hermite cubic interpolation for precise zero-crossing location, simultaneous multi-event support, Stop/Continue actions.
Dense output — Hermite cubic interpolation between steps for continuous trajectory evaluation.
Step observer — Callback after each accepted step for custom recording or analysis.
GPU batch propagation — Thousands of parallel integrations via wgpu compute shaders.

When to use RKF78

RKF78 excels on smooth, non-stiff ODE systems where high accuracy is needed:

Orbital mechanics and trajectory propagation
Oscillators, wave equations, vibration analysis
Celestial mechanics and N-body problems
Chemical kinetics (non-stiff regimes)
Control system simulation
Any dy/dt = f(t, y) where you need better than 4th-order accuracy

If your problem is stiff (e.g., fast chemical reactions, circuit transients with widely-separated time scales), use an implicit method instead.

Quick Start

use rkf78::{Rkf78, OdeSystem, Tolerances, IntegrationConfig};

// Define your ODE system: y'' + y = 0 (harmonic oscillator)
struct Oscillator { omega: f64 }

impl OdeSystem<f64, 2> for Oscillator {
    fn rhs(&self, _t: f64, y: &[f64; 2], dydt: &mut [f64; 2]) {
        dydt[0] = y[1];
        dydt[1] = -self.omega * self.omega * y[0];
    }
}

let sys = Oscillator { omega: 1.0 };
let tol = Tolerances::new(1e-12, 1e-12);
let mut solver = Rkf78::new(tol);

let y0 = [1.0, 0.0];  // y(0) = 1, y'(0) = 0
let config = IntegrationConfig::new(0.0, 10.0, 0.1);
let (tf, yf) = solver.integrate(&sys, &config, &y0).unwrap();

Event Detection

Detect when a user-defined function crosses zero during integration.

use rkf78::{EventFunction, EventConfig, EventDirection, IntegrationResult, IntegrationConfig};

struct ThresholdCrossing { value: f64 }

impl EventFunction<f64, 2> for ThresholdCrossing {
    fn eval(&self, _t: f64, y: &[f64; 2]) -> f64 {
        y[0] - self.value
    }
}

let event = ThresholdCrossing { value: 0.5 };
let config = EventConfig {
    direction: EventDirection::Falling,
    ..Default::default()
};

let int_config = IntegrationConfig::new(0.0, 10.0, 0.1);
let (result, _collected) = solver.integrate_to_event(&sys, &event, &config, &int_config, &y0).unwrap();
match result {
    IntegrationResult::Event(ev) => println!("Event at t = {:.6}", ev.t),
    IntegrationResult::Completed { t, .. } => println!("No event, reached t = {}", t),
}

For simultaneous events, implement MultiEventFunction<T, N, M> and use integrate_with_multi_events().

Events can also be configured with EventAction::Continue to record all crossings without stopping.

Dense Output

Record the full trajectory and evaluate at arbitrary times:

let (tf, yf, solution) = solver.integrate_dense(&sys, &config, &y0).unwrap();

// Evaluate anywhere in [t0, tf] via Hermite cubic interpolation
let y_mid = solution.eval(5.0).unwrap();
let dy_mid = solution.eval_derivative(5.0).unwrap();

Step Observer

Inspect each accepted step during integration:

use rkf78::StepObserver;

struct MaxErrorTracker(f64);

impl StepObserver<f64, 2> for MaxErrorTracker {
    fn on_step(&mut self, _t: f64, _y: &[f64; 2], _h: f64, error: f64) {
        self.0 = self.0.max(error);
    }
}

let mut tracker = MaxErrorTracker(0.0);
let (tf, yf) = solver.integrate_with_observer(&sys, &config, &y0, &mut tracker).unwrap();

Tolerance Selection

Precision Level	`atol`	`rtol`	Use Case
High	`1e-12`	`1e-12`	Reference solutions, precision analysis
Standard	`1e-10`	`1e-10`	General engineering analysis
Fast	`1e-6`	`1e-6`	Quick surveys, parameter sweeps

For mixed-unit state vectors (e.g., position in km and velocity in km/s), use per-component tolerances via Tolerances::with_components().

Algorithm Details

For a full explanation of the mathematics — Butcher tableau, error estimation, step-size control, and Brent's method — see docs/algorithm.md.

Building and Testing

cargo build            # Build the crate
cargo test             # Run all tests (~78 tests)
cargo test --features gpu  # Include GPU tests (requires GPU)
cargo bench            # Run criterion benchmarks
cargo clippy           # Lint
cargo fmt --check      # Check formatting
cargo run --example harmonic_oscillator  # Run an example

GPU Batch Propagation

With the gpu feature, RKF78 can propagate thousands of trajectories in parallel on the GPU using wgpu compute shaders. The GPU solver uses f32 precision — suitable for Monte Carlo studies, parameter sweeps, and trade studies where throughput matters more than last-digit precision.

Users supply their own WGSL force model, which is prepended to the RKF78 compute shader at pipeline creation time.

use rkf78::gpu::{GpuBatchPropagator, GpuIntegrationParams, GpuState};

let propagator = GpuBatchPropagator::new(force_model_wgsl)?;
let (final_states, statuses) = propagator.propagate_batch(&states, &params, &force_params)?;

See examples/gpu_two_body.rs for a complete example with a Keplerian force model.

Examples

Example	Description
Harmonic Oscillator	Basic integration with exact-solution comparison
Dense Output	Continuous trajectory evaluation via Hermite interpolation
Step Observer	Recording integration progress with a callback
Two-Body Orbit	Keplerian orbit with per-component tolerances
Event Detection	Zero-crossing detection with Stop and Continue actions
GPU Two-Body	GPU batch propagation of multiple orbits

See examples/README.md for a suggested reading order.

References

Fehlberg, E. (1968). "Classical Fifth-, Sixth-, Seventh-, and Eighth-Order Runge-Kutta Formulas with Stepsize Control". NASA TR R-287.
Hairer, E., Norsett, S.P., & Wanner, G. (1993). "Solving Ordinary Differential Equations I". Springer.
Brent, R.P. (1973). "Algorithms for Minimization without Derivatives". Prentice-Hall.

Development

This project was co-developed with Claude, an AI assistant by Anthropic.

License

Licensed under the MIT License.

Adding rkf78 library as a dependency

Adding `rkf78` library as a dependency