Merge pull request #1948 from HenryLangner/master

Added ```FineRKN5``` Algorithm

docs/src/dynamical/nystrom.md

@@ -5,6 +5,7 @@ IRKN3
 Nystrom4
 Nystrom4VelocityIndependent
 IRKN4
+FineRKN5
 Nystrom5VelocityIndependent
 DPRKN6
 DPRKN8

src/OrdinaryDiffEq.jl

@@ -437,7 +437,7 @@ export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog,
 
 export SplitEuler
 
-export Nystrom4, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent,
+export Nystrom4, FineRKN5, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent,
     IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7
 
 export ROCK2, ROCK4, RKC, IRKC, ESERK4, ESERK5, SERK2

src/alg_utils.jl

@@ -502,6 +502,7 @@ alg_order(alg::SofSpa10) = 10
 
 alg_order(alg::IRKN3) = 3
 alg_order(alg::Nystrom4) = 4
+alg_order(alg::FineRKN5) = 5
 alg_order(alg::Nystrom4VelocityIndependent) = 4
 alg_order(alg::IRKN4) = 4
 alg_order(alg::Nystrom5VelocityIndependent) = 5

src/algorithms.jl

@@ -5251,6 +5251,31 @@ E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equatio
 """
 struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end
 
+"""
+    FineRKN5()
+
+A 5th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. 
+Can only be used with fixed time steps.
+
+This method requires that the acceleration equation needs to be independent of the velocity
+otherwise it results in convergence order loss.
+
+## References
+```
+@article{fine1987low,
+  title={Low order practical {R}unge-{K}utta-{N}ystr{\"o}m methods},
+  author={Fine, Jerry Michael},
+  journal={Computing},
+  volume={38},
+  number={4},
+  pages={281--297},
+  year={1987},
+  publisher={Springer}
+}
+```
+"""
+struct FineRKN5 <: OrdinaryDiffEqPartitionedAlgorithm end
+
 """
     Nystrom4VelocityIdependent
 

src/caches/rkn_caches.jl

@@ -36,6 +36,47 @@ end
 
 # alg_cache(alg::Nystrom4,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = Nystrom4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
 
+@cache struct FineRKN5Cache{uType, rateType, reducedRateType, TabType} <:
+              OrdinaryDiffEqMutableCache
+    u::uType
+    uprev::uType
+    fsalfirst::rateType
+    k2::reducedRateType
+    k3::reducedRateType
+    k4::reducedRateType
+    k5::reducedRateType
+    k6::reducedRateType
+    k7::reducedRateType
+    k::rateType
+    tmp::uType
+    tab::TabType
+end
+
+function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits},
+                   ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+                   dt, reltol, p, calck,
+                   ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    reduced_rate_prototype = rate_prototype.x[2]
+    tab = FineRKN5ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+    k1 = zero(rate_prototype)
+    k2 = zero(reduced_rate_prototype)
+    k3 = zero(reduced_rate_prototype)
+    k4 = zero(reduced_rate_prototype)
+    k5 = zero(reduced_rate_prototype)
+    k6 = zero(reduced_rate_prototype)
+    k7 = zero(reduced_rate_prototype)
+    k = zero(rate_prototype)
+    tmp = zero(u)
+    FineRKN5Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k, tmp, tab)
+end
+
+function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits},
+                   ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+                   dt, reltol, p, calck,
+                   ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    FineRKN5ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+end
+
 @cache struct Nystrom4VelocityIndependentCache{uType, rateType, reducedRateType} <:
               OrdinaryDiffEqMutableCache
     u::uType

src/perform_step/rkn_perform_step.jl

@@ -4,14 +4,14 @@
 ## y₁ = y₀ + hy'₀ + h²∑b̄ᵢk'ᵢ
 ## y'₁ = y'₀ + h∑bᵢk'ᵢ
 
