Add the Makazaga-Murua MM5GEE global-error-estimating solver

Adds the DOPRI5-based order-5 scheme with cheap global error estimation
of Makazaga and Murua (BIT 43, 2003) as another tableau in the GLEE
general linear framework (SciML/GlobalDiffEq.jl#6, the coefficient-
complete published representative of the Dormand-Duckers-Prince RK
triple family). Stages 1-7 are the standard DOPRI5 stages, three extra
stages propagate an order-6 companion solution, and a generic
solution-stage FSAL detection reuses stage 7's evaluation for fsallast,
giving exactly 9 function evaluations per step (verified by an nf
accounting test). The tableau was verified with exact rational
arithmetic against the order and independency conditions.

Verified locally on the unstable Prince42 problem: solution order 5,
estimate accuracy order 6, est/true-error ratio 1.04, 9.01 f-evals/step.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

Author: ChrisRackauckas

docs/src/globalerrorcontrol/GlobalDiffEq.md

@@ -18,9 +18,10 @@ using GlobalDiffEq
 
   - For a solution accompanied by a running, asymptotically correct estimate of
     its global error at every time point, use the global-error-estimating
-    solvers [`GLEE24`](@ref) or [`GLEE35`](@ref) (Constantinescu 2016). These
-    cost only a few extra stages per step over a plain method of the same order
-    and require nothing beyond the right-hand side `f`.
+    solvers [`GLEE24`](@ref), [`GLEE35`](@ref) (Constantinescu 2016), or the
+    Dormand-Prince-based [`MM5GEE`](@ref) (Makazaga and Murua 2003). These cost
+    only a few extra stages per step over a plain method of the same order and
+    require nothing beyond the right-hand side `f`.
   - To *control* the endpoint global error to a tolerance `gtol`, wrap any
     adaptive solver in [`GlobalAdjoint`](@ref) (adjoint-based, for endpoint
     functionals; requires SciMLSensitivity and QuadGK to be loaded). It solves
@@ -51,6 +52,7 @@ errs = global_error_estimate(sol)  # global error estimate at every sol.t
 GLEE23
 GLEE24
 GLEE35
