Summary
ESDIRK659L2SA fails to integrate any problem in adaptive mode — it runs to MaxIters (~10⁶ steps) and returns a wrong answer. The method works correctly with fixed dt, so the RK core is fine; the bug is in the embedded error estimator coefficients (btilde).
Reproduction
using OrdinaryDiffEqSDIRK
prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
# adaptive: fails
sol = solve(prob, ESDIRK659L2SA(), reltol = 1e-6, abstol = 1e-6)
@show sol.retcode # MaxIters
@show sol.stats.naccept # 999979
@show sol.u[end] # 0.575... (wrong; exp(-1) = 0.3679)
# fixed dt: correct
sol2 = solve(prob, ESDIRK659L2SA(), dt = 0.1, adaptive = false)
@show sol2.u[end] # 0.36787944118... ✓
Root cause
A valid embedded pair has btilde = b - b̂ with sum(b) = sum(b̂) = 1, so sum(btilde) = 0. ESDIRK659L2SA violates this:
| Method |
sum(btilde) |
| ESDIRK54I8L2SA |
5e-16 ✓ |
| ESDIRK436L2SA2 |
-1e-17 ✓ |
| ESDIRK437L2SA |
4e-17 ✓ |
| ESDIRK547L2SA2 |
0.0 ✓ |
| ESDIRK659L2SA |
-1.69 ✗ |
sum(bi) = 1.0 is correct (hence fixed-dt works), but the btilde coefficients in ESDIRK659L2SATableau (lib/OrdinaryDiffEqSDIRK/src/sdirk_tableaus.jl) imply sum(b̂) ≈ 2.69. The error estimate is therefore meaningless, the step controller can never satisfy the tolerance, and dt collapses until MaxIters.
Notes
- Pre-existing: the bad coefficients date to the original "Added SDIRK solvers" commit (
e32dd5a842), long before the recent tableau-form migration. The migration copied them verbatim.
- Fix requires the correct embedded
b̂ row from the reference (Kennedy & Carpenter ESDIRK6(5)9). Once corrected, sum(btilde) should be ~0 and adaptive solves should work.
- Suggest adding a tableau sanity test asserting
sum(btilde) ≈ 0 (and sum(bi) ≈ 1) for every adaptive SDIRK/ESDIRK method to catch this class of error.
cc @ChrisRackauckas
Summary
ESDIRK659L2SAfails to integrate any problem in adaptive mode — it runs toMaxIters(~10⁶ steps) and returns a wrong answer. The method works correctly with fixeddt, so the RK core is fine; the bug is in the embedded error estimator coefficients (btilde).Reproduction
Root cause
A valid embedded pair has
btilde = b - b̂withsum(b) = sum(b̂) = 1, sosum(btilde) = 0. ESDIRK659L2SA violates this:sum(bi) = 1.0is correct (hence fixed-dt works), but thebtildecoefficients inESDIRK659L2SATableau(lib/OrdinaryDiffEqSDIRK/src/sdirk_tableaus.jl) implysum(b̂) ≈ 2.69. The error estimate is therefore meaningless, the step controller can never satisfy the tolerance, anddtcollapses untilMaxIters.Notes
e32dd5a842), long before the recent tableau-form migration. The migration copied them verbatim.b̂row from the reference (Kennedy & Carpenter ESDIRK6(5)9). Once corrected,sum(btilde)should be ~0 and adaptive solves should work.sum(btilde) ≈ 0(andsum(bi) ≈ 1) for every adaptive SDIRK/ESDIRK method to catch this class of error.cc @ChrisRackauckas