-const NystromCCDefaultInitialization = Union{Nystrom4ConstantCache,
-    Nystrom4VelocityIndependentConstantCache,
-    Nystrom5VelocityIndependentConstantCache,
-    IRKN3ConstantCache, IRKN4ConstantCache,
-    DPRKN4ConstantCache, DPRKN5ConstantCache,
-    DPRKN6FMConstantCache, DPRKN8ConstantCache,
-    DPRKN12ConstantCache, ERKN4ConstantCache,
-    ERKN5ConstantCache, ERKN7ConstantCache}
+const NystromCCDefaultInitialization = Union{Nystrom4ConstantCache, FineRKN5ConstantCache,
+   Nystrom4VelocityIndependentConstantCache,
+   Nystrom5VelocityIndependentConstantCache,
+   IRKN3ConstantCache, IRKN4ConstantCache,
+   DPRKN4ConstantCache, DPRKN5ConstantCache,
+   DPRKN6FMConstantCache, DPRKN8ConstantCache,
+   DPRKN12ConstantCache, ERKN4ConstantCache,
+   ERKN5ConstantCache, ERKN7ConstantCache}
 
 function initialize!(integrator, cache::NystromCCDefaultInitialization)
     integrator.kshortsize = 2
@@ -25,14 +25,14 @@ function initialize!(integrator, cache::NystromCCDefaultInitialization)
     integrator.fsalfirst = ArrayPartition((kdu, ku))
 end
 
-const NystromDefaultInitialization = Union{Nystrom4Cache,
-    Nystrom4VelocityIndependentCache,
-    Nystrom5VelocityIndependentCache,
-    IRKN3Cache, IRKN4Cache,
-    DPRKN4Cache, DPRKN5Cache,
-    DPRKN6FMCache, DPRKN8Cache,
-    DPRKN12Cache, ERKN4Cache,
-    ERKN5Cache, ERKN7Cache}
+const NystromDefaultInitialization = Union{Nystrom4Cache, FineRKN5Cache,
+   Nystrom4VelocityIndependentCache,
+   Nystrom5VelocityIndependentCache,
+   IRKN3Cache, IRKN4Cache,
+   DPRKN4Cache, DPRKN5Cache,
+   DPRKN6FMCache, DPRKN8Cache,
+   DPRKN12Cache, ERKN4Cache,
+   ERKN5Cache, ERKN7Cache}
 
 function initialize!(integrator, cache::NystromDefaultInitialization)
     @unpack fsalfirst, k = cache
@@ -122,6 +122,113 @@ end
     integrator.stats.nf2 += 1
 end
 
