Add backward-difference stiff interpolation for BDF solvers

[ChrisRackauckas-Claude]

Summary

• Replace the default Hermite (cubic) interpolation fallback with native Newton backward difference interpolation for FBDF, QNDF, and DFBDF solvers
• Uses the formula p(Θ) = y₁ + Σ φⱼ(Θ-1) * k[j] where backward differences come from the solver's own history points, giving interpolation accuracy consistent with the solver order
• Supports both function values (Val{0}) and first derivatives (Val{1})
• Follows the same pattern established by the Rosenbrock stiff interpolation (has_stiff_interpolation trait, _ode_interpolant dispatch)

Files modified (6)

• OrdinaryDiffEqBDF.jl — Added imports and includes for the 3 new files
• alg_utils.jl — Added has_stiff_interpolation trait for QNDF, FBDF, DFBDF
• bdf_caches.jl — Added dense field to QNDFCache and FBDFCache for aliasing integrator.k
• dae_caches.jl — Added dense field to DFBDFCache
• bdf_perform_step.jl — Modified initialize! (kshortsize = max_order) and perform_step! (store backward differences in integrator.k) for QNDF and FBDF
• dae_perform_step.jl — Same for DFBDF

Files created (3)

• interp_func.jl — SciMLBase.interp_summary overrides for all BDF cache types
• bdf_interpolants.jl — Newton backward difference interpolation (Val{0} and Val{1})
• stiff_addsteps.jl — No-op _ode_addsteps! (backward diffs precomputed in perform_step!)

Test plan

• BDF Regression Tests (8/8 pass + autodiff)
• BDF Convergence Tests (72/72 pass)
• DAE Convergence Tests (20/20 pass)
• DAE AD Tests (16/16 pass)
• DAE Event Tests (pass)
• DAE Initialization Tests (pass)
• Inference Tests (pass)
• Dense interpolation accuracy ~1e-10 for QNDF and FBDF
• First derivative accuracy ~1e-10
• saveat with interpolation works correctly
• Runic formatting applied

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Update: Lagrange interpolation + Dense tests

Changes in this commit:

FBDF/DFBDF interpolation switched to direct Lagrange formula:

• bdf_interpolants.jl: Split into QNDF (Newton backward diff) and FBDF/DFBDF (Lagrange) paths
• Direct Lagrange formula: L_j(Θ) = Π_{m≠j} (σ+m)/(m-j) where denominators are constant integers
• This avoids barycentric division singularities that cause NaN with ForwardDiff dual numbers
• Same polynomial as calc_Lagrange_interp in bdf_utils.jl

DFBDF calck blocks updated:

• dae_perform_step.jl: Store past solution values (uprev, u_corrector columns) in k instead of backward differences, matching FBDF pattern

Dense output tests added:

• ode_dense_tests.jl: Added QNDF and FBDF regression tests with empirically verified tolerances
  • QNDF: 1D ~0.009, 2D ~0.01, deriv ~1e-7
  • FBDF: 1D ~0.046, 2D ~0.076, deriv ~1e-7

BDF not added to special_interps.jl:

• That test compares adaptive (MutableCache) vs non-adaptive (ConstantCache) interpolation consistency
• For multi-step methods, interpolation data (D matrix, u_history) depends on cumulative integration history
• With adaptive=false, BDF order stays at 1 (no adaptation), yielding fundamentally different interpolation
• This is inherent to multi-step methods — all other methods in that test are single-step

All existing OrdinaryDiffEqBDF tests pass.