diff --git a/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl b/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl index ca6328e98c..162f8bd514 100644 --- a/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl +++ b/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl @@ -80,7 +80,7 @@ PrecompileTools.@compile_workload begin solver_list = nothing end -export Vern6, Vern7, Vern8, Vern9 -export AutoVern6, AutoVern7, AutoVern8, AutoVern9 +export Vern6, Vern7, Vern8, Vern9, RKV87 +export AutoVern6, AutoVern7, AutoVern8, AutoVern9, RKV87 end diff --git a/lib/OrdinaryDiffEqVerner/src/alg_utils.jl b/lib/OrdinaryDiffEqVerner/src/alg_utils.jl index f38d2ec768..50b54947bf 100644 --- a/lib/OrdinaryDiffEqVerner/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqVerner/src/alg_utils.jl @@ -12,4 +12,10 @@ alg_stability_size(alg::Vern7) = 4.6400 alg_stability_size(alg::Vern8) = 5.8641 alg_stability_size(alg::Vern9) = 4.4762 +# ── Verner RKV87-IIa traits ──────────────────────────────────────── +alg_order(::RKV87) = 8 # primary solution order +embedded_order(::RKV87) = 7 # will matter once adaptivity is added +isfsal(::RKV87) = false # pair is NOT FSAL +has_dense_output(::RKV87)= false # we haven’t wired the interpolant yet + SciMLBase.has_lazy_interpolation(alg::Union{Vern6, Vern7, Vern8, Vern9}) = true diff --git a/lib/OrdinaryDiffEqVerner/src/algorithms.jl b/lib/OrdinaryDiffEqVerner/src/algorithms.jl index e1fbb1fb16..6ee36dd9e1 100644 --- a/lib/OrdinaryDiffEqVerner/src/algorithms.jl +++ b/lib/OrdinaryDiffEqVerner/src/algorithms.jl @@ -85,6 +85,27 @@ function Vern8(stage_limiter!, step_limiter! = trivial_limiter!; lazy = true) Vern8(stage_limiter!, step_limiter!, False(), lazy) end +# ------------------------------------------------------------------ +# Verner RKV87-IIa : 13-stage explicit 8 (7) Runge–Kutta pair +# Fixed-step implementation (no adaptivity or dense output yet) +# ------------------------------------------------------------------ + +@doc explicit_rk_docstring( + "Verner’s “efficient” 13-stage 8/7 Runge–Kutta method (RKV87-IIa).\n" * + "This variant is wired **for fixed time steps**. Use it via\n\n" * + "```julia\nsol = solve(prob, RKV87(), dt = h)\n```", + "RKV87") +Base.@kwdef struct RKV87{StageLimiter, StepLimiter, Thread} <: + OrdinaryDiffEqConstantAlgorithm + stage_limiter!::StageLimiter = trivial_limiter! + step_limiter!::StepLimiter = trivial_limiter! + thread::Thread = False() +end +TruncatedStacktraces.@truncate_stacktrace RKV87 3 + +# Register the algorithm with the dispatcher +add_algorithm(RKV87(), RKV87ConstantCache, Vern87IIaTableau) + @doc explicit_rk_docstring( "Verner's “Most Efficient” 9/8 Runge-Kutta method. (lazy9th order interpolant).", "Vern9", @@ -113,6 +134,7 @@ function Vern9(stage_limiter!, step_limiter! = trivial_limiter!; lazy = true) Vern9(stage_limiter!, step_limiter!, False(), lazy) end + """ Automatic switching algorithm that can switch between the (non-stiff) `Vern6()` and `stiff_alg`. @@ -153,3 +175,5 @@ To gain access to stiff algorithms you might have to install additional librarie such as `OrdinaryDiffEqRosenbrock`. """ AutoVern9(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern9(lazy = lazy), alg; kwargs...) + + diff --git a/lib/OrdinaryDiffEqVerner/src/verner_caches.jl b/lib/OrdinaryDiffEqVerner/src/verner_caches.jl index 9a3444d101..e87c0b6974 100644 --- a/lib/OrdinaryDiffEqVerner/src/verner_caches.jl +++ b/lib/OrdinaryDiffEqVerner/src/verner_caches.jl @@ -263,9 +263,12 @@ struct Vern9ConstantCache <: OrdinaryDiffEqConstantCache lazy::Bool end -function alg_cache(alg::Vern9, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - Vern9ConstantCache(alg.lazy) +# ------------------------------------------------------------------ +# Cache for Verner RKV87-IIa (13 stage values + 1 work vector) +# ------------------------------------------------------------------ +@cache RKV87ConstantCache begin + k1::K ; k2::K ; k3::K ; k4::K ; k5::K ; k6::K ; k7::K + k8::K ; k9::K ; k10::K ; k11::K ; k12::K ; k13::K + tmp::K # work vector for low-storage kernel end + diff --git a/lib/OrdinaryDiffEqVerner/src/verner_rk_perform_step.jl b/lib/OrdinaryDiffEqVerner/src/verner_rk_perform_step.jl index 898be5af5d..55aa0b29c3 100644 --- a/lib/OrdinaryDiffEqVerner/src/verner_rk_perform_step.jl +++ b/lib/OrdinaryDiffEqVerner/src/verner_rk_perform_step.jl @@ -1013,6 +1013,18 @@ end end end +# ------------------------------------------------------------------ +# RKV87-IIa wrappers +# ------------------------------------------------------------------ + +initialize!(integrator, cache::RKV87ConstantCache) = nothing + +# The tableau we built is an ExplicitRKTableau with 13 stages +perform_step!(integrator, + cache::RKV87ConstantCache, + tab::ExplicitRKTableau{13}) = + verner_lowstorage_step!(integrator, cache, tab) + function initialize!(integrator, cache::Vern9Cache) @unpack k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15, k16 = cache @unpack k = integrator diff --git a/lib/OrdinaryDiffEqVerner/src/verner_tableaus.jl b/lib/OrdinaryDiffEqVerner/src/verner_tableaus.jl index 5edc4f4a10..f93a7780b5 100644 --- a/lib/OrdinaryDiffEqVerner/src/verner_tableaus.jl +++ b/lib/OrdinaryDiffEqVerner/src/verner_tableaus.jl @@ -1,3 +1,5 @@ +using StaticArrays +module OrdinaryDiffEqVerner ## Vern6 struct Vern6ExtraStages{T, T2} c10::T2 @@ -2306,6 +2308,83 @@ function Vern8Tableau(T, T2) extra, interp) end +# -------------------------------------------------------------------- +# Verner RKV87-IIa : 13-stage explicit 8 (7) pair (efficient family) +# -------------------------------------------------------------------- +function Vern87IIaTableau(::Type{T}, ::Type{T2}) where {T, T2} + c = @SVector T[ + 0.0, + 0.092662, + 0.131223036175, + 0.196834554263, + 0.427173, + 0.485972, + 0.161915, + 0.985468, + 0.96269773486, + 0.99626, + 0.997947, + 1.0, + 1.0 + ] + + A = @SMatrix T[ + 0 0 0 0 0 0 0 0 0 0 0 0 0; + 0.092662 0 0 0 0 0 0 0 0 0 0 0 0; + 0.0383074654825 0.0929155706929 0 0 0 0 0 0 0 0 0 0 0; + 0.0492086385658 0 0.147625915697 0 0 0 0 0 0 0 0 0 0; + 0.27430760857 0 -0.93198872031 1.08485411174 0 0 0 0 0 0 0 0 0; + 0.0646185297094 0 0 0.268762921337 0.152590548954 0 0 0 0 0 0 0 0; + 0.0718915581977 0 0 0.312986472775 0.140861770313 0.0559558792414 0 0 0 0 0 0 0; + 0.159677185234 0 0 0.693559173307 -0.316843836922 0.420825919995 -0.0847501736448 0 0 0 0 0 0; + 0.0591917679674 0 0 0 -2.23377950868 0.390944379945 -0.936828383822 58.4572000999 0.00589430972373 0 0 0 0; + -6.68986189932 0 0 -81.4427100405 13.36778825698 -4.47077763842 80.2332139216 -0.0131363833621 0.0117437830322 0 0 0 0; + -6.7888419558 0 0 -82.6563985593 13.5997392187 -4.65073689329 82.1513831957 -0.0135453129954 0.0121304908312 -0.000312288856455 0 0 0; + 0.0514137795362 0 0 0 0 0.318131086582 0.458946584938 -0.400240655913 0.0569851840034 -0.0010554370029 0.0000499398326837 0 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 + ] + + b = @SVector T[ + 0.0462554315971, + 0.0, + 0.0, + 0.0, + 0.0, + 0.370666616552, + 0.259044082455, + -679.984146817, + 49.8916112904, + 10271.2352221373, + -14782.1966063569, + 5141.37795361606, + 0.0 + ] + + b̂ = @SVector T[ # embedded order-7 weights + 0.0463850423437, + 0.0, + 0.0, + 0.0, + 0.0, + 0.372576768158, + 0.258568549512, + -147.495076759, + 23.8436271264, + 347.426416673, + 0.0, + 0.0, + -223.452497401 + ] + + return ExplicitRKTableau(c, A, b, b̂) +end + +# Fast path for CompiledFloats (optional; fine to add later) +function Vern87IIaTableau(::Type{<:CompiledFloats}, ::Type{<:CompiledFloats}) + Vern87IIaTableau(Float64, Float64) +end + + ## Vern9 struct Vern9ExtraStages{T, T2} c17::T2 @@ -3891,4 +3970,82 @@ end a1610, a1611, a1612, a1613, b1, b8, b9, b10, b11, b12, b13, b14, b15, btilde1, btilde8, btilde9, btilde10, btilde11, btilde12, btilde13, btilde14, btilde15, btilde16) + + +end + +# -------------------------------------------------------------------- +# Verner RKV87-IIa : 13-stage explicit 8 (7) pair (efficient family) +# -------------------------------------------------------------------- +function Vern87IIaTableau(::Type{T}, ::Type{T2}) where {T, T2} + c = @SVector T[ + 0.0, + 0.092662, + 0.131223036175, + 0.196834554263, + 0.427173, + 0.485972, + 0.161915, + 0.985468, + 0.96269773486, + 0.99626, + 0.997947, + 1.0, + 1.0 + ] + + A = @SMatrix T[ + 0 0 0 0 0 0 0 0 0 0 0 0 0; + 0.092662 0 0 0 0 0 0 0 0 0 0 0 0; + 0.0383074654825 0.0929155706929 0 0 0 0 0 0 0 0 0 0 0; + 0.0492086385658 0 0.147625915697 0 0 0 0 0 0 0 0 0 0; + 0.27430760857 0 -0.93198872031 1.08485411174 0 0 0 0 0 0 0 0 0; + 0.0646185297094 0 0 0.268762921337 0.152590548954 0 0 0 0 0 0 0 0; + 0.0718915581977 0 0 0.312986472775 0.140861770313 0.0559558792414 0 0 0 0 0 0 0; + 0.159677185234 0 0 0.693559173307 -0.316843836922 0.420825919995 -0.0847501736448 0 0 0 0 0 0; + 0.0591917679674 0 0 0 -2.23377950868 0.390944379945 -0.936828383822 58.4572000999 0.00589430972373 0 0 0 0; + -6.68986189932 0 0 -81.4427100405 13.36778825698 -4.47077763842 80.2332139216 -0.0131363833621 0.0117437830322 0 0 0 0; + -6.7888419558 0 0 -82.6563985593 13.5997392187 -4.65073689329 82.1513831957 -0.0135453129954 0.0121304908312 -0.000312288856455 0 0 0; + 0.0514137795362 0 0 0 0 0.318131086582 0.458946584938 -0.400240655913 0.0569851840034 -0.0010554370029 0.0000499398326837 0 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 + ] + + b = @SVector T[ + 0.0462554315971, + 0.0, + 0.0, + 0.0, + 0.0, + 0.370666616552, + 0.259044082455, + -679.984146817, + 49.8916112904, + 10271.2352221373, + -14782.1966063569, + 5141.37795361606, + 0.0 + ] + + b̂ = @SVector T[ # embedded order-7 weights + 0.0463850423437, + 0.0, + 0.0, + 0.0, + 0.0, + 0.372576768158, + 0.258568549512, + -147.495076759, + 23.8436271264, + 347.426416673, + 0.0, + 0.0, + -223.452497401 + ] + + return ExplicitRKTableau(c, A, b, b̂) +end + +# Fast path for CompiledFloats (optional; fine to add later) +function Vern87IIaTableau(::Type{<:CompiledFloats}, ::Type{<:CompiledFloats}) + Vern87IIaTableau(Float64, Float64) end