+@muladd function perform_step!(integrator, cache::FineRKN5ConstantCache,
+                               repeat_step = false)
+    @unpack t, dt, f, p = integrator
+    duprev, uprev = integrator.uprev.x
+    @unpack c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a43, a51, a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar73, abar74, abar75, abar76, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, bbar6 = cache
+    k1 = integrator.fsalfirst.x[1]
+
+    ku = uprev + (c2 * dt * duprev) + (dt^2 * a21 * k1)
+    kdu = duprev + dt * (abar21 * k1)
+
+    k2 = f.f1(kdu, ku, p, t + dt * c2)
+    ku = uprev + dt * (c3 * duprev + dt * (a31 * k1 + a32 * k2))
+    kdu = duprev + dt * (abar31 * k1 + abar32 * k2)
+
+    k3 = f.f1(kdu, ku, p, t + dt * c3)
+    ku = uprev + dt * (c4 * duprev + dt * (a41 * k1 + a43 * k3)) # a42 = 0
+    kdu = duprev + dt * (abar41 * k1 + abar42 * k2 + abar43 * k3)
+
+    k4 = f.f1(kdu, ku, p, t + dt * c4)
+    ku = uprev + dt * (c5 * duprev + dt * (a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4))
+    kdu = duprev + dt * (abar51 * k1 + abar52 * k2 + abar53 * k3 + abar54 * k4)
+
+    k5 = f.f1(duprev, ku, p, t + dt * c5)
+    ku = uprev +
+         dt * (c6 * duprev + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4)) # a65 = 0
+    kdu = duprev +
+          dt * (abar61 * k1 + abar62 * k2 + abar63 * k3 + abar64 * k4 + abar65 * k5)
+
+    k6 = f.f1(duprev, ku, p, t + dt * c6)
+    ku = uprev +
+         dt * (c7 * duprev +
+               dt * (a71 * k1 + a73 * k3 + a74 * k4 + a75 * k5)) # a72 = a76 = 0
+    kdu = duprev +
+          dt * (abar71 * k1 + abar73 * k3 + abar74 * k4 + abar75 * k5 +
+                abar76 * k6) # abar72 = 0
+
+    k7 = f.f1(duprev, ku, p, t + dt * c7)
+    u = uprev + dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)) # no b6, b7 
+    du = duprev + dt * (bbar1 * k1 + bbar3 * k3 + bbar4 * k4 + bbar5 * k5 + bbar6 * k6) # no b2, b7
+
+    integrator.u = ArrayPartition((du, u))
+    integrator.fsallast = ArrayPartition((f.f1(du, u, p, t + dt), f.f2(du, u, p, t + dt)))
+    integrator.stats.nf += 4
+    integrator.stats.nf2 += 1
+    integrator.k[1] = integrator.fsalfirst
+    integrator.k[2] = integrator.fsallast
+end
+
+@muladd function perform_step!(integrator, cache::FineRKN5Cache, repeat_step = false)
+    @unpack t, dt, f, p = integrator
+    du, u = integrator.u.x
+    duprev, uprev = integrator.uprev.x
+    @unpack tmp, fsalfirst, k2, k3, k4, k5, k6, k7, k = cache
+    @unpack c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a43, a51, a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar73, abar74, abar75, abar76, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, bbar6 = cache.tab
+
+    kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
+    k1 = integrator.fsalfirst.x[1]
+
+    @.. broadcast=false ku=uprev + dt * (c2 * duprev + dt * (a21 * k1))
+    @.. broadcast=false kdu=duprev + dt * (abar21 * k1)
+
+    f.f1(k2, kdu, ku, p, t + dt * c2)
+    @.. broadcast=false ku=uprev + dt * (c3 * duprev + dt * (a31 * k1 + a32 * k2))
+    @.. broadcast=false kdu=duprev + dt * (abar31 * k1 + abar32 * k2)
+
+    f.f1(k3, kdu, ku, p, t + dt * c3)
+    @.. broadcast=false ku=uprev +
+                           dt * (c4 * duprev + dt * (a41 * k1 + a43 * k3)) # a42 = 0
+    @.. broadcast=false kdu=duprev + dt * (abar41 * k1 + abar42 * k2 + abar43 * k3)
+
+    f.f1(k4, kdu, ku, p, t + dt * c4)
+    @.. broadcast=false ku=uprev +
+                           dt *
+                           (c5 * duprev + dt * (a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4))
+    @.. broadcast=false kdu=duprev +
+                            dt * (abar51 * k1 + abar52 * k2 + abar53 * k3 + abar54 * k4)
+
+    f.f1(k5, kdu, ku, p, t + dt * c5)
+    @.. broadcast=false ku=uprev +
+                           dt * (c6 * duprev +
+                            dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4)) # a65 = 0
+    @.. broadcast=false kdu=duprev +
+                            dt * (abar61 * k1 + abar62 * k2 + abar63 * k3 + abar64 * k4 +
+                             abar65 * k5)
+
+    f.f1(k6, kdu, ku, p, t + dt * c6)
+    @.. broadcast=false ku=uprev +
+                           dt * (c7 * duprev +
+                            dt * (a71 * k1 + a73 * k3 + a74 * k4 + a75 * k5)) # a72 = a76 = 0
+
+    @.. broadcast=false kdu=duprev +
+                            dt * (abar71 * k1 + abar73 * k3 + abar74 * k4 +
+                             abar75 * k5 + abar76 * k6) # abar72 = 0
+
+    f.f1(k7, kdu, ku, p, t + dt * c7)
+    @.. broadcast=false u=uprev +
+                          dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5))
+    @.. broadcast=false du=duprev +
+                           dt *
+                           (bbar1 * k1 + bbar3 * k3 + bbar4 * k4 + bbar5 * k5 + bbar6 * k6)
+
+    f.f1(k.x[1], du, u, p, t + dt)
+    f.f2(k.x[2], du, u, p, t + dt)
+    integrator.stats.nf += 1
+    integrator.stats.nf2 += 1
+end
+
 @muladd function perform_step!(integrator, cache::Nystrom4VelocityIndependentConstantCache,
     repeat_step = false)
     @unpack t, dt, f, p = integrator