BDF solvers return Unstable at valid zero-force equilibrium (du=0)

[1-Bart-1]

Describe the bug 🐞

BDF-family solvers (FBDF, QNDF, ABDF2) return Unstable when stepping from a valid zero-force equilibrium where all derivatives are zero (du = [0, 0, ...]). The RHS evaluates correctly with no NaN, but the first adaptive step fails. Explicit solvers (e.g. Tsit5) handle the same system correctly.

A tiny perturbation to the parameters (changing a wind vector from [0,0,0] to [1e-6,0,0]) makes all BDF solvers succeed, even though the resulting force is negligible (~1e-12 N).

Expected behavior

A system at valid equilibrium (du = 0) should be steppable by all solvers. The solution should remain at the equilibrium state.

Minimal Reproducible Example 👇

Mass on a spring-damper with aerodynamic drag. Initial conditions are at exact equilibrium: rest length equals actual distance, zero velocity, zero wind.

# MWE: BDF solvers return Unstable at zero-force equilibrium
#
# BDF-family solvers (FBDF, QNDF, ABDF2) return Unstable
# when stepping from a valid zero-force equilibrium (du = 0).
# Explicit solvers (Tsit5) handle it correctly.
#
# System: mass on spring-damper with aerodynamic drag.
# At init: zero wind, zero preload → all forces zero → du = 0.
# A tiny perturbation (wind = 1e-6 m/s) makes BDF succeed.
#
# Results:
#   Solver  wind=zero   wind=1e-6
#   FBDF    Unstable    Success
#   QNDF    Unstable    Success
#   ABDF2   Unstable    Success
#   Tsit5   Success     Success

using ModelingToolkit, OrdinaryDiffEqBDF, OrdinaryDiffEqCore,
      OrdinaryDiffEqTsit5, StaticArrays, LinearAlgebra
using ModelingToolkit: t_nounits as t, D_nounits as D

# ── Wind parameter via registered array symbolic ──────────────
const MVec3 = MVector{3, Float64}

mutable struct WindSettings
    wind_vec::MVec3
end

get_wind(s::WindSettings) = s.wind_vec
@register_array_symbolic get_wind(s::WindSettings) begin
    size = (3,)
    eltype = Float64
end

# ── Build system: mass on spring with aero drag ───────────────
@parameters p_ws::WindSettings
@variables begin
    pos(t)[1:3]
    vel(t)[1:3]
    wind_vec(t)[1:3]
    seg_vec(t)[1:3]
    seg_len(t)
    unit_vec(t)[1:3]
    spring_force(t)[1:3]
    va(t)[1:3]
    va_perp(t)[1:3]
    drag(t)[1:3]
    total_force(t)[1:3]
end

anchor = [0.0, 0.0, 5.0]
k = 50.0         # spring stiffness [N/m]
c = 5.0          # damping [Ns/m]
l0 = 3.0         # rest length [m] (= actual distance, no preload)
m = 1.0          # point mass [kg]
rho = 1.225      # air density [kg/m³]
cd_tether = 0.958
diam = 0.005     # tether diameter [m]

eqs = [
    wind_vec ~ get_wind(p_ws)
    seg_vec ~ anchor - pos
    seg_len ~ norm(collect(seg_vec))
    unit_vec ~ seg_vec / seg_len
    spring_force ~ (k * (seg_len - l0) + c * (vel ⋅ unit_vec)) * unit_vec
    va ~ wind_vec - vel
    va_perp ~ va - (va ⋅ unit_vec) * unit_vec
    drag ~ (0.5 * rho * cd_tether * norm(collect(va)) *
            seg_len * diam) * va_perp
    total_force ~ spring_force + drag
    D(pos) ~ vel
    D(vel) ~ total_force / m
]

@named sys = System(eqs, t)
ssys = structural_simplify(sys)

# ── Build problem ────────────────────────────────────────────
ws = WindSettings(MVec3(0.0, 0.0, 0.0))
defs = Dict(pos[1] => 0.0, pos[2] => 0.0, pos[3] => 2.0,
            vel[1] => 0.0, vel[2] => 0.0, vel[3] => 0.0,
            p_ws => ws)
prob = ODEProblem(ssys, defs, (0.0, 0.05))

# ── Solver matrix ────────────────────────────────────────────
solvers = [
    "FBDF"            => FBDF(),
    "QNDF"            => QNDF(),
    "ABDF2"           => ABDF2(),
    "Tsit5"           => Tsit5(),
]

wind_cases = [
    "zero" => MVec3(0.0, 0.0, 0.0),
    "1e-6" => MVec3(1e-6, 0.0, 0.0),
]

set_p = ModelingToolkit.setp(ssys, ssys.p_ws)

# Print header
print(rpad("Solver", 8))
for (w, _) in wind_cases
    print(rpad("wind=$w", 16))
end
println()
println("-"^72)

