Merge remote-tracking branch 'origin/master'

Author: ChrisRackauckas

src/OrdinaryDiffEq.jl

@@ -152,5 +152,5 @@ module OrdinaryDiffEq
   export SplitEuler
 
   export Nystrom4, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent,
-         IRKN3, IRKN4, DPRKN6, DPRKN8, DPRKN12, ERKN4
+         IRKN3, IRKN4, DPRKN6, DPRKN8, DPRKN12, ERKN4, ERKN5
 end # module

src/alg_utils.jl

@@ -102,6 +102,7 @@ isfsal(alg::DPRKN6) = true
 isfsal(alg::DPRKN8) = true
 isfsal(alg::DPRKN12) = true
 isfsal(alg::ERKN4) = true
+isfsal(alg::ERKN5) = true
 
 fsal_typeof(alg::OrdinaryDiffEqAlgorithm,rate_prototype) = typeof(rate_prototype)
 #fsal_typeof(alg::LawsonEuler,rate_prototype) = Vector{typeof(rate_prototype)}
@@ -269,6 +270,7 @@ alg_order(alg::DPRKN6) = 6
 alg_order(alg::DPRKN8) = 8
 alg_order(alg::DPRKN12) = 12
 alg_order(alg::ERKN4) = 4
+alg_order(alg::ERKN5) = 5
 
 alg_order(alg::Midpoint) = 2
 alg_order(alg::GenericIIF1) = 1

src/algorithms.jl

@@ -133,6 +133,7 @@ struct DPRKN6 <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct DPRKN8 <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct DPRKN12 <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct ERKN4 <: OrdinaryDiffEqAdaptiveAlgorithm end
+struct ERKN5 <: OrdinaryDiffEqAdaptiveAlgorithm end
 
 ################################################################################
 

src/caches/rkn_caches.jl

