Fix RK45 dense-output FSAL corruption and restore local extrapolation

[jtravs]

Two independent defects in the DOPRI5 integrator, as two separate commits because their blast radii are very different: the first leaves the stepped trajectory bit-identical, the second changes the result of every simulation.

1. 7747c93 — the FSAL move ran before the dense output was consumed

step! performed the FSAL move (k7 of the completed step becomes k1 of the next) inside its accept branch, but solve consumes that step's dense output after step! returns — both the output=true saves and the interpolant closure handed to stepfun. Every interpolated save therefore evaluated the dense-output polynomial with k7 sitting in the k1 slot, and the k1 weight is not zero (interpC column 1, b1(1) = 35/384), so each value carried a spurious dt*b1(σ)*(k7 − k1) = O(h²) term. The dense output was 2nd order instead of 4th.

Not a corner case: Luna.run leaves step_on unset by default and Output.GridCondition saves via yfun(ts) with ts < tn, so every saved z-slice of a default run is an interpolated value. Runs that pass step_on aligned to the save grid were unaffected — those saves are exact step endpoints.

The fix moves it to the top of step!, guarded by the previous step's s.ok. step! is evaluate!'s only caller, so this is one site with one guard, and it stays ahead of evaluate!(::PreconStepper)'s rebase of ks[1] into the new anchor frame — the copy must still happen in the old frame, exactly as before.

2. 1d7e4c3 — locextrap=true propagated the lower-order solution

dopri.jl had the two weight vectors of the Dormand–Prince RK5(4)7M pair labelled the wrong way round. b5, commented "Weights for 5th order method", held [5179/57600, …], which is the embedded 4th-order vector; b4 held [35/384, …], the 5th-order one. step! propagated b5, so locextrap=true — the default, and the only setting anything in the repo ever uses — ran the lower-order solution of the pair.

Told apart by the quadrature order conditions Σᵢ bᵢcᵢᵖ = 1/(p+1) in exact rational arithmetic: [35/384, …] satisfies them through p = 4; [5179/57600, …] gives Σbc⁴ = 0.19974 ≠ 1/5, so it cannot be 5th order. [35/384, …] is also the final Butcher row B[6] and the σ=1 value of the dense-output polynomial — both consequences of DOPRI5 being FSAL, neither true of the other vector — and it is what scipy's RK45 and MATLAB's ode45 propagate.

The commit swaps the values so the names mean what they say and selects between them in step!. errest is unchanged (b5 - b4 became b4 - b5: the same numbers, still (4th) − (5th)), so the error estimate and the PI controller's 1/5 exponent are untouched. Both branches now form yn explicitly rather than one of them relying on the stage-6 accumulation evaluate! happens to leave behind, so the result no longer depends on the RHS leaving its input array alone.

Two things fall out. FSAL becomes exact — stage 7 is the RHS evaluated at the solution that is now actually propagated. And the dense output becomes a continuous extension of the propagated solution: at the endpoint the interpolant reproduces yn to round-off instead of differing by the full embedded error estimate.

Measurements