for (solver_name, solver) in solvers
    print(rpad(solver_name, 8))
    for (_, wind_val) in wind_cases
        local integ
        ws.wind_vec = wind_val
        integ = OrdinaryDiffEqCore.init(prob, solver;
            abstol=0.01, reltol=0.01,
            save_on=false, save_everystep=false)
        set_p(integ, ws)
        OrdinaryDiffEqCore.reinit!(integ; reinit_dae=true)
        OrdinaryDiffEqCore.step!(integ, 0.01, true)
        print(rpad(string(integ.sol.retcode), 16))
    end
    println()
end

Error & Stacktrace ⚠️

No exception — the solver returns Unstable retcode silently:

Solver  wind=zero       wind=1e-6       
------------------------------------------------------------------------
FBDF    Unstable        Success         
QNDF    Unstable        Success         
ABDF2   Unstable        Success         
Tsit5   Success         Success         

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

(mwes) pkg> st
Status `~/Code/Kite/SymbolicAWEModels.jl/mwes/Project.toml`
  [13f3f980] CairoMakie v0.15.9
  [5ae59095] Colors v0.13.1
  [aaaaaaaa] ControlSystemsBase v1.20.3
  [e9467ef8] GLMakie v0.13.9
  [a98d9a8b] Interpolations v0.16.2
  [b964fa9f] LaTeXStrings v1.4.0
  [ee78f7c6] Makie v0.24.9
  [961ee093] ModelingToolkit v11.20.0
  [2774e3e8] NLsolve v4.5.1
  [6ad6398a] OrdinaryDiffEqBDF v1.24.0
  [bbf590c4] OrdinaryDiffEqCore v3.29.0
  [b1df2697] OrdinaryDiffEqTsit5 v1.11.0
  [90137ffa] StaticArrays v1.9.18
  [10745b16] Statistics v1.11.1
  [0c5d862f] Symbolics v7.18.1

• Output of versioninfo()

julia> versioninfo()
Julia Version 1.12.5
Commit 5fe89b8ddc1 (2026-02-09 16:05 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × 11th Gen Intel(R) Core(TM) i7-11850H @ 2.50GHz
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, tigerlake)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 16 virtual cores)

[ChrisRackauckas-Claude]

Root cause: norm(va) produces a NaN Jacobian at va = 0

I reproduced this and traced it to a non-differentiable point in the RHS, not a bug in the BDF solvers.

What's happening

The drag equation contains norm(collect(va)) where va = wind_vec - vel. At the initial state with zero wind and zero velocity, va == [0,0,0], and the 2-norm is non-differentiable at the origin.

BDF solvers build a Newton system W = M/(γ·dt) - J, where J = ∂f/∂u is computed by ForwardDiff by default. The ForwardDiff rule for norm(x) is x / norm(x), which evaluates to 0/0 = NaN at x = 0. That NaN propagates into J, then W, then the Newton increment dz = W \ b, then ndz = ‖dz‖ — which is NaN. In nlsolve!, !isfinite(ndz) marks the solve as Divergence, the step is rejected, dt is halved, and this repeats until dt reaches subnormal floats. check_error then returns ReturnCode.Unstable.

Tsit5 works because it never builds a Jacobian. A tiny wind perturbation works because it moves va off the origin so the derivative is well-defined. Confirmation via ForwardDiff on the actual RHS:

Jacobian at u0 (wind = zero):
6×6 Matrix{Float64}:
 NaN  NaN  NaN  NaN  NaN  NaN
 NaN  NaN  NaN  NaN  NaN  NaN
 NaN  NaN  NaN  NaN  NaN  NaN
 1.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0

With wind = 1e-6, the same Jacobian is fully finite.

Workarounds

Two things fix it for your MWE without changing the physics meaningfully:

1. Smooth the norm (preferred — fixes the model):

smooth_norm(v) = sqrt(sum(abs2, v) + 1e-12^2)
...
drag ~ (0.5 * rho * cd_tether * smooth_norm(collect(va)) * seg_len * diam) * va_perp

After this change, all four solvers return Success on wind=zero.

2. Switch Jacobian to finite differences (keeps the non-smooth model):

using ADTypes
FBDF(autodiff = AutoFiniteDiff())

FiniteDiff perturbs u by ε and never evaluates the symbolic va/‖va‖ form, so it doesn't hit the 0/0. I confirmed FBDF, QNDF, and ABDF2 all succeed with AutoFiniteDiff() on the original model.

Should OrdinaryDiffEq do anything here?

This is a user-model issue (aerodynamic drag with ‖v_a‖ is a classic non-smooth point at the origin), but the failure mode is unusually opaque: ReturnCode.Unstable with no indication that the Jacobian was the problem. A cheap improvement would be detecting !isfinite entries in W/J on the very first Newton step and surfacing a diagnostic like "Jacobian evaluation produced NaN/Inf — often caused by non-smooth terms (e.g. norm(x) at x=0) at the initial state. Consider smoothing the RHS or using autodiff=AutoFiniteDiff()." Happy to send a PR for that if it sounds useful.