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

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

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

References

1. Fehlberg, E. (1968). "Classical Fifth-, Sixth-, Seventh-, and Eighth-Order Runge-Kutta Formulas with Stepsize Control". NASA TR R-287.
2. Hairer, E., Norsett, S.P., & Wanner, G. (1993). "Solving Ordinary Differential Equations I". Springer.
3. 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.