Benchmarking Stiff DAEs: 8 Problems I Added to SciMLBenchmarks.jl

Most of my GSoC work splits into two halves: making solvers faster (in OrdinaryDiffEq.jl) and proving which solvers are faster (in SciMLBenchmarks.jl). This post is about the second half — the eight differential-algebraic equation (DAE) benchmarks I added.

Why DAEs are hard

An ODE is u' = f(u, t). A DAE mixes differential equations with algebraic constraints:

M u' = f(u, t)        # M is singular → some equations are pure constraints

The index measures how many times you have to differentiate the constraints to recover an ODE. Index-1 is gentle; index-2 and index-3 problems (mechanical systems with position/velocity/acceleration constraints) are where solvers earn their keep — and where they quietly fall over.

The problems I contributed

I worked through classic problems from the IVP Test Set, the standard suite for stiff/DAE solvers:

• Andrews' squeezing mechanism — a 7-body planar linkage, index-3.
• Car Axis — index-3 multibody.
• Wheelset — a railway wheel-rail contact model.
• Water Tube System — a 49-dimensional index-2 hydraulic network.
• Slider-crank (Simeon 1998) — index-2.
• Charge Pump — an index-2 electrical circuit.
• Two-Bit Adding Unit — a 350-variable digital-circuit DAE.
• Fekete problem — constrained optimization on a sphere.

Three formulations each

The interesting part: each problem went in with multiple formulations so we can compare solver families fairly:

1. Mass-matrix form M u' = f — for solvers that accept a singular M.
2. Residual (fully implicit) form g(u', u, t) = 0 — for IDA / DAE-native solvers.
3. ModelingToolkit (MTK) — symbolic, so the framework can simplify and generate Jacobians automatically.

# Mass-matrix form
f = ODEFunction(rhs!; mass_matrix = M)
prob = ODEProblem(f, u0, tspan, p)

# ...vs the same physics as a fully-implicit residual
prob_dae = DAEProblem(residual!, du0, u0, tspan, p)

Getting all three to agree on the same trajectory (to a tight tolerance) is its own debugging adventure — a sign error in the residual form shows up as a solver that "converges" to the wrong answer.

What benchmarks actually teach you

A benchmark isn't just a number. Building these taught me:

• Consistent initial conditions matter. For index-2/3 DAEs you can't pick u0 freely — it has to satisfy the hidden constraints, or the solver starts in an invalid state.
• Work-precision diagrams (error vs wall-clock) are the only honest way to compare solvers — a method that's fast at low accuracy can be hopeless at high accuracy.
• Reproducibility is a feature. Every benchmark has to run start-to-finish in CI on a fresh machine.

These eight problems now run in the SciML suite, helping decide which solver gets recommended for stiff, constrained systems. That's a small but real contribution to every scientist who calls solve.