Harsh Singh

Full-Stack Developer & Open Source Contributor

I build performant web apps and contribute to scientific computing — React on the front, Julia & Go under the hood.

H
All posts
·7 min read

Implementing IMEX and Multirate Integrators in OrdinaryDiffEq.jl

  • Julia
  • SciML
  • Numerical Methods
  • ODE

A big chunk of my OrdinaryDiffEq.jl work has been adding IMEX and multirate integrators. Both come from the same observation: a single system often contains terms with wildly different character, and forcing one integration strategy on all of them is wasteful.

IMEX: implicit where you must, explicit where you can

In many PDE discretizations you can split the right-hand side:

u' = f_explicit(u) + f_implicit(u)

The stiff part (say, diffusion) needs an implicit solve for stability; the non-stiff part (advection, reactions) is cheaper done explicitly. IMEX Runge–Kutta methods do exactly that in one scheme.

I implemented several families:

  • ARS schemes (ARS222, ARS232, ARS443) — Ascher–Ruuth–Spiteri.
  • BHR(5,5,3) — Boscarino–Russo.
  • The IMEX-SSP (strong-stability-preserving) family.

The unifying idea was an ESDIRKIMEXTableau abstraction so all of these share one generic perform_step! instead of each method hand-rolling its own stepping code.

# One generic stepper, driven by a tableau, instead of N copies
function perform_step!(integrator, cache::ESDIRKIMEXCache)
    # explicit stages + implicit Newton solves, read from the tableau
end

Multirate: different clocks for different terms

Multirate methods go further: when part of your system evolves much faster than the rest, you take small steps for the fast part and large steps for the slow part. I added:

  • MRI-GARK (Multirate Infinitesimal Generalized-structure ARK) — ERK22a/22b, Sandu 2019.
  • MIS2 — Multirate Infinitesimal Step (Wensch–Knoth–Galant).
  • MRAB — a multirate Adams–Bashforth method.
  • MREEF — a multirate extrapolated explicit Euler.

The fun (and the hard part) is the slow–fast coupling: the slow method produces a "forcing" function that the fast integrator sees as it sub-steps. Get the accumulation of that forcing wrong and you silently lose an order of accuracy — which only a careful convergence-order test catches.

Why this matters

These methods aren't academic curiosities — they're what make large stiff simulations (climate, combustion, circuits) tractable. Contributing them taught me to:

  • Read papers as specifications. A Butcher tableau in a 2009 paper becomes a tested method in a library.
  • Trust convergence tests over intuition. "It looks right" is not a proof; an order-of-accuracy plot is.
  • Design for the next method. The best PR isn't one new solver — it's an abstraction that makes the next ten solvers easy.

Designed & built by Harsh Singh · singhharsh.in