Merge branch 'master' into myb/bump

src/OrdinaryDiffEq.jl

@@ -149,7 +149,7 @@ module OrdinaryDiffEq
          SSPRK54, SSPRK104, RK4, ExplicitRK, OwrenZen3, OwrenZen4, OwrenZen5,
          BS3, BS5, CarpenterKennedy2N54,
          DP5, DP5Threaded, Tsit5, DP8, Vern6, Vern7, Vern8, TanYam7, TsitPap8,
-         Vern9,Feagin10, Feagin12, Feagin14, CompositeAlgorithm
+         Vern9,Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5
 
   export ImplicitEuler, ImplicitMidpoint, Trapezoid, TRBDF2, SDIRK2, Kvaerno3,
          KenCarp3, Cash4, Hairer4, Hairer42, SSPSDIRK2, Kvaerno4, Kvaerno5,

src/alg_utils.jl

@@ -118,6 +118,7 @@ alg_order(alg::EPIRK5P2) = 5
 alg_order(alg::EXPRB53s3) = 5
 alg_order(alg::SplitEuler) = 1
 alg_order(alg::ETD2) = 2
+alg_order(alg::Anas5) = 5
 
 alg_order(alg::SymplecticEuler) = 1
 alg_order(alg::VelocityVerlet) = 2

src/algorithms.jl

@@ -37,6 +37,11 @@ struct Heun <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct Ralston <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct Midpoint <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct RK4 <: OrdinaryDiffEqAdaptiveAlgorithm end
+struct Anas5{T} <: OrdinaryDiffEqAlgorithm
+  w::T		
+end
+Base.@pure Anas5(; w=1) = Anas5(w)
+
 struct OwrenZen3 <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct OwrenZen4 <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct OwrenZen5 <: OrdinaryDiffEqAdaptiveAlgorithm end

src/caches/low_order_rk_caches.jl

@@ -504,3 +504,39 @@ function alg_cache(alg::DP5Threaded,u,rate_prototype,uEltypeNoUnits,uBottomEltyp
   tab = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
   DP5ThreadedCache(u,uprev,k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp,tab)
 end
+
+struct Anas5Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
+  u::uType
+  uprev::uType
+  k1::rateType
+  k2::rateType
+  k3::rateType
+  k4::rateType
+  k5::rateType
+  k6::rateType
+  k7::rateType
+  utilde::uArrayType
+  tmp::uType
+  atmp::uEltypeNoUnits
+  tab::TabType
+end
+
+u_cache(c::Anas5Cache) = (c.atmp,c.utilde)
+du_cache(c::Anas5Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7)
+
+function alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
+  tab = Anas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
+  k1 = zeros(rate_prototype)
+  k2 = zeros(rate_prototype)
+  k3 = zeros(rate_prototype)
+  k4 = zeros(rate_prototype)
+  k5 = zeros(rate_prototype)
+  k6 = zeros(rate_prototype)
+  k7 = zeros(rate_prototype)
+  utilde = similar(u,indices(u))
+  atmp = similar(u,uEltypeNoUnits)
+  tmp = similar(u)
+  Anas5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp,tab)
+end
+
+alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Anas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

src/integrators/integrator_utils.jl

@@ -270,7 +270,7 @@ function loopfooter!(integrator)
       # integrator.EEst has unitless type of integrator.t
       if typeof(integrator.EEst)<: AbstractFloat && !isempty(integrator.opts.tstops)
         tstop = top(integrator.opts.tstops)
-        abs(ttmp - tstop) < 10eps(typeof(integrator.t))*oneunit(integrator.t) ?
+        abs(ttmp - tstop) < 10eps(integrator.t/oneunit(integrator.t))*oneunit(integrator.t) ?
                                   (integrator.t = tstop) : (integrator.t = ttmp)
       else
         integrator.t = ttmp
@@ -283,7 +283,7 @@ function loopfooter!(integrator)
     # integrator.EEst has unitless type of integrator.t
     if typeof(integrator.EEst)<: AbstractFloat && !isempty(integrator.opts.tstops)
       tstop = top(integrator.opts.tstops)
-      abs(ttmp - tstop) < 10eps(typeof(integrator.t))*oneunit(integrator.t) ?
+      abs(ttmp - tstop) < 10eps(integrator.t/oneunit(integrator.t))*oneunit(integrator.t) ?
                                   (integrator.t = tstop) : (integrator.t = ttmp)
     else
       integrator.t = ttmp

src/perform_step/fixed_timestep_perform_step.jl

@@ -392,3 +392,90 @@ end
 
   f( k,  u, p, t+dt)
 end