@@ -294,3 +294,34 @@ function alg_cache(alg::ERKN4,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev
   tmp = similar(u)
   ERKN4Cache(u,uprev,k1,k2,k3,k4,k,utilde,tmp,atmp,tab)
 end
+
+struct ERKN5Cache{uType,uArrayType,rateType,reducedRateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
+  u::uType
+  uprev::uType
+  fsalfirst::rateType
+  k2::reducedRateType
+  k3::reducedRateType
+  k4::reducedRateType
+  k::rateType
+  utilde::uArrayType
+  tmp::uType
+  atmp::uEltypeNoUnits
+  tab::TabType
+end
+
+u_cache(c::ERKN5Cache) = (c.atmp,c.utilde)
+du_cache(c::ERKN5Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.k)
+
+function alg_cache(alg::ERKN5,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,::Type{Val{true}})
+  reduced_rate_prototype = rate_prototype.x[2]
+  tab = ERKN5ConstantCache(real(uEltypeNoUnits),real(tTypeNoUnits))
+  k1 = zeros(rate_prototype)
+  k2 = zeros(reduced_rate_prototype)
+  k3 = zeros(reduced_rate_prototype)
+  k4 = zeros(reduced_rate_prototype)
+  k  = zeros(rate_prototype)
+  utilde = similar(u,indices(u))
+  atmp = similar(u,uEltypeNoUnits)
+  tmp = similar(u)
+  ERKN5Cache(u,uprev,k1,k2,k3,k4,k,utilde,tmp,atmp,tab)
+end

src/dense/interpolants.jl

@@ -576,7 +576,13 @@ end
   b12Θ = @evalpoly(Θ, 0,    0, r122, r123, r124, r125, r126)
 
   if out == nothing
-    return y₀[idxs] + dt*(k[1][idxs]*b1Θ + k[4][idxs]*b4Θ + k[5][idxs]*b5Θ + k[6][idxs]*b6Θ + k[7][idxs]*b7Θ + k[8][idxs]*b8Θ + k[9][idxs]*b9Θ + k[10][idxs]*b10Θ + k[11][idxs]*b11Θ + k[12][idxs]*b12Θ)
+    if idxs == nothing
+      # return @. y₀ + dt*(k[1]*b1Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ)
+      return y₀ + dt*(k[1]*b1Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ)
+    else
+      # return @. y₀[idxs] + dt*(k[1][idxs]*b1Θ + k[4][idxs]*b4Θ + k[5][idxs]*b5Θ + k[6][idxs]*b6Θ + k[7][idxs]*b7Θ + k[8][idxs]*b8Θ + k[9][idxs]*b9Θ + k[10][idxs]*b10Θ + k[11][idxs]*b11Θ + k[12][idxs]*b12Θ)
+      return y₀[idxs] + dt*(k[1][idxs]*b1Θ + k[4][idxs]*b4Θ + k[5][idxs]*b5Θ + k[6][idxs]*b6Θ + k[7][idxs]*b7Θ + k[8][idxs]*b8Θ + k[9][idxs]*b9Θ + k[10][idxs]*b10Θ + k[11][idxs]*b11Θ + k[12][idxs]*b12Θ)
+    end
   elseif idxs == nothing
     #@. out = y₀ + dt*(k[1]*b1Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ)
     @inbounds for i in eachindex(out)

src/integrators/rkn_integrators.jl

@@ -7,9 +7,9 @@
 const NystromDefaultInitialization = Union{Nystrom4Cache,
                                            Nystrom4VelocityIndependentCache,
                                            Nystrom5VelocityIndependentCache,
-                                           IRKN3Cache, IRKN4Cache,
+                                           IRKN3Cache,  IRKN4Cache,
                                            DPRKN8Cache, DPRKN12Cache,
-                                           ERKN4Cache}
+                                           ERKN4Cache,  ERKN5Cache}
 
 function initialize!(integrator,cache::NystromDefaultInitialization)
   @unpack fsalfirst,k = cache
@@ -492,3 +492,40 @@ end
     integrator.EEst = integrator.opts.internalnorm(atmp)
   end
 end
+
+@muladd function perform_step!(integrator,cache::ERKN5Cache,repeat_step=false)
+  @unpack t,dt,f = integrator
+  u,du = integrator.u.x
+  uprev,duprev = integrator.uprev.x
+  @unpack tmp,atmp,fsalfirst,k2,k3,k4,k,utilde = cache
+  @unpack c1, c2, c3, a21, a31, a32, a41, a42, a43, b1, b2, b3, b4, bp1, bp2, bp3, bp4, btilde1, btilde2, btilde3, btilde4 = cache.tab
+  ku, kdu = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
+  uidx = eachindex(integrator.uprev.x[2])
+  k1 = fsalfirst.x[2]
+
+  @. ku = uprev + dt*(c1*duprev + dt*a21*k1)
+
+  f.f2(t+dt*c1,ku,du,k2)
+  @. ku = uprev + dt*(c2*duprev + dt*(a31*k1 + a32*k2))
+
+  f.f2(t+dt*c2,ku,du,k3)
+  @. ku = uprev + dt*(c3*duprev + dt*(a41*k1 + a42*k2 + a43*k3))
+
+  f.f2(t+dt*c3,ku,du,k4)
+  @tight_loop_macros for i in uidx
+    @inbounds u[i]  = uprev[i] + dt*(duprev[i] + dt*(b1 *k1[i] + b2 *k2[i] + b3 *k3[i] + b4 *k4[i]))
+    @inbounds du[i] = duprev[i]                + dt*(bp1*k1[i] + bp2*k2[i] + bp3*k3[i] + bp4*k4[i])
+  end
+
+  f.f1(t+dt,u,du,k.x[1])
+  f.f2(t+dt,u,du,k.x[2])
+  if integrator.opts.adaptive
+    uhat, duhat = utilde.x
+    dtsq = dt^2
+    @tight_loop_macros for i in uidx
+      @inbounds uhat[i]  = dtsq*(btilde1*k1[i] + btilde2*k2[i] + btilde3*k3[i] + btilde4*k4[i])
+    end
+    calculate_residuals!(atmp, uhat, integrator.uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol)
+    integrator.EEst = integrator.opts.internalnorm(atmp)
+  end
+end

src/tableaus/rkn_tableaus.jl

@@ -227,6 +227,88 @@ ERKN4ConstantCache(
                   T(0.0016835016835016834))
 end
 
+struct ERKN5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
+  c1::T2
+  c2::T2
+  c3::T2
+  a21::T
+  a31::T
+  a32::T
+  a41::T
+  a42::T
+  a43::T
+  b1::T
+  b2::T
+  b3::T
+  b4::T
+  bp1::T # bp denotes bprime
+  bp2::T
+  bp3::T
+  bp4::T
+  btilde1::T
+  btilde2::T
+  btilde3::T
+  btilde4::T
+  # bptilde1::T
+  # bptilde2::T
+  # bptilde3::T
+  # bptilde4::T
+end
+
+function ERKN5ConstantCache(T::Type,T2::Type)
+  c1 = T2(1//2)
+  c2 = T2(19//70)
+  c3 = T2(44//51)
+  a21 = T(1//8)
+  a31 = T(2907//343000)
+  a32 = T(1216//42875)
+  a41 = T(6624772//Int64(128538819))
+  a42 = T(6273905//Int64(54121608))
+  a43 = T(Int64(210498365)//Int64(1028310552))
+  b1 = T(479//5016)
+  b2 = T(235//1776)
+  b3 = T(145775//641744)
+  b4 = T(309519//6873416)
+  btilde1 = T(479//5016 - 184883//2021250)
+  btilde2 = T(235//1776 - 411163//3399375)
+  btilde3 = T(145775//641744 - 6//25)
+  btilde4 = T(309519//6873416 - 593028//Int64(12464375))
+  bp1 = b1
+  bp2 = T(235//888)
+  bp3 = T(300125//962616)
+  bp4 = T(2255067//6873416)
+  # bptilde1 = T(0)
+  # bptilde2 = T(0)
+  # bptilde3 = T(0)
+  # bptilde4 = T(0)
+  ERKN5ConstantCache(c1, c2, c3, a21, a31, a32, a41, a42, a43, b1, b2, b3, b4, bp1, bp2, bp3, bp4, btilde1, btilde2, btilde3, btilde4)
+end
+
+Base.@pure function ERKN5ConstantCache{T<:CompiledFloats,T2<:CompiledFloats}(::Type{T},::Type{T2})
+  ERKN5ConstantCache(
+                     T2(0.5),
+                     T2(0.2714285714285714),
+                     T2(0.8627450980392157),
+                     T(0.125),
+                     T(0.008475218658892128),
+                     T(0.028361516034985424),
+                     T(0.051539076300366506),
+                     T(0.11592236875149756),
+                     T(0.20470310704348388),
+                     T(0.09549441786283891),
+                     T(0.13231981981981983),
+                     T(0.22715444164651324),
+                     T(0.04503132067082801),
+                     T(0.09549441786283891),
+                     T(0.26463963963963966),
+                     T(0.3117806061814888),
+                     T(0.32808533631603265),
+                     T(0.004024782736060931),
+                     T(0.011367291781577495),
+                     T(-0.012845558353486749),
+                     T(-0.0025465161641516788))
+end
+
 struct DPRKN6ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
   c1::T2
   c2::T2

test/partitioned_methods_tests.jl

@@ -144,6 +144,9 @@ sim = test_convergence(dts,prob_big,DPRKN12(),dense_errors=true)
 sim = test_convergence(dts,prob_big,ERKN4(),dense_errors=true)
 @test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1
 @test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
+sim = test_convergence(dts,prob_big,ERKN5(),dense_errors=true)
+@test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1
+@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
 
 # Adaptive methods regression test
 sol = solve(prob, DPRKN6())
@@ -152,8 +155,10 @@ sol = solve(prob, DPRKN8())
 @test length(sol.u) < 13
 sol = solve(prob, DPRKN12())
 @test length(sol.u) < 9
-sol = solve(prob, ERKN4())
-@test length(sol.u) < 15
+sol = solve(prob, ERKN4(),reltol=1e-8)
+@test length(sol.u) < 38
+sol = solve(prob, ERKN5(),reltol=1e-8)
+@test length(sol.u) < 29
 
 f = function (t,u,du)
   du.x[1] .= u.x[2]