Add Ralston4 explicit Runge–Kutta method to OrdinaryDiffEqLowOrderRK

Author: Harsh Singh

lib/OrdinaryDiffEqLowOrderRK/src/OrdinaryDiffEqLowOrderRK.jl

@@ -45,6 +45,7 @@ export Euler, SplitEuler, Heun, Ralston, Midpoint, RK4,
     BS3, OwrenZen3, OwrenZen4, OwrenZen5, BS5,
     DP5, Anas5, RKO65, FRK65, RKM, MSRK5, MSRK6,
     PSRK4p7q6, PSRK3p5q4, PSRK3p6q5, Stepanov5, SIR54,
-    Alshina2, Alshina3, Alshina6, AutoDP5
+    Alshina2, Alshina3, Alshina6, AutoDP5,
+    Ralston4
 
 end

lib/OrdinaryDiffEqLowOrderRK/src/alg_utils.jl

@@ -14,6 +14,7 @@ alg_order(alg::RKO65) = 5
 alg_order(alg::FRK65) = 6
 alg_order(alg::RK4) = 4
 alg_order(alg::RKM) = 4
+alg_order(alg::Ralston4) = 4
 alg_order(alg::MSRK5) = 5
 alg_order(alg::MSRK6) = 6
 alg_order(alg::PSRK4p7q6) = 4

lib/OrdinaryDiffEqLowOrderRK/src/algorithms.jl

@@ -348,6 +348,28 @@ function RKM(stage_limiter!, step_limiter! = trivial_limiter!)
     return RKM(stage_limiter!, step_limiter!, False())
 end
 