Fixed-step model problem with an exact solution (y' = i(a + |y|²)y, h = 0.1 and 0.025):

	master	after 7747c93	after 1d7e4c3
dense-output order (plain / precon)	2.00 / 2.00	4.90 / 4.99	4.90 / 4.99
interpolant slope error at step start	4.0e-1 / 2.5e-2	2.2e-5 / 2.0e-7	2.0e-7 / 2.0e-7
endpoint order	4.22 / 4.01	4.22 / 4.01	5.00 / 5.29
FSAL residual |ks[1] − f(t,y)|	4.0e-1	2.3e-5	0
interpolant vs yn at the endpoint	1.4e-2	2.1e-5	1.9e-14

N=5 soliton problem at rtol=1e-8 (exact conserved energy 1638400):

• 7747c93: stepped trajectory bit-identical — step count and hashes of the z, dz and final-field sequences all unchanged, both steppers. The summed energy of 51 interpolated saves moves from 8.35581116e7 to 8.35584000e7, against an exact 8.3558400113e7.
• 1d7e4c3: endpoint energy error 2.2e-3 → 1.1e-4 in 21944 → 14328 steps (plain) and 4.6e-1 → 7.6e-3 in 5426 steps (precon). 20–60× more accurate at the same tolerance, and 35% fewer steps for the plain stepper — the old FSAL inconsistency was inflating the error estimate too.

Tests

162 new lines in test/test_rk45.jl, on a 1-element problem that is exact for both steppers and for the preconditioned split (|y| is conserved), with min_dt == max_dt == h pinning every step:

• white-box: ks[1] must survive a step! and change on the next one;
• the interpolant's slope at the step start must be f(y(tₙ)), not f at the far end;
• convergence order of the pure dense-output error, measured against the exact solution continued from each step's own start so the trajectory's accumulated error can't mask it;
• the order conditions above, in exact rational arithmetic, pinning which vector is which permanently;
• endpoint convergence order, FSAL exactness, and interpolant/endpoint continuity;
• guards that the steps really were pinned to h, plus a repeat at rtol=1e-14 where every step is accepted only by the min_dt clamp, exercising the steplims!-forces-ok=true path the deferred move must preserve.

Reviewer notes

• Full suite: 4351 pass, 12 broken (pre-existing), 0 failures, all 31 test files. No literal needed re-baselining, including test_gnlse.jl's argmax comparisons at rtol=1e-14.
• Verified for the preconditioned path specifically (the one Luna.run uses): all kᵢ live in the frame anchored at the step start; the deferred copy stays ahead of prop!(ks[1], t, tn) so it is still rebased in the old frame; prop!_maybe and interpolate's trailing prop!(out, t, ti) still transform saved output to the lab frame on every path.
• The window/filter surgery in Luna.run's stepfun (Eω .*= ωwin, twin) mutates s.yn in place after each step, so the FSAL k1 is evaluated at the un-windowed endpoint. Pre-existing and unchanged by either commit — commit 2 slightly reduces the residual inconsistency.
• #2 supersedes the rationale for step_on only partly: the dense output is still 4th order while the endpoints are now 5th, so landing weak-signal saves on step endpoints is still worthwhile.
• Branched off master at e2d1a3d; the three commits upstream has added since touch Polarisation.jl, Scans.jl and the workflows only, so there is no overlap. The full suite above was run at that base.

🤖 Generated with Claude Code

[Copilot]

Pull request overview

Fixes two independent correctness defects in the Dormand–Prince RK45 (DOPRI5) integrator: (1) FSAL state was being moved too early and corrupting dense output, and (2) the 5th/4th order weight vectors were swapped, causing locextrap=true (the default) to propagate the lower-order solution.

Changes:

• Defer the FSAL k7 → k1 move to the start of the next step! call (guarded by the previous step’s ok) so dense output is evaluated with the correct stage data.
• Correct the DOPRI5 weight-vector labeling (b5 = 5th order, b4 = embedded 4th order) and explicitly form the propagated solution based on locextrap.
• Add targeted regression tests covering dense-output order, interpolant slope at step start, FSAL move timing/exactness, and rational order-condition checks.

Reviewed changes

Copilot reviewed 3 out of 3 changed files in this pull request and generated no comments.

File	Description
test/test_rk45.jl	Adds regression tests that detect FSAL/dense-output corruption and permanently pin the correct b-vector ordering and convergence properties.
src/RK45.jl	Defers the FSAL move until after dense output is consumed and selects the propagated solution via the corrected b5/b4 weights.
src/dopri.jl	Swaps and documents the correct 5th/4th order Dormand–Prince weights and updates the error-estimate vector accordingly.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[jtravs]

@chrisbrahms this is quite a significant bug. I doubt it has affected much of our work, but it is probably what underlines several other issues, such as #203 It came up on my ModelPNPS code as it was this that was causing weird errors when not landing exactly on z grid points. This PR seems to fix everything, but does need careful review (I am doing additonal offline testing).

[chrisbrahms]

Huh, good find. Sort of surprised this ever worked then.

I had a quick go at repeating one of my current scans because it's quick, and found negligible differences. But that's just one example.

Do you think this could be the reason that DiffEq performed so much better in your earlier tests (#400)?

[jtravs]

Yes, I think it is a very good chance.

This also makes little difference for me on general propagation. The interpolated error was still order 2, and often we had small step sizes anyway. Claude noticed it because of this pathological situation where I am sensitive to a very weak field in the presence of very strong ones on the TG-FROG sims. Even then it is minor. But it is better to be correct.

[jtravs]

That being said, while the reasoning is pretty clear, this was mostly Claude's work, so I think we need to do more manual sanity checks before merging. I am testing with more aggressive errors (higher order solitons etc.). I haven't looked up the original DOPRI5 tables though.

[chrisbrahms]

Yeah I somehow turned them around originally. From Dormand and Prince:

So the qth order is the higher-order solution (here q=5, p=4)

In their notation, a hat indicates the qth-order formula

The coefficients for that solution are b^hat

And b^hat is:

[jtravs]

Example with fundamental soliton:

A big improvement. I suspect we don't normally care about such small errors, but this comes free from performance impact.

[chrisbrahms]

LGTM and ready to merge. Can you add a brief summary to #415 ?