After index reduction the Fekete index-3 DAE goes Unstable at t ~ 4 of 1000, even with exact initialization; hand-written index-2 form of the same model reaches t = 1000

[ChrisRackauckas]

What happens

The Fekete-points index-3 DAE (20 particles on the unit sphere, pairwise Coulomb repulsion +
damping, 20 position-level constraints |p_i|² = 1) is handed to structural_simplify, which
performs the index reduction and returns a square 140-equation system. Every route that gets
past t = 0 then loses stability at t ≈ 4 of a [0, 1000] span — 0.4% of the requested
integration. The solver, the initialization algorithm and the tolerance do not matter; the
failure lands in t ∈ [4.0, 4.6] every time.

Meanwhile the same model, hand-written as an index-2 mass-matrix ODE and started from the same
initial data, integrates to t = 1000 with Rodas5P and lands on the correct answer
(λ → −4.75, max ‖p_i‖² − 1 = 2.2e-16, matching the Fortran RADAU5 reference of the IVP test
set). So the drift that index reduction is supposed to control is not being controlled.

This is independent of #5011 (the over-prescribed initialization on the same model). This
issue is about what happens after initialization is made exact.

Minimal reproducer

Same model file as #5011. Only ModelingToolkit and OrdinaryDiffEq are needed.

