Documentation | Examples | GitHub | Crates.io | PyPI

---

A library of numerical methods for solving initial value problems (IVPs)
for Rust and Python.

---

This library provides a pure Rust implementation of SciPy's solve_ivp functionality with a typed builder API for first-order, second-order, and Hamiltonian systems. It is also available as a Python package with a SciPy-compatible solve_ivp interface.

Features

Currently implemented solvers:

• DOP853: An 8th order Dormand-Prince method with step-size control and dense output.
• DOPRI5: A 5th order Dormand-Prince method with step-size control and dense output.
• RK4: The classic 4th order Runge-Kutta method with fixed step-size and cubic Hermite interpolation for dense output.
• RK23: A 3rd order Runge-Kutta method with 2nd order error estimate for step-size control.
• LSODA: An automatic Adams/BDF switching multistep method for problems that may change stiffness.
• Radau: A 5th order implicit Runge-Kutta method of Radau IIA type with step-size control and dense output.
• BDF: A variable-order (1 to 5) Backward Differentiation Formula method for stiff ODEs with adaptive step-size control and dense output.
• Symplectic methods: Fixed-step structured solvers for separable Hamiltonian and second-order systems, including Symplectic Euler, Velocity Verlet, Ruth 3, and Yoshida 4.

Installation

Rust

cargo add ivp

Python

pip install ivp-rs

Example Usage

from ivp import solve_ivp

def exponential_decay(t, y):
    return -0.5 * y

# Solve the ODE
sol = solve_ivp(exponential_decay, (0, 10), [1.0], method='RK45', rtol=1e-6, atol=1e-9)

print(f"Final time: {sol.t[-1]}")
print(f"Final state: {sol.y[:, -1]}")

use ivp::prelude::*;

struct ExponentialDecay;

impl FirstOrderSystem for ExponentialDecay {
    fn derivative(&self, _t: f64, y: &[f64], dydt: &mut [f64]) {
        dydt[0] = -0.5 * y[0];
    }
}

fn main() {
    let decay = ExponentialDecay;
    let y0 = [1.0];

    let sol = Ivp::first_order(&decay, 0.0, 10.0, &y0)
        .method(Method::DOPRI5)
        .rtol(1e-6)
        .atol(1e-9)
        .solve()
        .unwrap();
    
    println!("Final time: {}", sol.t.last().unwrap());
    println!("Final state: {:?}", sol.y.last().unwrap());
}

More complete examples:

• Rust: docs/rust-examples.md
• Python: docs/python-examples.md

Symplectic Usage

For structured mechanics problems, use the dedicated second-order or Hamiltonian builders instead of the Rust first-order Ivp::first_order(...) path.

use ivp::prelude::*;

struct HarmonicOscillator;

impl SecondOrderSystem for HarmonicOscillator {
    fn acceleration(&self, _t: f64, q: &[f64], a: &mut [f64]) {
        a[0] = -q[0];
    }
}

fn main() {
    let q0 = [1.0];
    let v0 = [0.0];

    let sol = Ivp::second_order(&HarmonicOscillator, 0.0, 20.0, &q0, &v0)
        .method(SymplecticMethod::VelocityVerlet)
        .step_size(0.05)
        .solve()
        .unwrap();

    println!("Final state: {:?}", sol.y.last().unwrap());
}

Python can also route these methods through solve_ivp. Second-order problems can use a plain acceleration callback, and Hamiltonian problems can use a callback pair:

import numpy as np
from ivp import solve_ivp

def acceleration(t, q):
    return -np.asarray(q, dtype=float)

sol = solve_ivp(
    acceleration,
    (0.0, 20.0),
    [1.0, 0.0, -0.5, 0.0, 1.0, 0.25],  # [q..., v...]
    method="VelocityVerlet",
    step_size=0.05,
    dense_output=True,
)

print(sol.sol(0.5))

def position_derivative(t, p):
    return np.array([p[0]])

def momentum_derivative(t, q):
    return np.array([-q[0]])

sol = solve_ivp(
    (position_derivative, momentum_derivative),
    (0.0, 20.0),
    [1.0, 0.0],
    method="Yoshida4",
    step_size=0.05,
)

On the Python side, invalid callback shapes now raise normal Python exceptions with explicit messages. For example, a symplectic callback returning the wrong number of values raises ValueError, and an incomplete Hamiltonian callback pair raises TypeError.

Current limitation: event handling is not yet supported for symplectic methods.