Skip to content
24 changes: 24 additions & 0 deletions .github/workflows/draft-pdf.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
name: Draft PDF
on: [push]

jobs:
paper:
runs-on: ubuntu-latest
name: Paper Draft
steps:
- name: Checkout
uses: actions/checkout@v4
- name: Build draft PDF
uses: openjournals/openjournals-draft-action@master
with:
journal: joss
# This should be the path to the paper within your repo.
paper-path: paper/paper.md
- name: Upload
uses: actions/upload-artifact@v4
with:
name: paper
# This is the output path where Pandoc will write the compiled
# PDF. Note, this should be the same directory as the input
# paper.md
path: paper/paper.pdf
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ The following ODE solvers are available in diffsol

1. A variable order Backwards Difference Formulae (BDF) solver, suitable for stiff problems and singular mass matrices. The basic algorithm is derived in [(Byrne & Hindmarsh, 1975)](#1), however this particular implementation follows that implemented in the Matlab routine ode15s [(Shampine & Reichelt, 1997)](#4) and the SciPy implementation [(Virtanen et al., 2020)](#5), which features the NDF formulas for improved stability
2. A Singly Diagonally Implicit Runge-Kutta (SDIRK or ESDIRK) solver, suitable for moderately stiff problems and singular mass matrices. Two different butcher tableau are provided, TR-BDF2 [(Hosea & Shampine, 1996)](#2) and ESDIRK34 [(Jørgensen et al., 2018)](#3), or users can supply their own.
3. A variable order Explict Runge-Kutta (ERK) solver, suitable for non-stiff problems. One butcher tableau is provided, the 4th order TSIT45 [(Tsitouras, 2011)](#5), or users can supply their own.
3. A variable order Explicit Runge-Kutta (ERK) solver, suitable for non-stiff problems. One butcher tableau is provided, the 4th order TSIT45 [(Tsitouras, 2011)](#5), or users can supply their own.

All solvers feature:

Expand Down
2 changes: 2 additions & 0 deletions paper/.gitignore
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
*.pdf
jats
180 changes: 180 additions & 0 deletions paper/paper.bib
Original file line number Diff line number Diff line change
@@ -0,0 +1,180 @@
@software{cranelift,
author = {{Bytecode Alliance}},
title = {cranelift},
url = {https://wasmtime.dev/},
version = {0.115.1},
date = {2025-03-13},
}

@software{nalgebra,
author = {{Dimforge}},
title = {nalgebra},
url = {https://github.com/dimforge/nalgebra},
version = {0.33.2},
date = {2025-03-13},
}

@software{faer,
author = {Sarah El Kazdadi},
title = {faer-rs},
url = {https://faer-rs.github.io/},
version = {0.1.0},
date = {2025-03-13},
}

@article{virtanen2020scipy,
title={SciPy 1.0: fundamental algorithms for scientific computing in Python},
author={Virtanen, Pauli and Gommers, Ralf and Oliphant, Travis E and Haberland, Matt and Reddy, Tyler and Cournapeau, David and Burovski, Evgeni and Peterson, Pearu and Weckesser, Warren and Bright, Jonathan and others},
journal={Nature methods},
volume={17},
number={3},
pages={261--272},
year={2020},
publisher={Nature Publishing Group US New York},
doi={10.1038/s41592-019-0686-2}
}

@article{shampine1997matlab,
title={The matlab ode suite},
author={Shampine, Lawrence F and Reichelt, Mark W},
journal={SIAM journal on scientific computing},
volume={18},
number={1},
pages={1--22},
year={1997},
publisher={SIAM},
doi = {10.1137/s1064827594276424}
}

@article{gardner2022sundials,
title = {Enabling new flexibility in the {SUNDIALS} suite of nonlinear and differential/algebraic equation solvers},
author = {Gardner, David J and Reynolds, Daniel R and Woodward, Carol S and Balos, Cody J},
journal = {ACM Transactions on Mathematical Software (TOMS)},
publisher = {ACM},
volume = {48},
number = {3},
pages = {1--24},
year = {2022},
doi = {10.1145/3539801}
}

@article{DifferentialEquations.jl-2017,
author = {Rackauckas, Christopher and Nie, Qing},
doi = {10.5334/jors.151},
journal = {The Journal of Open Research Software},
keywords = {Applied Mathematics},
note = {Exported from https://app.dimensions.ai on 2019/05/05},
number = {1},
pages = {},
title = {DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia},
url = {https://app.dimensions.ai/details/publication/pub.1085583166 and http://openresearchsoftware.metajnl.com/articles/10.5334/jors.151/galley/245/download/},
volume = {5},
year = {2017}
}

@phdthesis{kidger2021on,
title={{O}n {N}eural {D}ifferential {E}quations},
author={Patrick Kidger},
year={2021},
school={University of Oxford},
doi={h10.48550/arXiv.2202.02435}
}

@article{byrne1975polyalgorithm,
title={A polyalgorithm for the numerical solution of ordinary differential equations},
author={Byrne, George D. and Hindmarsh, Alan C.},
journal={ACM Transactions on Mathematical Software (TOMS)},
volume={1},
number={1},
pages={71--96},
year={1975},
publisher={ACM New York, NY, USA},
doi={10.1145/355626.355636}
}

@article{bank1985transient,
title={Transient simulation of silicon devices and circuits},
author={Bank, Randolph E and Coughran, William M and Fichtner, Wolfgang and Grosse, Eric H and Rose, Donald J and Smith, R Kent},
journal={IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems},
volume={4},
number={4},
pages={436--451},
year={1985},
publisher={IEEE},
doi={10.1109/T-ED.1985.22232}
}

@article{hosea1996analysis,
title={Analysis and implementation of TR-BDF2},
author={Hosea, ME and Shampine, LF},
journal={Applied Numerical Mathematics},
volume={20},
number={1-2},
pages={21--37},
year={1996},
publisher={Elsevier},
doi={10.1016/0168-9274(95)00115-8}
}

@article{jorgensen2018family,
title={A family of ESDIRK integration methods},
author={J{\o}rgensen, John Bagterp and Kristensen, Morten Rode and Thomsen, Per Grove},
journal={arXiv preprint arXiv:1803.01613},
year={2018},
doi= {10.48550/arXiv.1803.01613}
}

@inproceedings{moses22enzyme,
author = {Moses, William S. and Narayanan, Sri Hari Krishna and Paehler, Ludger and Churavy, Valentin and Schanen, Michel and H\"{u}ckelheim, Jan and Doerfert, Johannes and Hovland, Paul},
title = {Scalable Automatic Differentiation of Multiple Parallel Paradigms through Compiler Augmentation},
year = {2022},
isbn = {9784665454445},
publisher = {IEEE Press},
booktitle = {Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis},
articleno = {60},
numpages = {18},
keywords = {automatic differentiation, tasks, OpenMP, compiler, Julia, parallel, Enzyme, C++, RAJA, hybrid parallelization, MPI, distributed, LLVM},
location = {Dallas, Texas},
series = {SC '22},
doi = {10.1109/SC41404.2022.00065}
}

@inproceedings{lattner2004llvm,
title={LLVM: A compilation framework for lifelong program analysis \& transformation},
author={Lattner, Chris and Adve, Vikram},
booktitle={International symposium on code generation and optimization, 2004. CGO 2004.},
pages={75--86},
year={2004},
organization={IEEE},
doi={10.1109/CGO.2004.1281665}
}

@article{tsitouras2011runge,
title={Runge--Kutta pairs of order 5 (4) satisfying only the first column simplifying assumption},
author={Tsitouras, Ch},
journal={Computers \& mathematics with applications},
volume={62},
number={2},
pages={770--775},
year={2011},
publisher={Elsevier},
doi={10.1016/j.camwa.2011.06.002}
}

@software{PyO3_Project_and_Contributors,
author = {{PyO3 Project and Contributors}},
license = {["Apache-2.0", "MIT"]},
url = {https://github.com/PyO3},
title = {{PyO3}}
}

@article{fritzson2020OpenModelica,
author = {Fritzson, Peter and Pop, Adrian and Abdelhak, Karim and Ashgar, Adeel and Bachmann, Bernhard and Braun, Willi and Bouskela, Daniel and Braun, Robert and Buffoni, Lena and Casella, Francesco and Castro, Rodrigo and Franke, Rüdiger and Fritzson, Dag and Gebremedhin, Mahder and Heuermann, Andreas and Lie, Bernt and Mengist, Alachew and Mikelsons, Lars and Moudgalya, Kannan and Ochel, Lennart and Palanisamy, Arunkumar and Ruge, Vitalij and Schamai, Wladimir and Sjölund, Martin and Thiele, Bernhard and Tinnerholm, John and Östlund, Per},
doi = {10.4173/mic.2020.4.1},
journal = {Modeling, Identification and Control},
number = {4},
pages = {241--295},
title = {{The OpenModelica Integrated Environment for Modeling, Simulation, and Model-Based Development}},
volume = {41},
year = {2020}
}
79 changes: 79 additions & 0 deletions paper/paper.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,79 @@
---
title: '`diffsol`: Rust crate for solving differential equations'
tags:
- rust
- scientific computing
- solver
- ordinary differential equations
- differential algebraic equations
- runge-kutta
- backward differentiation formula
authors:
- name: Martin Robinson
orcid: 0000-0002-1572-6782
corresponding: true
affiliation: 1
- name: Alex Allmont
orcid: 0009-0001-4162-0180
affiliation: 1


affiliations:
- name: Oxford Research Software Engineering Group, Doctoral Training Centre, University of Oxford, Oxford, UK
index: 1
ror: 052gg0110
date: 6 Nov 2025
bibliography: paper.bib
---

# Summary

Ordinary Differential Equations (ODEs) are powerful tools for modelling a wide range of physical systems. Unlike purely data-driven models, ODEs can be based on the underlying physics, biology, or chemistry of the system being modelled, making them particularly useful for predicting the behaviour of a system under conditions that have not been observed. ODEs can be used to model everything from the motion of planets to the spread of infectious diseases.

[`diffsol`](https://github.com/martinjrobins/diffsol) is a Rust crate for solving ordinary differential equations (ODEs) or semi-explicit differential algebraic equations (DAEs). It can solve equations in the following form:

$$
M \frac{dy}{dt} = f(t, y, p)
$$

where $y$ is the state of the system, $p$ are a set of parameters, $t$ is time, $f(t, y, p)$ is a function that describes how the state of the system changes over time and $M$ is an optional and possibly singular mass matrix. The solution to an ODE is a function $y(t)$ that satisfies the ODE and any initial conditions.

The equations (e.g. $f(t, y, p)$) can be provided by the user either using Rust code, or a custom Domain Specific Language (DSL) called `DiffSL`. `DiffSL` uses automatic differentiation using Enzyme [@moses22enzyme] to calculate the necessary gradients, and JIT compilation (using either `LLVM` [@lattner2004llvm] or `Cranelift` [@cranelift]) to generate efficient native code at runtime. The DSL is ideal for using `diffsol` from a higher-level language like Python or R while still maintaining similar performance to pure rust. Diffsol currently provides Python bindings through the [`pydiffsol`](https://github.com/alexallmont/pydiffsol) package, with further language bindings planned.

ODE solvers require linear algebra containers (e.g. vectors, matrices), operators and linear solvers. `diffsol` allows users to choose both dense and sparse matrices and solvers from the `nalgebra` [@nalgebra] or `faer` [@faer] crates, and uses a trait-based approach to allow other linear algebra libraries to be added at a later date.

# Statement of need

ODE solvers have a long history in scientific computing, and many libraries currently exist. Some notable examples include `scipy.integrate.odeint` [@virtanen2020scipy] in Python, `ode45` [@shampine1997matlab] in MATLAB, and the `Sundials` suite of solvers [@gardner2022sundials] in C. Rust is a systems programming language that is gaining popularity in the scientific computing community due to its performance, safety, and ease of use. There is currently no ODE solver library written in Rust that provides the same level of functionality as these other libraries, and this is the gap that `diffsol` aims to fill.

ODE solvers written in lower-level languages like C, Fortran or Rust offer significant performance benefits. However, these solvers are often more difficult to wrap and use in higher-level languages like Python or MATLAB, primarily because users must supply their equations in the language of the solver. `diffsol` solves this issue by providing its own custom `DiffSL` DSL which is JIT compiled to efficient native code at run-time, meaning that users from a higher-level language like Python or R can specify their equations using a simple string-based format while still maintaining similar performance to pure Rust. Two other popular ODE solvers that take advantage of JIT compilation are `DifferentialEquations.jl` [@DifferentialEquations.jl-2017] in Julia, and `diffrax` [@kidger2021on] in Python. However, both these packages compile the entire solver as well as the equations, which is a significant amount of code. `diffSol` only compiles the equations, meaning that it has a significantly smaller "time-to-first-plot" for users. Another popular differential equations solver package utilising a DSL is OpenModelica [@fritzson2020OpenModelica]. Wrappers to this package in higher-level languages like Python rely on messaging to a separate OpenModelica server, which can be slow and more complicated to set up. In contrast, `diffsol` can be integrated directly into higher-level languages using language bindings and linking to a single shared library, see for example the `pydiffsol` Python bindings discussed below.

# Features

The following solvers are available in `diffsol`:

1. A variable order Backwards Difference Formulae (BDF) solver, suitable for stiff problems and singular mass matrices. The basic algorithm is derived in [@byrne1975polyalgorithm], however this particular implementation follows that implemented in the Matlab routine ode15s [@shampine1997matlab] and the SciPy implementation [@virtanen2020scipy], which features the NDF formulas for improved stability
2. A Singly Diagonally Implicit Runge-Kutta (SDIRK or ESDIRK) solver, suitable for moderately stiff problems and singular mass matrices. Two different butcher tableau are provided, TR-BDF2 [@bank1985transient; @hosea1996analysis] and ESDIRK34 [@jorgensen2018family], or users can supply their own.
3. A variable order Explicit Runge-Kutta (ERK) solver, suitable for non-stiff problems. One butcher tableau is provided, the 4th order TSIT45 [@tsitouras2011runge], or users can supply their own.

All solvers feature:

- Linear algebra containers and linear solvers from the `nalgebra` or `faer` crates, including both dense and sparse matrix support.
- Adaptive step-size control to given relative and absolute tolerances. Tolerances can be set separately for the main equations, quadrature of the output function, and sensitivity analysis.
- Dense output, interpolating to times provided by the user.
- Event handling, stopping when a given condition $g_e(t, y, p)$ is met or at a specific time.
- Numerical quadrature of an optional output $g_o(t, y, p)$ function over time.
- Forward sensitivity analysis, calculating the gradient of an output function or the solver states $y$ with respect to the parameters $p$.
- Adjoint sensitivity analysis, calculating the gradient of cost function $G(p)$ with respect to the parameters $p$. The cost function can be the integral of a continuous output function $g(t, y, p)$ or a sum of a set of discrete functions $g_i(t_i, y, p)$ at time points $t_i$.

# Bindings in higher-level languages

[`pydiffsol`](https://github.com/alexallmont/pydiffsol) provides Python bindings to `diffsol` using the `PyO3` [@PyO3_Project_and_Contributors] crate. It allows users to define ODEs in Python using the `DiffSL` DSL, and solve them using the Rust `diffsol` library. `pydiffsol` aims to provide a simple and easy-to-use interface for solving ODEs in Python, while still maintaining the performance benifits of using Rust under the hood.

The goal is to develop further [bindings to other higher-level languages](https://github.com/martinjrobins/diffsol/issues/131), including R, JAX and JavaScript, exploiting the `DiffSL` DSL to provide high performance while maintaining ease of use and positioning the core `diffsol` library as a widely used cross-language, cross-platform, high-performance ODE solver.

# Acknowledgements

We gratefully acknowledge the support of all the past and future contributors to the `diffsol` project, for their advice, enthusiasm, bug reports and code. In particular, we would like to thank the authors of the `pharmsol` crate, Julian Otalvaro and Markus Hovd.

# References