+MM5GEE
 global_error_estimate
 ```
 

lib/GlobalDiffEq/Project.toml

@@ -1,7 +1,7 @@
 name = "GlobalDiffEq"
 uuid = "1d72d19b-84cc-4cb7-a099-7cbdb9ccc67c"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "1.4.0"
+version = "1.5.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"

lib/GlobalDiffEq/README.md

@@ -15,6 +15,8 @@ provides:
   - `GLEE23`, `GLEE24`, `GLEE35`: explicit general linear methods with built-in
     global error estimation (Constantinescu 2016), which propagate the solution
     together with an asymptotically correct estimate of its global error.
+  - `MM5GEE`: the Makazaga-Murua (2003) Dormand-Prince-based order-5 scheme
+    with a cheap built-in global error estimate.
   - `GlobalAdjoint`: adjoint-based a posteriori endpoint error estimation and
     control (Cao and Petzold 2004), available as a package extension when
     SciMLSensitivity and QuadGK are loaded.

lib/GlobalDiffEq/src/GlobalDiffEq.jl

@@ -22,7 +22,7 @@ include("glee/caches.jl")
 include("glee/perform_step.jl")
 
 export GlobalAdjoint, GlobalRichardson, adjoint_error_estimate
-export GLEE23, GLEE24, GLEE35, global_error_estimate
+export GLEE23, GLEE24, GLEE35, MM5GEE, global_error_estimate
 
 @setup_workload begin
     # Simple test ODE: exponential decay du/dt = -u

lib/GlobalDiffEq/src/glee/algorithms.jl

@@ -58,13 +58,33 @@ $(_GLEE_DOCS_SHARED)
 """
 struct GLEE35 <: AbstractGLEEAlgorithm end
 
+"""
+    MM5GEE()
+
+Fifth-order Dormand-Prince-based scheme with cheap global error estimation
+(Makazaga and Murua, BIT Numerical Mathematics 43, 2003). Propagates the
+standard DOPRI5 solution together with an order-6 companion solution through
+three extra stages (9 function evaluations per step in total), whose
+difference is an asymptotically correct global error estimate. This is the
+recommended method of the Runge-Kutta triple family of Dormand, Duckers and
+Prince, as it is the coefficient-complete published scheme of that type.
+
+$(_GLEE_DOCS_SHARED)
+
+  - J. Makazaga and A. Murua, New Runge-Kutta based schemes for ODEs with
+    cheap global error estimation, BIT Numerical Mathematics 43 (2003).
+"""
+struct MM5GEE <: AbstractGLEEAlgorithm end
+
 SciMLBase.alg_order(::GLEE23) = 2
 SciMLBase.alg_order(::GLEE24) = 2
 SciMLBase.alg_order(::GLEE35) = 3
+SciMLBase.alg_order(::MM5GEE) = 5
 
 OrdinaryDiffEqCore.alg_adaptive_order(alg::AbstractGLEEAlgorithm) =
     SciMLBase.alg_order(alg)
 
 _glee_tableau_for(::GLEE23, ::Type{T}, ::Type{T2}) where {T, T2} = GLEE23Tableau(T, T2)
 _glee_tableau_for(::GLEE24, ::Type{T}, ::Type{T2}) where {T, T2} = GLEE24Tableau(T, T2)
 _glee_tableau_for(::GLEE35, ::Type{T}, ::Type{T2}) where {T, T2} = GLEE35Tableau(T, T2)
+_glee_tableau_for(::MM5GEE, ::Type{T}, ::Type{T2}) where {T, T2} = MM5GEETableau(T, T2)

lib/GlobalDiffEq/src/glee/tableaus.jl

@@ -142,3 +142,48 @@ function GLEE35Tableau(::Type{T}, ::Type{T2}) where {T, T2}
     Uyε, Byε = _yytilde_to_yeps(U, B, Z)
     return _glee_tableau(T, T2, A, Uyε, Byε, 0, 3)
 end
+
+"""
+    MM5GEETableau(T, T2)
+
+Tableau of the Makazaga-Murua global-error-estimating scheme built on the
+Dormand-Prince 5(4) method (Makazaga and Murua, BIT Numerical Mathematics 43,
+2003, Table 4.1), rewritten in the y-ε general linear form: stages 1-7 are the
+standard DOPRI5 stages (stage 7 the FSAL stage), stages 8-10 propagate the
+order-6 companion solution `ȳ = y + ε` through the mixing coefficients
+`U[i,2] = 1 - μ_i`, and the second output row is `b̄ - b`.
+"""
+function MM5GEETableau(::Type{T}, ::Type{T2}) where {T, T2}
+    R(num, den) = big(num) // big(den)
+    Z = R(0, 1)
+    b = [
+        R(35, 384), Z, R(500, 1113), R(125, 192), R(-2187, 6784), R(11, 84), Z,
+        Z, Z, Z,
+    ]
+    A = [
+        Z Z Z Z Z Z Z Z Z Z;
+        R(1, 5) Z Z Z Z Z Z Z Z Z;
+        R(3, 40) R(9, 40) Z Z Z Z Z Z Z Z;
+        R(44, 45) R(-56, 15) R(32, 9) Z Z Z Z Z Z Z;
+        R(19372, 6561) R(-25360, 2187) R(64448, 6561) R(-212, 729) Z Z Z Z Z Z;
+        R(9017, 3168) R(-355, 33) R(46732, 5247) R(49, 176) R(-5103, 18656) Z Z Z Z Z;
+        R(35, 384) Z R(500, 1113) R(125, 192) R(-2187, 6784) R(11, 84) Z Z Z Z;
+        R(26251126, 75292183) R(-30511879, 68834945) R(11490887, 155205387) R(700737845, 174891007) R(-5336, 941) R(5735, 1214) R(-2507, 898) Z Z Z;
+        R(-126276029, 115017392) R(153409379, 49308629) R(-107711621, 48274693) R(-675136779, 64711289) R(559269939, 36928210) R(-669687859, 52442748) R(193952703, 25738526) R(169021117, 130072535) Z Z;
+        R(89178409, 82486612) R(-275044175, 99029299) R(115406143, 68971088) R(140298385, 24130572) R(-344040692, 42025591) R(121333564, 17575013) R(-190380249, 47005513) R(-12078143, 165601005) R(56747365, 92317949) Z
+    ]
+    # U[:, 2] = 1 - μ_i: the weight of the companion register ȳ = y + ε in each
+    # stage; μ_i = 1 for the DOPRI5 stages 1-7.
+    one_minus_μ = [
+        Z, Z, Z, Z, Z, Z, Z,
+        R(140719960, 143529893), R(941, 896), R(92493035, 95359057),
+    ]
+    U = hcat(fill(R(1, 1), 10), one_minus_μ)
+    b̄ = [
+        R(56696811, 789712427), Z, R(-47431484, 279691831), R(72791025, 357831874),
+        R(17490085, 349505178), R(-66245097, 563676842), R(-24, 611),
+        R(40757463, 82884629), R(33159666, 111811519), R(42422453, 199331202),
+    ]
+    B = permutedims(hcat(b, b̄ .- b))
+    return _glee_tableau(T, T2, A, U, B, 0, 5)
+end

lib/GlobalDiffEq/test/glee_tests.jl

@@ -42,12 +42,15 @@ end
             # to the true error one order faster than the solution converges
             @test est_order ≈ p + 1 atol = 0.2
         end
+        sol_order, est_order = glee_convergence(prob, MM5GEE(), 2.0 .^ -(2:5))
+        @test sol_order ≈ 5 atol = 0.35
+        @test est_order > 5.5
     end
 end
 
 @testset "GLEE estimate tracks the true error" begin
     prob = ODEProblem(f_prince!, [0.0], prince_tspan)
-    for alg in (GLEE23(), GLEE24(), GLEE35())
+    for alg in (GLEE23(), GLEE24(), GLEE35(), MM5GEE())
         sol = solve(prob, alg; dt = 1 / 128, adaptive = false)
         err = prince_exact(prince_tspan[2]) - sol.u[end].x[1][1]
         est = global_error_estimate(sol)[end][1]
@@ -57,7 +60,7 @@ end
 
 @testset "GLEE adaptive stepping" begin
     prob = ODEProblem(f_prince!, [0.0], prince_tspan)
-    for alg in (GLEE24(), GLEE35())
+    for alg in (GLEE24(), MM5GEE())
         sol = solve(prob, alg; abstol = 1.0e-8, reltol = 1.0e-8)
         @test SciMLBase.successful_retcode(sol)
         err = prince_exact(prince_tspan[2]) - sol.u[end].x[1][1]
@@ -78,6 +81,14 @@ end
     @test global_error_estimate(sol)[end][1] / err ≈ 1 rtol = 0.1
 end
 
+@testset "MM5GEE function evaluation count" begin
+    prob = ODEProblem(f_prince!, [0.0], prince_tspan)
+    sol = solve(prob, MM5GEE(); dt = 1 / 32, adaptive = false)
+    # 9 f-evals per step (FSAL first stage + solution stage reuse) plus the
+    # initialization evaluation
+    @test sol.stats.nf == 9 * (length(sol.t) - 1) + 1
+end
+
 @testset "GLEE argument validation" begin
     prob = ODEProblem(f_prince!, [0.0], prince_tspan)
     partitioned_prob = ODEProblem(
@@ -96,4 +107,5 @@ end
     @test SciMLBase.alg_order(GLEE23()) == 2
     @test SciMLBase.alg_order(GLEE24()) == 2
     @test SciMLBase.alg_order(GLEE35()) == 3
+    @test SciMLBase.alg_order(MM5GEE()) == 5
 end