From 85c4bd7b0a081e5aa36517b118820443f53fe089 Mon Sep 17 00:00:00 2001 From: Lucas Booher Date: Tue, 26 Jun 2018 02:20:27 -0700 Subject: [PATCH 1/9] Implementing the Anas5 method --- src/OrdinaryDiffEq.jl | 2 +- src/alg_utils.jl | 1 + src/algorithms.jl | 1 + src/caches/low_order_rk_caches.jl | 36 +++++++ .../fixed_timestep_perform_step.jl | 72 ++++++++++++++ src/tableaus/low_order_rk_tableaus.jl | 97 +++++++++++++++++++ 6 files changed, 208 insertions(+), 1 deletion(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index c8e71e99cf..148db9e325 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -146,7 +146,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, diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 11d6a12155..6fae0dc35e 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -110,6 +110,7 @@ alg_order(alg::ETDRK4) = 4 alg_order(alg::HochOst4) = 4 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 diff --git a/src/algorithms.jl b/src/algorithms.jl index 8b3ee94bcb..891ccaf302 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -37,6 +37,7 @@ struct Heun <: OrdinaryDiffEqAdaptiveAlgorithm end struct Ralston <: OrdinaryDiffEqAdaptiveAlgorithm end struct Midpoint <: OrdinaryDiffEqAdaptiveAlgorithm end struct RK4 <: OrdinaryDiffEqAdaptiveAlgorithm end +struct Anas5 <: OrdinaryDiffEqAlgorithm end struct OwrenZen3 <: OrdinaryDiffEqAdaptiveAlgorithm end struct OwrenZen4 <: OrdinaryDiffEqAdaptiveAlgorithm end struct OwrenZen5 <: OrdinaryDiffEqAdaptiveAlgorithm end diff --git a/src/caches/low_order_rk_caches.jl b/src/caches/low_order_rk_caches.jl index 2f814dd6cf..f45b376132 100644 --- a/src/caches/low_order_rk_caches.jl +++ b/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)) diff --git a/src/perform_step/fixed_timestep_perform_step.jl b/src/perform_step/fixed_timestep_perform_step.jl index 6155850070..5874f304ec 100644 --- a/src/perform_step/fixed_timestep_perform_step.jl +++ b/src/perform_step/fixed_timestep_perform_step.jl @@ -392,3 +392,75 @@ end f( k, u, p, t+dt) end + +function initialize!(integrator, cache::Anas5ConstantCache) + integrator.kshortsize = 6 + 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,a51,a52,a53,a54,a61,a62,a63,a64,a65,c1,c2,c3,c4,c5,c6,b1,b2,b3,b4,b5,b6 = cache + 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+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+b2*k2+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 = 6 + 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,a65,c1,c2,c3,c4,c5,c6,b1,b2,b3,b4,b5,b6,b7 = cache.tab + @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]+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]+b2*k2[i]+b3*k3[i]+b4*k4[i]+b5*k5[i]+b6*k6[i]) + end + f(k7, u, p, t+dt) +end diff --git a/src/tableaus/low_order_rk_tableaus.jl b/src/tableaus/low_order_rk_tableaus.jl index aad5fa4785..191ac55daa 100644 --- a/src/tableaus/low_order_rk_tableaus.jl +++ b/src/tableaus/low_order_rk_tableaus.jl @@ -1386,3 +1386,100 @@ 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 + c1::T2 + c2::T2 + c3::T2 + c4::T2 + c5::T2 + c6::T2 + b1::T + b2::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(-5.19) + a62 = T(5.0) + a63 = T(1.89) + a64 = T(-2.295) + a65 = T(1) + c1 = T2(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) + b2 = T(0.0) + 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,c1,c2,c3,c4,c5,c6,b1,b2,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) + a51 = T(-11277/8000) + a52 = T(171/80) + a53 = T(-459/2000) + a54 = T(3213/8000) + a61 = T(-519/100) + a62 = T(5) + a63 = T(189/100) + a64 = T(-459/200) + a65 = T(1) + c1 = T2(0) + c2 = T2(1/10) + c3 = T2(1/3) + c4 = T2(2/3) + c5 = T2(9/10) + c6 = T2(1) + b1 = T(23/216) + b2 = T(0) + 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,c1,c2,c3,c4,c5,c6,b1,b2,b3,b4,b5,b6) +end From c75df142599bdeae0d1cbf81239b0442cc4e5c55 Mon Sep 17 00:00:00 2001 From: Lucas Booher Date: Thu, 28 Jun 2018 20:00:17 -0700 Subject: [PATCH 2/9] Continuation on the implementation of Anas5 --- src/algorithms.jl | 4 +- .../fixed_timestep_perform_step.jl | 9 ++- src/tableaus/low_order_rk_tableaus.jl | 60 +++++++++---------- test/ode/ode_convergence_tests.jl | 2 + 4 files changed, 39 insertions(+), 36 deletions(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 891ccaf302..fb0b326e4a 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -37,7 +37,9 @@ struct Heun <: OrdinaryDiffEqAdaptiveAlgorithm end struct Ralston <: OrdinaryDiffEqAdaptiveAlgorithm end struct Midpoint <: OrdinaryDiffEqAdaptiveAlgorithm end struct RK4 <: OrdinaryDiffEqAdaptiveAlgorithm end -struct Anas5 <: OrdinaryDiffEqAlgorithm end +struct Anas5{T} <: OrdinaryDiffEqAlgorithm + w::T +end struct OwrenZen3 <: OrdinaryDiffEqAdaptiveAlgorithm end struct OwrenZen4 <: OrdinaryDiffEqAdaptiveAlgorithm end struct OwrenZen5 <: OrdinaryDiffEqAdaptiveAlgorithm end diff --git a/src/perform_step/fixed_timestep_perform_step.jl b/src/perform_step/fixed_timestep_perform_step.jl index 5874f304ec..9f6d2639a4 100644 --- a/src/perform_step/fixed_timestep_perform_step.jl +++ b/src/perform_step/fixed_timestep_perform_step.jl @@ -409,14 +409,19 @@ 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,a51,a52,a53,a54,a61,a62,a63,a64,a65,c1,c2,c3,c4,c5,c6,b1,b2,b3,b4,b5,b6 = cache + @unpack a21,a31,a32,a41,a42,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)) 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+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+b2*k2+b3*k3+b4*k4+b5*k5+b6*k6) + 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; diff --git a/src/tableaus/low_order_rk_tableaus.jl b/src/tableaus/low_order_rk_tableaus.jl index 191ac55daa..77ba342fcd 100644 --- a/src/tableaus/low_order_rk_tableaus.jl +++ b/src/tableaus/low_order_rk_tableaus.jl @@ -1408,14 +1408,12 @@ struct Anas5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache a63::T a64::T a65::T - c1::T2 c2::T2 c3::T2 c4::T2 c5::T2 c6::T2 b1::T - b2::T b3::T b4::T b5::T @@ -1438,48 +1436,44 @@ Base.@pure function Anas5ConstantCache{T<:CompiledFloats,T2<:CompiledFloats}(::T a63 = T(1.89) a64 = T(-2.295) a65 = T(1) - c1 = T2(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) - b2 = T(0.0) 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,c1,c2,c3,c4,c5,c6,b1,b2,b3,b4,b5,b6) + 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) - a51 = T(-11277/8000) - a52 = T(171/80) - a53 = T(-459/2000) - a54 = T(3213/8000) - a61 = T(-519/100) - a62 = T(5) - a63 = T(189/100) - a64 = T(-459/200) - a65 = T(1) - c1 = T2(0) - c2 = T2(1/10) - c3 = T2(1/3) - c4 = T2(2/3) - c5 = T2(9/10) - c6 = T2(1) - b1 = T(23/216) - b2 = T(0) - 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,c1,c2,c3,c4,c5,c6,b1,b2,b3,b4,b5,b6) + 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(-519//100) + a62 = T(5//1) + a63 = T(189//100) + a64 = T(-459//200) + a65 = T(1//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 diff --git a/test/ode/ode_convergence_tests.jl b/test/ode/ode_convergence_tests.jl index 0d08c8620c..602c09b097 100644 --- a/test/ode/ode_convergence_tests.jl +++ b/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=5)) + @test abs(sim160.𝒪est[:final]-5) < testTol println("Stiff Solvers") From 2511f505f7711d2fc8613b2fa39601d965a77508 Mon Sep 17 00:00:00 2001 From: Lucas Booher Date: Thu, 28 Jun 2018 22:19:31 -0700 Subject: [PATCH 3/9] Finishing up the implementation. --- src/algorithms.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index fb0b326e4a..687427207c 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -38,8 +38,10 @@ struct Ralston <: OrdinaryDiffEqAdaptiveAlgorithm end struct Midpoint <: OrdinaryDiffEqAdaptiveAlgorithm end struct RK4 <: OrdinaryDiffEqAdaptiveAlgorithm end struct Anas5{T} <: OrdinaryDiffEqAlgorithm - w::T + w::T end +Base.@pure Anas5(; w=1) = Anas5(w) + struct OwrenZen3 <: OrdinaryDiffEqAdaptiveAlgorithm end struct OwrenZen4 <: OrdinaryDiffEqAdaptiveAlgorithm end struct OwrenZen5 <: OrdinaryDiffEqAdaptiveAlgorithm end From 8606aa3a2ea291a5020397857d81a165c250ff12 Mon Sep 17 00:00:00 2001 From: Lucas Booher Date: Thu, 28 Jun 2018 23:46:56 -0700 Subject: [PATCH 4/9] Fixing some mistakes --- src/perform_step/fixed_timestep_perform_step.jl | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/perform_step/fixed_timestep_perform_step.jl b/src/perform_step/fixed_timestep_perform_step.jl index 9f6d2639a4..c9ba785ca8 100644 --- a/src/perform_step/fixed_timestep_perform_step.jl +++ b/src/perform_step/fixed_timestep_perform_step.jl @@ -394,7 +394,7 @@ end end function initialize!(integrator, cache::Anas5ConstantCache) - integrator.kshortsize = 6 + integrator.kshortsize = 7 integrator.k = typeof(integrator.k)(integrator.kshortsize) integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal @@ -409,12 +409,12 @@ 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,a51,a52,a53,a54,a61,a62,a63,a64,c2,c3,c4,c5,c6,b1,b3,b4,b5,b6 = 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 ## 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)) + 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)) 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) @@ -429,7 +429,7 @@ end end function initialize!(integrator, cache::Anas5Cache) - integrator.kshortsize = 6 + 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; @@ -443,7 +443,11 @@ end @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,a65,c1,c2,c3,c4,c5,c6,b1,b2,b3,b4,b5,b6,b7 = cache.tab + @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)) @tight_loop_macros for i in uidx @inbounds tmp[i] = uprev[i]+a21*k1[i] end @@ -465,7 +469,7 @@ end end f(k6, tmp, p, t+c6*dt) @tight_loop_macros for i in uidx - @inbounds u[i] = uprev[i]+dt*(b1*k1[i]+b2*k2[i]+b3*k3[i]+b4*k4[i]+b5*k5[i]+b6*k6[i]) + @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 From b5b4bcd3914b038d688e11458a38da6ff9764310 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 10 Jul 2018 21:22:33 -0700 Subject: [PATCH 5/9] eps scaling in time cutoffs --- src/integrators/integrator_utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/integrators/integrator_utils.jl b/src/integrators/integrator_utils.jl index 538697136e..89c02246ae 100644 --- a/src/integrators/integrator_utils.jl +++ b/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 From 12fcc1b865e7a425d09738b7adf62e3c7efd0d80 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 11 Jul 2018 02:41:03 -0700 Subject: [PATCH 6/9] Fixing mistakes related to the constants a61 to a64 Fixing error in perform step. This has greatly increased the order. Update ode_convergence_tests.jl Remove DS_Store --- src/perform_step/fixed_timestep_perform_step.jl | 14 ++++++++++---- src/tableaus/low_order_rk_tableaus.jl | 16 ++++++++-------- test/ode/ode_convergence_tests.jl | 4 ++-- 3 files changed, 20 insertions(+), 14 deletions(-) diff --git a/src/perform_step/fixed_timestep_perform_step.jl b/src/perform_step/fixed_timestep_perform_step.jl index c9ba785ca8..8de2f1c382 100644 --- a/src/perform_step/fixed_timestep_perform_step.jl +++ b/src/perform_step/fixed_timestep_perform_step.jl @@ -414,11 +414,14 @@ end 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)) + 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+k3), p, t+c4*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) @@ -447,7 +450,10 @@ end 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)) + 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 @@ -457,7 +463,7 @@ end 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]+k3[i]) + @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 diff --git a/src/tableaus/low_order_rk_tableaus.jl b/src/tableaus/low_order_rk_tableaus.jl index 77ba342fcd..059970573a 100644 --- a/src/tableaus/low_order_rk_tableaus.jl +++ b/src/tableaus/low_order_rk_tableaus.jl @@ -1431,11 +1431,11 @@ Base.@pure function Anas5ConstantCache{T<:CompiledFloats,T2<:CompiledFloats}(::T a52 = T(2.1375) a53 = T(-0.2295) a54 = T(0.401625) - a61 = T(-5.19) + a61 = T(-4.0) a62 = T(5.0) - a63 = T(1.89) - a64 = T(-2.295) - a65 = T(1) + a63 = T(0.0) + a64 = T(0.0) + a65 = T(0.0) c2 = T2(0.1) c3 = T2(0.3333333333333333) c4 = T2(0.6666666666666667) @@ -1460,11 +1460,11 @@ Base.@pure function Anas5ConstantCache{T,T2}(::Type{T},::Type{T2}) a52 = T(171//80) a53 = T(-459//2000) a54 = T(3213//8000) - a61 = T(-519//100) + a61 = T(-4//1) a62 = T(5//1) - a63 = T(189//100) - a64 = T(-459//200) - a65 = T(1//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) diff --git a/test/ode/ode_convergence_tests.jl b/test/ode/ode_convergence_tests.jl index 602c09b097..4582a42ed4 100644 --- a/test/ode/ode_convergence_tests.jl +++ b/test/ode/ode_convergence_tests.jl @@ -51,8 +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=5)) - @test abs(sim160.𝒪est[:final]-5) < testTol + sim160 = test_convergence(dts,prob,Anas5(w=2)) + @test abs(sim160.𝒪est[:l2]-4) < testTol println("Stiff Solvers") From 1126ced4177c4d41d3849fa12674ec7de861624a Mon Sep 17 00:00:00 2001 From: Unknown Date: Fri, 13 Jul 2018 02:22:22 +0530 Subject: [PATCH 7/9] Change save_start to take saveat into consideration --- src/solve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solve.jl b/src/solve.jl index cfd582115c..8195909dbf 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -19,7 +19,7 @@ function init{algType<:OrdinaryDiffEqAlgorithm,recompile_flag}( 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 = isempty(saveat) ? true : prob.tspan[1] in saveat,save_end = true, callback=nothing, dense = save_everystep && !(typeof(alg) <: FunctionMap), calck = (callback != nothing && callback != CallbackSet()) || # Empty callback From be73bb131337141bd9cbe289c0ba986048c2ab16 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 12 Jul 2018 14:39:16 -0700 Subject: [PATCH 8/9] Update ode_convergence_tests.jl --- test/ode/ode_convergence_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/ode/ode_convergence_tests.jl b/test/ode/ode_convergence_tests.jl index 041634a47e..75e2d6b790 100644 --- a/test/ode/ode_convergence_tests.jl +++ b/test/ode/ode_convergence_tests.jl @@ -52,7 +52,7 @@ for i = 1:2 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) < testTol + @test abs(sim160.𝒪est[:l2]-4) < 2*testTol println("Stiff Solvers") From 29aeafce652f9f1df8f62f41d5b12f254cec52c2 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 12 Jul 2018 19:58:17 -0700 Subject: [PATCH 9/9] update saveat tests and fix saveat behavior --- src/solve.jl | 6 +++-- test/ode/ode_saveat_tests.jl | 52 +++++++++++++++++++----------------- 2 files changed, 31 insertions(+), 27 deletions(-) diff --git a/src/solve.jl b/src/solve.jl index 8195909dbf..e68be5ffee 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -19,9 +19,11 @@ function init{algType<:OrdinaryDiffEqAlgorithm,recompile_flag}( d_discontinuities= eltype(prob.tspan)[], save_idxs = nothing, save_everystep = isempty(saveat), - save_timeseries = nothing,save_start = isempty(saveat) ? true : prob.tspan[1] in saveat,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 diff --git a/test/ode/ode_saveat_tests.jl b/test/ode/ode_saveat_tests.jl index 49de600c30..f25685cd2c 100644 --- a/test/ode/ode_saveat_tests.jl +++ b/test/ode/ode_saveat_tests.jl @@ -3,65 +3,67 @@ using OrdinaryDiffEq, DiffEqProblemLibrary, Base.Test @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])