+@doc explicit_rk_docstring(
+    "4th order Ralston method with minimum truncation error for 4-stage explicit Runge-Kutta.",
+    "Ralston4",
+    references = "@article{ralston1962runge,
+    title={Runge-Kutta methods with minimum error bounds},
+    author={Ralston, Anthony},
+    journal={Mathematics of Computation},
+    volume={16},
+    number={80},
+    pages={431--437},
+    year={1962}
+    }"
+)
+Base.@kwdef struct Ralston4{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+end
+function Ralston4(stage_limiter!, step_limiter! = trivial_limiter!)
+    return Ralston4(stage_limiter!, step_limiter!, False())
+end
+
 @doc explicit_rk_docstring(
     "5th order method.", "MSRK5",
     references = "Misha Stepanov - https://arxiv.org/pdf/2202.08443.pdf : Figure 3."

lib/OrdinaryDiffEqLowOrderRK/src/low_order_rk_caches.jl

@@ -1665,3 +1665,47 @@ function alg_cache(
         alg.step_limiter!, alg.thread
     )
 end
+
+@cache struct Ralston4Cache{uType, rateType, TabType, StageLimiter, StepLimiter, Thread} <:
+    OrdinaryDiffEqMutableCache
+    u::uType
+    uprev::uType
+    fsalfirst::rateType
+    k2::rateType
+    k3::rateType
+    k4::rateType
+    k::rateType
+    tmp::uType
+    tab::TabType
+    stage_limiter!::StageLimiter
+    step_limiter!::StepLimiter
+    thread::Thread
+end
+
+function alg_cache(
+        alg::Ralston4, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{false}, verbose
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    return Ralston4ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+end
+
+function alg_cache(
+        alg::Ralston4, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{true}, verbose
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    tab = Ralston4ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+    fsalfirst = zero(rate_prototype)
+    k2 = zero(rate_prototype)
+    k3 = zero(rate_prototype)
+    k4 = zero(rate_prototype)
+    k = zero(rate_prototype)
+    tmp = zero(u)
+    return Ralston4Cache(
+        u, uprev, fsalfirst, k2, k3, k4, k, tmp, tab, alg.stage_limiter!,
+        alg.step_limiter!, alg.thread
+    )
+end

lib/OrdinaryDiffEqLowOrderRK/src/low_order_rk_perform_step.jl

@@ -2209,3 +2209,65 @@ function perform_step!(integrator, cache::Alshina6Cache, repeat_step = false)
     integrator.fsallast = k7
     return nothing
 end
+
+function initialize!(integrator, cache::Ralston4ConstantCache)
+    integrator.kshortsize = 2
+    integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
+    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+    integrator.fsallast = zero(integrator.fsalfirst)
+    integrator.k[1] = integrator.fsalfirst
+    return integrator.k[2] = integrator.fsallast
+end
+
+@muladd function perform_step!(
+        integrator, cache::Ralston4ConstantCache, repeat_step = false
+    )
+    (; t, dt, uprev, u, f, p) = integrator
+    (; a21, a31, a32, a41, a42, a43, c2, c3, b1, b2, b3, b4) = cache
+    k1 = integrator.fsalfirst
+    k2 = f(uprev + dt * a21 * k1, p, t + c2 * dt)
+    k3 = f(uprev + dt * (a31 * k1 + a32 * k2), p, t + c3 * dt)
+    k4 = f(uprev + dt * (a41 * k1 + a42 * k2 + a43 * k3), p, t + dt)
+    u = uprev + dt * (b1 * k1 + b2 * k2 + b3 * k3 + b4 * k4)
+    integrator.fsallast = f(u, p, t + dt)
+    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 4)
+    integrator.k[1] = integrator.fsalfirst
+    integrator.k[2] = integrator.fsallast
+    integrator.u = u
+end
+
+get_fsalfirstlast(cache::Ralston4Cache, u) = (cache.fsalfirst, cache.k)
+function initialize!(integrator, cache::Ralston4Cache)
+    (; fsalfirst, k) = cache
+    integrator.fsalfirst = fsalfirst
+    integrator.fsallast = k
+    integrator.kshortsize = 2
+    resize!(integrator.k, integrator.kshortsize)
+    integrator.k[1] = integrator.fsalfirst
+    integrator.k[2] = integrator.fsallast
+    integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
+    return OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+end
+
+@muladd function perform_step!(integrator, cache::Ralston4Cache, repeat_step = false)
+    (; t, dt, uprev, u, f, p) = integrator
+    (; fsalfirst, k2, k3, k4, k, tmp, tab, stage_limiter!, step_limiter!, thread) = cache
+    (; a21, a31, a32, a41, a42, a43, c2, c3, b1, b2, b3, b4) = tab
+    k1 = fsalfirst
+    @.. broadcast = false thread = thread tmp = uprev + dt * a21 * k1
+    stage_limiter!(tmp, integrator, p, t + c2 * dt)
+    f(k2, tmp, p, t + c2 * dt)
+    @.. broadcast = false thread = thread tmp = uprev + dt * (a31 * k1 + a32 * k2)
+    stage_limiter!(tmp, integrator, p, t + c3 * dt)
+    f(k3, tmp, p, t + c3 * dt)
+    @.. broadcast = false thread = thread tmp = uprev + dt * (a41 * k1 + a42 * k2 + a43 * k3)
+    stage_limiter!(tmp, integrator, p, t + dt)
+    f(k4, tmp, p, t + dt)
+    @.. broadcast = false thread = thread u = uprev + dt * (b1 * k1 + b2 * k2 + b3 * k3 + b4 * k4)
+    stage_limiter!(u, integrator, p, t + dt)
+    step_limiter!(u, integrator, p, t + dt)
+    f(k, u, p, t + dt)
+    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 4)
+    return nothing
+end

lib/OrdinaryDiffEqLowOrderRK/src/low_order_rk_tableaus.jl

@@ -1964,3 +1964,37 @@ function Alshina6ConstantCache(T, T2)
         b1, b5, b6, b7, c2, c3, c4, c5, c6, c7
     )
 end
+
+struct Ralston4ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache
+    a21::T
+    a31::T
+    a32::T
+    a41::T
+    a42::T
+    a43::T
+    c2::T2
+    c3::T2
+    b1::T
+    b2::T
+    b3::T
+    b4::T
+end
+
+function Ralston4ConstantCache(T, T2)
+    s5 = sqrt(convert(T, 5))
+    b1 = (263 + 24 * s5) / 1812
+    b4 = (30 - 4 * s5) / 123
+    b2 = (125 - 1000 * s5) / 3828
+    b3 = 1 - b1 - b2 - b4
+    c2 = convert(T2, 2 // 5)
+    c3 = convert(T2, (14 - 3 * s5) / 16)
+    a21 = convert(T, c2)
+    Ac3 = 2 / (3 * b3 * (2 + 3 * s5))
+    a32 = (5 // 2) * Ac3
+    a31 = c3 - a32
+    a43 = b3 * (2 + 3 * s5) / (16 * b4)
+    Ac4 = (1 // 6 - b3 * Ac3) / b4
+    a42 = (5 // 2) * (Ac4 - a43 * c3)
+    a41 = 1 - a42 - a43
+    return Ralston4ConstantCache(a21, a31, a32, a41, a42, a43, c2, c3, b1, b2, b3, b4)
+end

src/OrdinaryDiffEq.jl

@@ -204,12 +204,14 @@ using OrdinaryDiffEqLowOrderRK: Euler, SplitEuler, Heun, Ralston, Midpoint, RK4,
     BS3, OwrenZen3, OwrenZen4, OwrenZen5, BS5,
     DP5, Anas5, RKO65, FRK65, RKM, MSRK5, MSRK6,
     PSRK4p7q6, PSRK3p5q4, PSRK3p6q5, Stepanov5, SIR54,
-    Alshina2, Alshina3, Alshina6, AutoDP5
+    Alshina2, Alshina3, Alshina6, AutoDP5,
+    Ralston4
 export OrdinaryDiffEqLowOrderRK, Euler, SplitEuler, Heun, Ralston, Midpoint, RK4,
     BS3, OwrenZen3, OwrenZen4, OwrenZen5, BS5,
     DP5, Anas5, RKO65, FRK65, RKM, MSRK5, MSRK6,
     PSRK4p7q6, PSRK3p5q4, PSRK3p6q5, Stepanov5, SIR54,
-    Alshina2, Alshina3, Alshina6, AutoDP5
+    Alshina2, Alshina3, Alshina6, AutoDP5,
+    Ralston4
 
 import OrdinaryDiffEqFunctionMap: OrdinaryDiffEqFunctionMap
 using OrdinaryDiffEqFunctionMap: FunctionMap

test/integrators/step_limiter_test.jl

@@ -29,7 +29,7 @@ end
         Rosenbrock23, Rosenbrock32, ROS3P, Rodas3, Rodas23W, Rodas3P, Rodas4, Rodas42,
         Rodas4P, Rodas4P2, Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr,
         AdaptiveRadau, RadauIIA9, RadauIIA5, RadauIIA3, SIR54,
-        Alshina2, Alshina3, Heun, Ralston, Midpoint, RK4,
+        Alshina2, Alshina3, Heun, Ralston, Ralston4, Midpoint, RK4,
         OwrenZen3, OwrenZen4, OwrenZen5,
         BS3, DP5, Tsit5, DP8, TanYam7, TsitPap8, FRK65, PFRK87, BS5, Vern6, Vern7,
         Vern8, Vern9, QPRK98, SSPRKMSVS43, SSPRKMSVS32, SSPRK432, SSPRK43,

test/interface/float32.jl

@@ -11,7 +11,7 @@ tspan = (zero(NF), NF(1.0e3)) #integrate from t=0 to t = 1000
 ode_prob = ODEProblem(some_arbitrary_function!, u, tspan, params)
 
 for alg in [
-        Euler(), Midpoint(), Heun(), Ralston(), RK4(), SSPRK104(), SSPRK22(), SSPRK33(),
+        Euler(), Midpoint(), Heun(), Ralston(), Ralston4(), RK4(), SSPRK104(), SSPRK22(), SSPRK33(),
         SSPRK43(), SSPRK432(), BS3(), BS5(), DP5(), DP8(), Feagin10(), Feagin12(),
         Feagin14(), TanYam7(), Tsit5(), TsitPap8(), Vern6(), Vern7(), Vern8(), Vern9(),
     ]

test/interface/units_tests.jl

@@ -3,7 +3,7 @@ using LinearAlgebra, Test, ADTypes
 
 @testset "Algorithms" begin
     algs = [
-        Euler(), Midpoint(), Heun(), Ralston(), RK4(), SSPRK104(), SSPRK22(), SSPRK33(),
+        Euler(), Midpoint(), Heun(), Ralston(), Ralston4(), RK4(), SSPRK104(), SSPRK22(), SSPRK33(),
         SSPRK43(), SSPRK432(), BS3(), BS5(), DP5(), DP8(), Feagin10(), Feagin12(),
         Feagin14(), TanYam7(), Tsit5(), TsitPap8(), Vern6(), Vern7(), Vern8(), Vern9(),
     ]