Add arbitrary-order GaussLegendre FIRK method (draft) by Sreeram-Shankar · Pull Request #3470 · SciML/OrdinaryDiffEq.jl

Sreeram-Shankar · 2026-04-17T23:27:32Z

Add arbitrary-order Gauss-Legendre FIRK method

This is a draft implementation of arbitrary-order Gauss-Legendre
collocation methods as requested in the linked issue. With s stages
the method has classical order 2s and is symplectic, making it
suitable for Hamiltonian systems and problems requiring long-time
geometric integration.

What's implemented

GaussLegendre(num_stages=s) algorithm struct following the existing
Radau interface pattern
Tableau generation using FastGaussQuadrature.gausslegendre with
Lagrange basis integration for the A matrix
Tableau cache for Float64 at stages 2-5
Full Newton iteration on the real-valued sn × sn system
Basic embedded error estimate using s-1 GL quadrature weights for
adaptive stepping (placeholder, needs improvement)
Fixed dt works correctly and produces accurate solutions

Known issues / TODOs

Adaptive dt is overly conservative due to the placeholder error estimate
No stage decoupling yet, I've looked at the approach in the Antonana, Makazaga,
and Murua paper but haven't implemented it yet
No dense output / interpolation
W matrix construction could be more efficient

Testing

Basic out-of-place solve works correctly:

f(u, p, t) = [u[2], -u[1]]
prob = ODEProblem(f, [1.0, 0.0], (0.0, 10.0))
sol = solve(prob, GaussLegendre(num_stages=2), dt=0.1, adaptive=false)
# sol.u[end] ≈ [cos(10), -sin(10)] ✓

Happy to iterate on any of this based on feedback.

oscardssmith · 2026-04-17T23:48:12Z

Needs some tests, but otherwise, looks reasonable.

ranocha · 2026-04-18T03:26:12Z

+        end
+        e = b .- b_low
+    else
+        e = b


It looks like for s = 1, no real embedded method is implemented. Should an adaptive time integration for s = 1 be checked during initialization to throw an error in this case?

we probably should make s=2 be the minimum (otherwise we just have a trapezoid method)

If we use the general (all s) implementation, it could be nice to still allow s = 1 to compare the performance with the ImplicitMidpoint method.

That's fair. I think the lack of a natural embedded method may make integrating it difficult, but it is likely worth a try.

I changed the error estimator to a Richardson controller for now since embedding error estimators into the method in the way I tried requires extra stages and is nontrivial. I can look into natural embedded pairings for GL methods and implement that if you have any suggestions.

using Gauss-Lobatto quadrature to get 2s-2 order LobattoIIIA methods, which share the same nodes as GL plus the endpoints, could be an option

I think it would also be totally fine to just implement Gauss methods as non-adaptive methods for now.

Yes, one step to get them converging at the right order, then another to make them adaptive, is a very common thing to make the PRs saner to review.

Got it, I’ll scope this around the fixed-step Gauss methods so the review stays focused on correctness/order. DID also add convergence/order tests with some basic coverage for the Richardson controller.

ranocha · 2026-04-18T03:33:11Z

+    A = Matrix{T1}(undef, num_stages, num_stages)
+    for i in 1:num_stages
+        for j in 1:num_stages
+            nodes, weights = gausslegendre(2 * num_stages)


As far as I know, FastGaussQuadrature.jl only works with Float64. Would it make sense to implement the analytical coefficients for the known cases when higher precision is required?

With JuliaApproximation/FastGaussQuadrature.jl#149 that concern will be resolved.

ChrisRackauckas · 2026-04-19T11:27:48Z

Is there an advantage to having this in the same subpackage as the Radau methods?

oscardssmith · 2026-04-19T13:10:46Z

we could split them, but this is a FIRK method, it uses all the same deps, and this is only a single algorithm so the include and image size change should be minimal. IMO this seems like a good place for it.

Sreeram-Shankar · 2026-04-19T14:18:13Z

I’ve added convergence/order tests (checking 2s order for s = 2, 3 and a higher-order accuracy check for s = 4), plus some basic adaptive tests for the Richardson controller. Also integrated them into the existing FIRK problem tests and allocation tests.

Is there anything else you’d like to see test-wise? I was thinking about symplectic or A-stability checks.

Sreeram-Shankar · 2026-04-19T18:51:50Z

formatted with runic

…air (which can be brought back later) with a half-step based controller for now

…calls

oscardssmith · 2026-04-24T05:56:11Z

+        linsolve = nothing, precs = nothing,
+        extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
+        new_W_γdt_cutoff = 1 // 5,
+        controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,


we probably should use a PI controller until we get order adaptivity

oscardssmith · 2026-04-24T05:57:13Z

+    extra_keyword_description = extra_keyword_description,
+    extra_keyword_default = extra_keyword_default
+)
+struct GaussLegendre{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:


the type params here are missing the v7 updates

oscardssmith self-requested a review April 17, 2026 23:48

ranocha reviewed Apr 18, 2026

View reviewed changes

Sreeram-Shankar added 6 commits April 23, 2026 20:39

Fix merge conflicts and update for gauss methods

4baec54

Added experimental tag

9d6e918

Added Richardson based step size adaptivity, replacing the embedded p…

16290aa

…air (which can be brought back later) with a half-step based controller for now

Added tests for Guass Legendre and fixed issues in cache constructor …

0f7395f

…calls

Formatted with Runic

8ddd441

Update files to fix conflicts

9649936

Sreeram-Shankar force-pushed the add-gauss-legendre branch from 2f8be28 to 9649936 Compare April 24, 2026 02:40

oscardssmith reviewed Apr 25, 2026

View reviewed changes

Uh oh!

Conversation

Sreeram-Shankar commented Apr 17, 2026