fekete_model.jl (identical to the one in #5011)

# Fekete points on the sphere, 20 particles, index-3 constrained mechanical system.
# Standalone reduction of SciMLBenchmarks benchmarks/DAE/fekete.jmd.
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using LinearAlgebra

const N_ART = 20
const NEQN = 8 * N_ART          # 160 in the hand-written index-2 form
const ALPHA_DAMP = 0.5
const MTK = ModelingToolkit

# ---- hand-written RHS (used only for the consistent initial data and for the
#      mass-matrix comparison run) -------------------------------------------
function fekete_rhs!(dy, y, p, t)
    nart = N_ART
    T = eltype(dy)
    @inbounds for i in 1:nart
        lam_i = y[6*nart+i]; mu_i = y[7*nart+i]
        for k in 1:3
            dy[3*(i-1)+k] = y[3*nart+3*(i-1)+k] + 2*mu_i*y[3*(i-1)+k]
        end
        for k in 1:3
            pk = y[3*(i-1)+k]; qk = y[3*nart+3*(i-1)+k]
            force_k = -ALPHA_DAMP*qk + 2*lam_i*pk
            for j in 1:nart
                if j != i
                    rn = zero(T)
                    for m in 1:3
                        rn += (y[3*(i-1)+m] - y[3*(j-1)+m])^2
                    end
                    force_k += (pk - y[3*(j-1)+k]) / rn
                end
            end
            dy[3*nart+3*(i-1)+k] = force_k
        end
        phi_i = -one(T)
        for k in 1:3; phi_i += y[3*(i-1)+k]^2; end
        dy[6*nart+i] = phi_i
        gpq_i = zero(T)
        for k in 1:3; gpq_i += 2*y[3*(i-1)+k]*y[3*nart+3*(i-1)+k]; end
        dy[7*nart+i] = gpq_i
    end
    nothing
end

function fekete_init()
    y = zeros(NEQN)
    for i in 1:3
        a = 2pi*i/3 + pi/13; b = 3pi/8
        y[3*(i-1)+1] = cos(a)*cos(b); y[3*(i-1)+2] = sin(a)*cos(b); y[3*(i-1)+3] = sin(b)
    end
    for i in 4:10
        a = 2pi*(i-3)/7 + pi/29; b = pi/8
        y[3*(i-1)+1] = cos(a)*cos(b); y[3*(i-1)+2] = sin(a)*cos(b); y[3*(i-1)+3] = sin(b)
    end
    for i in 11:16
        a = 2pi*(i-10)/6 + pi/7; b = -2pi/15
        y[3*(i-1)+1] = cos(a)*cos(b); y[3*(i-1)+2] = sin(a)*cos(b); y[3*(i-1)+3] = sin(b)
    end
    for i in 17:20
        a = 2pi*(i-17)/4 + pi/17; b = -3pi/10
        y[3*(i-1)+1] = cos(a)*cos(b); y[3*(i-1)+2] = sin(a)*cos(b); y[3*(i-1)+3] = sin(b)
    end
    yprime = similar(y)
    fekete_rhs!(yprime, y, nothing, 0.0)
    for i in 1:N_ART
        s = 0.0
        for j in 1:3; s += y[3*(i-1)+j]*yprime[3*N_ART+3*(i-1)+j]; end
        y[6*N_ART+i] = -s/2
    end
    y
end

const y0 = fekete_init()
const tspan = (0.0, 1000.0)

# ---- MTK index-3 model, every unknown given a default ----------------------
function build_sys_defaults()
    ps = Vector{Num}(undef, 3*N_ART); qs = Vector{Num}(undef, 3*N_ART)
    ls = Vector{Num}(undef, N_ART)
    for i in 1:N_ART
        for k in 1:3
            idx = 3*(i-1)+k
            ps[idx] = only(@variables $(Symbol("p$(i)_$(k)"))(t) = y0[idx])
            qs[idx] = only(@variables $(Symbol("q$(i)_$(k)"))(t) = 0.0)
        end
        ls[i] = only(@variables $(Symbol("lam$(i)"))(t) = 0.0)
    end
    eqs = Equation[]
    for idx in 1:3*N_ART; push!(eqs, D(ps[idx]) ~ qs[idx]); end
    for i in 1:N_ART, k in 1:3
        idx = 3*(i-1)+k
        coulomb = sum((ps[idx] - ps[3*(j-1)+k]) /
                      sum((ps[3*(i-1)+m] - ps[3*(j-1)+m])^2 for m in 1:3)
                      for j in 1:N_ART if j != i)
        push!(eqs, D(qs[idx]) ~ -ALPHA_DAMP*qs[idx] + 2*ls[i]*ps[idx] + coulomb)
    end
    for i in 1:N_ART
        push!(eqs, sum(ps[3*(i-1)+k]^2 for k in 1:3) ~ 1)
    end
    @named sys_raw = ODESystem(eqs, t)
    structural_simplify(sys_raw)
end

# ---- same model, only the positions prescribed; q, lam supplied as guesses --
function build_sys_guesses()
    ps = Vector{Num}(undef, 3*N_ART); qs = Vector{Num}(undef, 3*N_ART)
    ls = Vector{Num}(undef, N_ART)
    for i in 1:N_ART
        for k in 1:3
            idx = 3*(i-1)+k
            ps[idx] = only(@variables $(Symbol("P$(i)_$(k)"))(t) = y0[idx])
            qs[idx] = only(@variables $(Symbol("Q$(i)_$(k)"))(t))     # no default
        end
        ls[i] = only(@variables $(Symbol("LAM$(i)"))(t))              # no default
    end
    eqs = Equation[]
    for idx in 1:3*N_ART; push!(eqs, D(ps[idx]) ~ qs[idx]); end
    for i in 1:N_ART, k in 1:3
        idx = 3*(i-1)+k
        coulomb = sum((ps[idx] - ps[3*(j-1)+k]) /
                      sum((ps[3*(i-1)+m] - ps[3*(j-1)+m])^2 for m in 1:3)
                      for j in 1:N_ART if j != i)
        push!(eqs, D(qs[idx]) ~ -ALPHA_DAMP*qs[idx] + 2*ls[i]*ps[idx] + coulomb)
    end
    for i in 1:N_ART
        push!(eqs, sum(ps[3*(i-1)+k]^2 for k in 1:3) ~ 1)
    end
    guesses = Dict{Any,Float64}()
    for v in qs; guesses[v] = 0.0; end
    for v in ls; guesses[v] = 0.0; end
    @named sys_raw_g = ODESystem(eqs, t)
    structural_simplify(sys_raw_g), guesses
end

A. Exact initialization, then integrate

build_sys_guesses() declares the velocities and multipliers without defaults and supplies them
as guesses, so only the 60 positions are prescribed. Initialization then solves exactly
(retcode = Success, ‖residual‖∞ = 8.5e-15) — see #5011 for why this matters. The
integration still fails:

sys_g, guesses = build_sys_guesses()
prob_g = ODEProblem(sys_g, [], tspan; guesses = guesses)

ig    = prob_g.f.initialization_data.initializeprob
isolg = solve(ig)                       # Success, residual 8.5e-15

for a in (FBDF(), Rodas5P())
    solve(prob_g, a; abstol = 1e-8, reltol = 1e-8, save_everystep = false, maxiters = Int(1e6))
end

unknowns prescribed as initial conditions : 60 / 140
initialization system                     : 80 equations, 100 unknowns
init solve retcode                        : Success
|init residual at solution|_inf           : 8.493206138382448e-15
FBDF     / positions-only ICs  retcode = Unstable     reached t = 4.0416 of 1000.0  (16.0 s)
Rodas5P  / positions-only ICs  retcode = Unstable     reached t = 4.3193 of 1000.0  (81.7 s)

with, from the integrator:

At t=4.041566968012109, dt was forced below floating point epsilon 8.881784197001252e-16.
Aborting. There is either an error in your model specification or the true solution is unstable
(or it cannot be represented in Float64 precision).

B. Original model, initialization forced

Same failure window when the over-prescribed model of #5011 is pushed past t = 0 with an
explicit initializealg:

sys_mtk = build_sys_defaults()
mtkprob = ODEProblem(sys_mtk, [], tspan)
for a in (Rodas5P(), FBDF()), ia in (BrownFullBasicInit(), ShampineCollocationInit())
    solve(mtkprob, a; abstol = 1e-8, reltol = 1e-8, save_everystep = false,
          maxiters = Int(1e6), initializealg = ia)
end

Rodas5P  / BrownFullBasicInit       retcode = Unstable     reached t = 4.0416 of 1000.0  (25.6 s)
Rodas5P  / ShampineCollocationInit  retcode = Unstable     reached t = 4.4175 of 1000.0  (17.3 s)
FBDF     / BrownFullBasicInit       retcode = Unstable     reached t = 4.5768 of 1000.0  (2.1 s)
FBDF     / ShampineCollocationInit  retcode = Unstable     reached t = 4.4175 of 1000.0  (4.1 s)

Six routes past t = 0 (two here × two solvers, plus the two in A), four solver/init
combinations, all failing in t ∈ [4.04, 4.58].

C. The same model, hand-written as an index-2 mass-matrix ODE

fekete_rhs! in the file above is the same system, written by hand at index 2 (positions,
velocities, λ and the extra multiplier μ; the constraint appears at position level and at
velocity level). Same initial data, same tolerances:

M = zeros(NEQN, NEQN); for i in 1:6*N_ART; M[i, i] = 1.0; end
mmprob = ODEProblem(ODEFunction(fekete_rhs!, mass_matrix = M), y0, tspan)
solve(mmprob, Rodas5P(); abstol = 1e-8, reltol = 1e-8, save_everystep = false, maxiters = Int(1e7))

|f(y0)|_inf over ALGEBRAIC rows (hand-written form) : 2.498001805406602e-16
Rodas5P  hand-written mass matrix  retcode = Success      reached t = 1000.0 of 1000.0  (13.2 s)  mean lambda = -4.75  max ||p_i|^2-1| = 2.220446049250313e-16

mean λ = −4.75 and max ‖p_i‖² − 1 = 2.2e-16 at t = 1000 are the correct answer — they match
the Fortran RADAU5 reference solution shipped with the IVP test set. So the model and the initial
data are fine; it is the index-reduced form that cannot be integrated.

(Caveat on this block: FBDF on the hand-written form as written above fails immediately at
t = 0, because this stripped-down script does not supply the analytical Jacobian that the
benchmark's version of mmprob carries. That is a property of this script, not evidence about
index reduction, so I am only quoting the Rodas5P line as the contrast.)

Expected

After structural_simplify performs the index reduction and initialization is solved to
8.5e-15, the integration should reach t = 1000 like the hand-written index-2 form of the
same model does. Reaching t ≈ 4 and then forcing dt below eps suggests the constraint
residual is growing without being controlled — but I have not instrumented the drift, so that
is a hypothesis, not a claim.

Version info

  [2b5f629d] DiffEqBase v7.18.1
  [961ee093] ModelingToolkit v11.39.0
  [7771a370] ModelingToolkitBase v1.67.0
  [8913a72c] NonlinearSolve v4.27.0
  [1dea7af3] OrdinaryDiffEq v7.6.0
  [0bca4576] SciMLBase v3.49.1
  [d1185830] SymbolicUtils v4.45.0
  [0c5d862f] Symbolics v7.36.0

Julia Version 1.11.8
Commit cf1da5e20e3 (2025-11-06 17:49 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 12 × Apple M2 Max
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m2)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)

(Wall times above are from a laptop that was carrying other load; the retcodes, system sizes,
residuals and t values are what matter and are reproducible.)

Where it was hit

benchmarks/DAE/fekete.jmd in SciMLBenchmarks.jl carried three formulations of this problem in
every work-precision diagram: hand-written index-2 mass-matrix ODE, hand-written DAE residual
form, and this MTK index-reduced form. The MTK formulation has been dropped from the sweeps in
SciML/SciMLBenchmarks.jl#1670 rather than published as flat, meaningless
curves, and that PR carries these diagnostics as a documented chunk. Restoring the third
formulation to the benchmark is the reason for this report.

Companion issue: #5011 — on the same reproducer, defaults on index-reduced-away states
become hard initial conditions and produce an overdetermined and inconsistent initialization
system. The two failures are independent: this one persists when initialization is exact.