New algorithm FineRKN4()

docs/src/dynamical/nystrom.md

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

src/OrdinaryDiffEq.jl

@@ -437,7 +437,7 @@ export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog,
 
 export SplitEuler
 
-export Nystrom4, FineRKN5, Nystrom4VelocityIndependent,
+export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent,
     Nystrom5VelocityIndependent,
     IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7
 

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::FineRKN4) = 4
 alg_order(alg::FineRKN5) = 5
 alg_order(alg::Nystrom4VelocityIndependent) = 4
 alg_order(alg::IRKN4) = 4

src/algorithms.jl

@@ -5251,6 +5251,28 @@ E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equatio
 """
 struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end
 
+"""
+    FineRKN4()
+
+A 4th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs.
+In particular, this method allows the acceleration equation to depend on the velocity.
+
+## 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 FineRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end
+
 """
     FineRKN5()
 

src/caches/rkn_caches.jl

@@ -36,6 +36,48 @@ 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 FineRKN4Cache{uType, rateType, reducedRateType, uNoUnitsType, TabType} <:
+              OrdinaryDiffEqMutableCache
+    u::uType
+    uprev::uType
+    fsalfirst::rateType
+    k2::reducedRateType
+    k3::reducedRateType
+    k4::reducedRateType
+    k5::reducedRateType
+    k::rateType
+    utilde::uType
+    tmp::uType
+    atmp::uNoUnitsType
+    tab::TabType
+end
+
+function alg_cache(alg::FineRKN4, 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 = FineRKN4ConstantCache(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)
+    k = zero(rate_prototype)
+    utilde = zero(u)
+    atmp = similar(u, uEltypeNoUnits)
+    recursivefill!(atmp, false)
+    tmp = zero(u)
+    FineRKN4Cache(u, uprev, k1, k2, k3, k4, k5, k, utilde, tmp, atmp, tab)
+end
+
+function alg_cache(alg::FineRKN4, u, rate_prototype, ::Type{uEltypeNoUnits},
+    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+    dt, reltol, p, calck,
+    ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    FineRKN4ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+end
+
 @cache struct FineRKN5Cache{uType, rateType, reducedRateType, uNoUnitsType, TabType} <:
               OrdinaryDiffEqMutableCache
     u::uType

src/perform_step/rkn_perform_step.jl

@@ -4,7 +4,8 @@
 ## y₁ = y₀ + hy'₀ + h²∑b̄ᵢk'ᵢ
 ## y'₁ = y'₀ + h∑bᵢk'ᵢ
 
-const NystromCCDefaultInitialization = Union{Nystrom4ConstantCache, FineRKN5ConstantCache,
+const NystromCCDefaultInitialization = Union{Nystrom4ConstantCache, FineRKN4ConstantCache,
+    FineRKN5ConstantCache,
     Nystrom4VelocityIndependentConstantCache,
     Nystrom5VelocityIndependentConstantCache,
     IRKN3ConstantCache, IRKN4ConstantCache,
@@ -25,7 +26,7 @@ function initialize!(integrator, cache::NystromCCDefaultInitialization)
     integrator.fsalfirst = ArrayPartition((kdu, ku))
 end
 
-const NystromDefaultInitialization = Union{Nystrom4Cache, FineRKN5Cache,
+const NystromDefaultInitialization = Union{Nystrom4Cache, FineRKN4Cache, FineRKN5Cache,
     Nystrom4VelocityIndependentCache,
     Nystrom5VelocityIndependentCache,
     IRKN3Cache, IRKN4Cache,
@@ -122,6 +123,115 @@ end
     integrator.stats.nf2 += 1
 end
 
+@muladd function perform_step!(integrator, cache::FineRKN4ConstantCache,
+    repeat_step = false)
+    @unpack t, dt, f, p = integrator
+    duprev, uprev = integrator.uprev.x
+    @unpack c2, c3, c4, c5, a21, a31, a32, a41, a43, a51,
+    a52, a53, a54, abar21, abar31, abar32, abar41, abar42, abar43, abar51,
+    abar52, abar53, abar54, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, btilde1, btilde3, btilde4, btilde5, bptilde1,
+    bptilde3, bptilde4, bptilde5 = cache
+    k1 = integrator.fsalfirst.x[1]
+
+    ku = uprev + dt * (c2 * duprev + dt * (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(kdu, ku, p, t + dt * c5)
+
+    u = uprev + dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)) # b2 = 0
+    du = duprev + dt * (bbar1 * k1 + bbar3 * k3 + bbar4 * k4 + bbar5 * k5) # bbar2 = 0
+
+    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 += 5
+    integrator.stats.nf2 += 1
+    integrator.k[1] = integrator.fsalfirst
+    integrator.k[2] = integrator.fsallast
+
+    if integrator.opts.adaptive
+        dtsq = dt^2
+        uhat = dtsq * (btilde1 * k1 + btilde3 * k3 + btilde4 * k4 + btilde5 * k5) # btilde2 = 0
+        duhat = dt * (bptilde1 * k1 + bptilde3 * k3 + bptilde4 * k4 + bptilde5 * k5) # bptilde2 = 0
+        utilde = ArrayPartition((duhat, uhat))
+        atmp = calculate_residuals(utilde, integrator.uprev, integrator.u,
+            integrator.opts.abstol, integrator.opts.reltol,
+            integrator.opts.internalnorm, t)
+        integrator.EEst = integrator.opts.internalnorm(atmp, t)
+    end
+end
+
+@muladd function perform_step!(integrator, cache::FineRKN4Cache, repeat_step = false)
+    @unpack t, dt, f, p = integrator
+    du, u = integrator.u.x
+    duprev, uprev = integrator.uprev.x
+    @unpack tmp, atmp, fsalfirst, k2, k3, k4, k5, k, utilde = cache
+    @unpack c2, c3, c4, c5, a21, a31, a32, a41, a43, a51,
+    a52, a53, a54, abar21, abar31, abar32, abar41, abar42, abar43, abar51,
+    abar52, abar53, abar54, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, btilde1, btilde3, btilde4, btilde5, bptilde1,
+    bptilde3, bptilde4, bptilde5 = cache.tab
+    kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
+    uidx = eachindex(integrator.uprev.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 u=uprev +
+                          dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)) # b2 = 0
+    @.. broadcast=false du=duprev +
+                           dt *
+                           (bbar1 * k1 + bbar3 * k3 + bbar4 * k4 + bbar5 * k5) # bbar2 = 0
+
+    f.f1(k.x[1], du, u, p, t + dt)
+    f.f2(k.x[2], du, u, p, t + dt)
+    integrator.stats.nf += 5
+    integrator.stats.nf2 += 1
+    if integrator.opts.adaptive
+        duhat, uhat = utilde.x
+        dtsq = dt^2
+        @.. broadcast=false uhat=dtsq *
+                                 (btilde1 * k1 + btilde3 * k3 + btilde4 * k4 +
+                                  btilde5 * k5) # btilde2 = 0
+        @.. broadcast=false duhat=dt *
+                                  (bptilde1 * k1 + bptilde3 * k3 + bptilde4 * k4 +
+                                   bptilde5 * k5) # bptilde2 = 0
+
+        calculate_residuals!(atmp, utilde, integrator.uprev, integrator.u,
+            integrator.opts.abstol, integrator.opts.reltol,
+            integrator.opts.internalnorm, t)
+        integrator.EEst = integrator.opts.internalnorm(atmp, t)
+    end
+end
+
 @muladd function perform_step!(integrator, cache::FineRKN5ConstantCache,
     repeat_step = false)
     @unpack t, dt, f, p = integrator

src/tableaus/rkn_tableaus.jl

@@ -1,3 +1,104 @@
+struct FineRKN4ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache
+    c1::T2
+    c2::T2
+    c3::T2
+    c4::T2
+    c5::T2
+    a21::T
+    a31::T
+    a32::T
+    a41::T
+    #a42::T
+    a43::T
+    a51::T
+    a52::T
+    a53::T
+    a54::T
+    abar21::T
+    abar31::T
+    abar32::T
+    abar41::T
+    abar42::T
+    abar43::T
+    abar51::T
+    abar52::T
+    abar53::T
+    abar54::T
+    b1::T
+    #b2::T
+    b3::T
+    b4::T
+    b5::T
+    bbar1::T
+    #bbar2::T
+    bbar3::T
+    bbar4::T
+    bbar5::T
+    btilde1::T
+    #btilde2::T
+    btilde3::T
+    btilde4::T
+    btilde5::T
+    bptilde1::T
+    #bptilde2::T
+    bptilde3::T
+    bptilde4::T
+    bptilde5::T
+end
+
+function FineRKN4ConstantCache(T::Type, T2::Type)
+    c1 = convert(T2, 1 // 1)
+    c2 = convert(T2, 2 // 9)
+    c3 = convert(T2, 1 // 3)
+    c4 = convert(T2, 3 // 4)
+    c5 = convert(T2, 1 // 1)
+    a21 = convert(T, 2 // 81)
+    a31 = convert(T, 1 // 36)
+    a32 = convert(T, 1 // 36)
+    a41 = convert(T, 9 // 128)
+    #a42 = convert(T, 0 // 1)
+    a43 = convert(T, 27 // 128)
+    a51 = convert(T, 11 // 60)
+    a52 = convert(T, -3 // 20)
+    a53 = convert(T, 9 // 25)
+    a54 = convert(T, 8 // 75)
+    abar21 = convert(T, 2 // 9)
+    abar31 = convert(T, 1 // 12)
+    abar32 = convert(T, 1 // 4)
+    abar41 = convert(T, 69 // 128)
+    abar42 = convert(T, -243 // 128)
+    abar43 = convert(T, 135 // 64)
+    abar51 = convert(T, -17 // 12)
+    abar52 = convert(T, 27 // 4)
+    abar53 = convert(T, -27 // 5)
+    abar54 = convert(T, 16 // 15)
+    b1 = convert(T, 19 // 180)
+    #b2 = convert(T, 0 // 1)
+    b3 = convert(T, 63 // 200)
+    b4 = convert(T, 16 // 225)
+    b5 = convert(T, 1 // 120)
+    bbar1 = convert(T, 1 // 9)
+    #bbar2 = convert(T, 0 // 1)
+    bbar3 = convert(T, 9 // 20)
+    bbar4 = convert(T, 16 // 45)
+    bbar5 = convert(T, 1 // 12)
+    btilde1 = convert(T, 25 // 1116)
+    #btilde2 = convert(T, 0 // 1) 
+    btilde3 = convert(T, -63 // 1240)
+    btilde4 = convert(T, 64 // 1395)
+    btilde5 = convert(T, -13 // 744)
+    bptilde1 = convert(T, 2 // 125)
+    #bptilde2 = convert(T, 0 // 1) 
+    bptilde3 = convert(T, -27 // 625)
+    bptilde4 = convert(T, 32 // 625)
+    bptilde5 = convert(T, -3 // 125)
+    FineRKN4ConstantCache(c1, c2, c3, c4, c5, a21, a31, a32, a41, a43, a51,
+        a52, a53, a54, abar21, abar31, abar32, abar41, abar42, abar43, abar51,
+        abar52, abar53, abar54, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, btilde1,
+        btilde3, btilde4, btilde5, bptilde1,
+        bptilde3, bptilde4, bptilde5)
+end
+
 struct FineRKN5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache
     c1::T2
     c2::T2