+
+function initialize!(integrator, cache::Anas5ConstantCache)
+  integrator.kshortsize = 7
+  integrator.k = typeof(integrator.k)(integrator.kshortsize)
+  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
+
+  # Avoid undefined entries if k is an array of arrays
+  integrator.fsallast = zero(integrator.fsalfirst)
+  integrator.k[1] = integrator.fsalfirst
+  @inbounds for i in 2:integrator.kshortsize-1
+    integrator.k[i] = zero(integrator.fsalfirst)
+  end
+  integrator.k[integrator.kshortsize] = integrator.fsallast
+end
+
+@muladd function perform_step!(integrator, cache::Anas5ConstantCache, repeat_step=false)
+  @unpack t,dt,uprev,u,f,p = integrator
+  @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,c2,c3,c4,c5,c6,b1,b3,b4,b5,b6 = cache
+  ## Note that c1 and b2 were 0.
+  w = integrator.alg.w
+  v = w*dt
+  ## Formula by Z.A. Anastassi, see the Anas5 caches in tableaus/low_order_rk_tableaus.jl for the full citation.
+  a65 = (-8000//1071)*(-a43*(v^5) + 6*tan(v)*(v^4) + 24*(v^3) - 72*tan(v)*(v^2) - 144*v + 144*tan(v))/((v^5)*(a43*tan(v)*v + 12 - 10*a43))
+  a61 += (-119//200)*a65
+  a63 += (189//100)*a65
+  a64 += (-459//200)*a65
+  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+2*k3), p, t+c4*dt)
+  k5 = f(uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4), p, t+c5*dt)
+  k6 = f(uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5), p, t+c6*dt)
+  u = uprev+dt*(b1*k1+b3*k3+b4*k4+b5*k5+b6*k6)
+  k7 = f(u, p, t+dt); integrator.fsallast = k7
+  integrator.k[1]=k1; integrator.k[2]=k2; integrator.k[3]=k3; integrator.k[4]=k4
+  integrator.k[5]=k5; integrator.k[6]=k6; integrator.k[7]=k7;
+  integrator.u = u
+end
+
+function initialize!(integrator, cache::Anas5Cache)
+  integrator.kshortsize = 7
+  resize!(integrator.k, integrator.kshortsize)
+  integrator.k[1]=cache.k1; integrator.k[2]=cache.k2;
+  integrator.k[3]=cache.k3; integrator.k[4]=cache.k4;
+  integrator.k[5]=cache.k5; integrator.k[6]=cache.k6;
+  integrator.k[7]=cache.k7;
+  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k7  # setup pointers
+  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
+end
+
+@muladd function perform_step!(integrator, cache::Anas5Cache, repeat_step=false)
+  @unpack t,dt,uprev,u,f,p = integrator
+  uidx = eachindex(integrator.uprev)
+  @unpack k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp = cache
+  @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,c2,c3,c4,c5,c6,b1,b3,b4,b5,b6 = cache.tab
+  w = integrator.alg.w
+  v = w*dt
+  ## Formula by Z.A. Anastassi, see the Anas5 caches in tableaus/low_order_rk_tableaus.jl for the full citation.
+  a65 = (-8000//1071)*(-a43*(v^5) + 6*tan(v)*(v^4) + 24*(v^3) - 72*tan(v)*(v^2) - 144*v + 144*tan(v))/((v^5)*(a43*tan(v)*v + 12 - 10*a43))
+  a61 += (-119//200)*a65
+  a63 += (189//100)*a65
+  a64 += (-459//200)*a65
+  @tight_loop_macros for i in uidx
+    @inbounds tmp[i] = uprev[i]+a21*k1[i]
+  end
+  f(k2, tmp, p, t+c2*dt)
+  @tight_loop_macros for i in uidx
+    @inbounds tmp[i] = uprev[i]+dt*(a31*k1[i]+a32*k2[i])
+  end
+  f(k3, tmp, p, t+c3*dt)
+  @tight_loop_macros for i in uidx
+    @inbounds tmp[i] = uprev[i]+dt*(a41*k1[i]+a42*k2[i]+2*k3[i])
+  end
+  f(k4, tmp, p, t+c4*dt)
+  @tight_loop_macros for i in uidx
+    @inbounds tmp[i] = uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
+  end
+  f(k5, tmp, p, t+c5*dt)
+  @tight_loop_macros for i in uidx
+    @inbounds tmp[i] = uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
+  end
+  f(k6, tmp, p, t+c6*dt)
+  @tight_loop_macros for i in uidx
+    @inbounds u[i] = uprev[i]+dt*(b1*k1[i]+b3*k3[i]+b4*k4[i]+b5*k5[i]+b6*k6[i])
+  end
+  f(k7, u, p, t+dt)
+end

src/solve.jl

@@ -19,9 +19,11 @@ function DiffEqBase.__init(
   d_discontinuities= eltype(prob.tspan)[],
   save_idxs = nothing,
   save_everystep = isempty(saveat),
-  save_timeseries = nothing,save_start = true,save_end = true,
+  save_timeseries = nothing,
+  save_start = save_everystep || isempty(saveat) || typeof(saveat) <: Number ? true : prob.tspan[1] in saveat,
+  save_end = save_everystep || isempty(saveat) || typeof(saveat) <: Number ? true : prob.tspan[2] in saveat,
   callback=nothing,
-  dense = save_everystep && !(typeof(alg) <: FunctionMap),
+  dense = save_everystep && !(typeof(alg) <: FunctionMap) && isempty(saveat),
   calck = (callback != nothing && callback != CallbackSet()) || # Empty callback
           (prob.callback != nothing && prob.callback != CallbackSet()) || # Empty prob.callback
           (!isempty(setdiff(saveat,tstops)) || dense), # and no dense output

src/tableaus/low_order_rk_tableaus.jl

@@ -1386,3 +1386,94 @@ function DP5_dense_bs(T)
   return b1,b3,b4,b5,b6,b7
 end
 =#
+
+"""
+An Optimized Runge-Kutta method for the solution of Orbital Problems
+by Z.A. Anastassi and T.E. Simos
+Journal of Computational and Applied Mathematics, Volume 175, Issue 1, 1 March 2005, Pages 1 to 9.
+"""
+struct Anas5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
+  a21::T
+  a31::T
+  a32::T
+  a41::T
+  a42::T
+  a43::T
+  a51::T
+  a52::T
+  a53::T
+  a54::T
+  a61::T
+  a62::T
+  a63::T
+  a64::T
+  a65::T
+  c2::T2
+  c3::T2
+  c4::T2
+  c5::T2
+  c6::T2
+  b1::T
+  b3::T
+  b4::T
+  b5::T
+  b6::T
+end
+
+Base.@pure function Anas5ConstantCache{T<:CompiledFloats,T2<:CompiledFloats}(::Type{T},::Type{T2})
+    a21 = T(0.1)
+    a31 = T(-0.2222222222222222)
+    a32 = T(0.5555555555555556)
+    a41 = T(3.111111111111111)
+    a42 = T(-4.444444444444444)
+    a43 = T(2.0)
+    a51 = T(-1.409625)
+    a52 = T(2.1375)
+    a53 = T(-0.2295)
+    a54 = T(0.401625)
+    a61 = T(-4.0)
+    a62 = T(5.0)
+    a63 = T(0.0)
+    a64 = T(0.0)
+    a65 = T(0.0)
+    c2 = T2(0.1)
+    c3 = T2(0.3333333333333333)
+    c4 = T2(0.6666666666666667)
+    c5 = T2(0.9)
+    c6 = T2(1)
+    b1 = T(0.1064814814814815)
+    b3 = T(0.4632352941176471)
+    b4 = T(0.1607142857142857)
+    b5 = T(0.3112356053532524)
+    b6 = T(-.04166666666666667)
+    Anas5ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,c2,c3,c4,c5,c6,b1,b3,b4,b5,b6)
+end
+
+Base.@pure function Anas5ConstantCache{T,T2}(::Type{T},::Type{T2})
+    a21 = T(1//10)
+    a31 = T(-2//9)
+    a32 = T(5//9)
+    a41 = T(28//9)
+    a42 = T(-40//9)
+    a43 = T(2//1)
+    a51 = T(-11277//8000)
+    a52 = T(171//80)
+    a53 = T(-459//2000)
+    a54 = T(3213//8000)
+    a61 = T(-4//1)
+    a62 = T(5//1)
+    a63 = T(0//1)
+    a64 = T(0//1)
+    a65 = T(0//1)
+    c2 = T2(1//10)
+    c3 = T2(1//3)
+    c4 = T2(2//3)
+    c5 = T2(9//10)
+    c6 = T2(1//1)
+    b1 = T(23//216)
+    b3 = T(63//136)
+    b4 = T(9//56)
+    b5 = T(1000//3213)
+    b6 = T(-1//24)
+    Anas5ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,c2,c3,c4,c5,c6,b1,b3,b4,b5,b6)
+end

test/ode/ode_convergence_tests.jl

@@ -51,6 +51,8 @@ for i = 1:2
   @test abs(sim105.𝒪est[:l2]-4) < testTol
   sim106 = test_convergence(dts,prob,VCABM5())
   @test abs(sim106.𝒪est[:l2]-5) < testTol
+  sim160 = test_convergence(dts,prob,Anas5(w=2))
+  @test abs(sim160.𝒪est[:l2]-4) < 2*testTol
 
   println("Stiff Solvers")
 

test/ode/ode_saveat_tests.jl

@@ -4,65 +4,67 @@ import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinea
 @testset "saveat Tests" begin
   prob = prob_ode_linear
 
-  sol =solve(prob,DP5(),dt=1//2^(2),save_everystep=false,dense=false)
-  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,dense=false,saveat=[1/2])
+  sol =solve(prob,DP5(),dt=1//2^(2),save_everystep=false)
+  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,saveat=[1/2])
 
-  @test symdiff(sol.t,sol2.t) == [1/2]
+  @test sort!(symdiff(sol.t,sol2.t)) == [0.0,1/2,1.0]
 
-  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,dense=false,saveat=1/2)
+  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,saveat=[0.0,1/2,1.0])
 
   @test symdiff(sol.t,sol2.t) == [1/2]
 
-  sol3=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,dense=false,saveat=[1/2],tstops=[1/2])
+  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,saveat=1/2)
+
+  @test sort!(symdiff(sol.t,sol2.t)) == [1/2]
+
+  sol3=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,saveat=[1/2],tstops=[1/2])
 
-  @test sol3.t == [0.0,0.5,1.00]
+  @test sol3.t == [1/2]
 
-  sol3=solve(prob,DP5(),dt=1//2^(2),saveat=[1/2],tstops=[1/2])
+  sol3=solve(prob,DP5(),dt=1//2^(2),saveat=[0.0,1/2,1.0],tstops=[1/2])
 
-  @test sol3.t == [0.0,0.5,1.00]
+  @test sol3.t == [0.0,1/2,1.0]
 
   sol3=solve(prob,DP5(),dt=1//2^(2),saveat=1/10,tstops=[1/2])
 
   @test sol3.t == collect(0.0:0.1:1.00)
 
-  #plot(0:1//2^(4):1,interpd)
-
-  sol =solve(prob,RK4(),dt=1/2^(2),save_everystep=true,dense=false)
-  sol2=solve(prob,RK4(),dt=1/2^(2),save_everystep=true,dense=false,saveat=[.125,.6,.61,.8])
+  sol =solve(prob,RK4(),dt=1/2^(2),save_everystep=true,adaptive=false)
+  sol2=solve(prob,RK4(),dt=1/2^(2),save_everystep=true,adaptive=false,saveat=[.125,.6,.61,.8])
 
   @test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]
 
-  sol =solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true,dense=false)
-  sol2=solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true,dense=false,saveat=[.125,.6,.61,.8])
+  sol =solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true)
+  sol2=solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true,saveat=[.125,.6,.61,.8])
 
   @test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]
 
-  sol =solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true,dense=false)
-  sol2=solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true,dense=false,saveat=[.125,.6,.61,.8])
+  sol =solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true)
+  sol2=solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true,saveat=[.125,.6,.61,.8])
 
   @test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]
 
   prob = prob_ode_2Dlinear
 
-  sol =solve(prob,DP5(),dt=1//2^(2),save_everystep=true,dense=false)
-  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=true,dense=false,saveat=[1/2])
+  sol =solve(prob,DP5(),dt=1//2^(2),save_everystep=true)
+  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=true,saveat=[0.0,1/2,1.0])
 
   @test symdiff(sol.t,sol2.t) == [1/2]
 
-  sol =solve(prob,RK4(),dt=1/2^(2),save_everystep=true,dense=false)
-  sol2=solve(prob,RK4(),dt=1/2^(2),save_everystep=true,dense=false,saveat=[.125,.6,.61,.8])
+  sol =solve(prob,RK4(),dt=1/2^(2),save_everystep=true,adaptive=false)
+  sol2=solve(prob,RK4(),dt=1/2^(2),save_everystep=true,adaptive=false,saveat=[0.0,.125,.6,.61,.8,1.0])
 
   @test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]
 
-  sol =solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true,dense=false)
-  sol2=solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true,dense=false,saveat=[.125,.6,.61,.8])
+  sol =solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true)
+  sol2=solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true,saveat=[.125,.6,.61,.8])
 
   @test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]
 
-  sol =solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true,dense=false)
-  sol2=solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true,dense=false,saveat=[.125,.6,.61,.8])
+  sol =solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=false)
+  sol2=solve(prob,GenericTrapezoid(),dt=1/2^(2),saveat=[.125,.6,.61,.8])
 
-  @test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]
+  @test sort!(symdiff(sol.t,sol2.t)) == [0.0,.125,.6,.61,.8,1.0]
 
   sol=solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true,dense=false,saveat=[0,.125,.6,.61,.8])