Merge pull request #714 from AlthausKonstantin/explicit_extrapolation

Adaptive explicit extrapolation

src/OrdinaryDiffEq.jl

@@ -227,7 +227,7 @@ module OrdinaryDiffEq
   export AutoSwitch, AutoTsit5, AutoDP5,
          AutoVern6, AutoVern7, AutoVern8, AutoVern9
 
-  export AitkenNeville
+  export AitkenNeville, ExtrapolationMidpointDeuflhard, ExtrapolationMidpointHairerWanner
 
   export KuttaPRK2p5
 

src/alg_utils.jl

@@ -116,6 +116,11 @@ get_current_alg_order(alg::QNDF,cache) = cache.order
 get_current_adaptive_order(alg::QNDF,cache) = cache.order
 get_current_adaptive_order(alg::OrdinaryDiffEqExtrapolationVarOrderVarStepAlgorithm,cache) = cache.cur_order
 get_current_alg_order(alg::OrdinaryDiffEqExtrapolationVarOrderVarStepAlgorithm,cache) = cache.cur_order
+get_current_alg_order(alg::ExtrapolationMidpointDeuflhard,cache) = 2(cache.n_curr + 1)
+get_current_adaptive_order(alg::ExtrapolationMidpointDeuflhard,cache) = 2cache.n_curr
+get_current_alg_order(alg::ExtrapolationMidpointHairerWanner,cache) = 2(cache.n_curr + 1)
+get_current_adaptive_order(alg::ExtrapolationMidpointHairerWanner,cache) = 2cache.n_curr
+
 
 #alg_adaptive_order(alg::OrdinaryDiffEqAdaptiveAlgorithm) = error("Algorithm is adaptive with no order")
 get_current_adaptive_order(alg::OrdinaryDiffEqAlgorithm,cache) = alg_adaptive_order(alg)
@@ -334,6 +339,8 @@ alg_order(alg::MEBDF2) = 2
 
 alg_maximum_order(alg) = alg_order(alg)
 alg_maximum_order(alg::CompositeAlgorithm) = maximum(alg_order(x) for x in alg.algs)
+alg_maximum_order(alg::ExtrapolationMidpointDeuflhard) = 2(alg.n_max+1)
+alg_maximum_order(alg::ExtrapolationMidpointHairerWanner) = 2(alg.n_max+1)
 
 alg_adaptive_order(alg::ExplicitRK) = alg.tableau.adaptiveorder
 alg_adaptive_order(alg::OrdinaryDiffEqAlgorithm) = alg_order(alg)-1
@@ -368,17 +375,23 @@ beta2_default(alg::FunctionMap) = 0
 beta2_default(alg::DP8) = 0//1
 beta2_default(alg::DP5) = 4//100
 beta2_default(alg::DP5Threaded) = 4//100
+beta2_default(alg::ExtrapolationMidpointDeuflhard) = 0//1
+beta2_default(alg::ExtrapolationMidpointHairerWanner) = 0//1
 
 beta1_default(alg::OrdinaryDiffEqAlgorithm,beta2) = 7//(10alg_order(alg))
 beta1_default(alg::FunctionMap,beta2) = 0
 beta1_default(alg::DP8,beta2) = typeof(beta2)(1//alg_order(alg)) - beta2/5
 beta1_default(alg::DP5,beta2) = typeof(beta2)(1//alg_order(alg)) - 3beta2/4
 beta1_default(alg::DP5Threaded,beta2) = typeof(beta2)(1//alg_order(alg)) - 3beta2/4
+beta1_default(alg::ExtrapolationMidpointDeuflhard,beta2) =  1//(2alg.n_init+1)
+beta1_default(alg::ExtrapolationMidpointHairerWanner,beta2) =  1//(2alg.n_init+1)
 
 gamma_default(alg::OrdinaryDiffEqAlgorithm) = 9//10
 gamma_default(alg::ESERK5) = 8//10
 gamma_default(alg::RKC) = 8//10
 gamma_default(alg::IRKC) = 8//10
+gamma_default(alg::ExtrapolationMidpointDeuflhard) = (1//4)^beta1_default(alg,beta2_default(alg))
+gamma_default(alg::ExtrapolationMidpointHairerWanner) = (65//100)^beta1_default(alg,beta2_default(alg))
 
 qsteady_min_default(alg::OrdinaryDiffEqAlgorithm) = 1
 qsteady_max_default(alg::OrdinaryDiffEqAlgorithm) = 1

src/algorithms.jl

@@ -46,6 +46,75 @@ struct AitkenNeville <: OrdinaryDiffEqExtrapolationVarOrderVarStepAlgorithm
 end
 AitkenNeville(;max_order=9,min_order=1,init_order=5) = AitkenNeville(max_order,min_order,init_order)
 
+struct ExtrapolationMidpointDeuflhard <: OrdinaryDiffEqExtrapolationVarOrderVarStepAlgorithm
+  n_min::Int # Minimal extrapolation order
+  n_init::Int # Initial extrapolation order
+  n_max::Int # Maximal extrapolation order
+  sequence_symbol::Symbol # Name of the subdividing sequence
+end
+function ExtrapolationMidpointDeuflhard(;min_extrapolation_order=1,init_extrapolation_order=5, max_extrapolation_order=10, sequence_symbol = :harmonic)
+  # Enforce 1 <=  min_extrapolation_order <= init_extrapolation_order <= max_extrapolation_order:
+  n_min = max(1,min_extrapolation_order)
+  n_init = max(n_min,init_extrapolation_order)
+  n_max = max(n_init,max_extrapolation_order)
+
+  # Warn user if orders have been changed
+  if (min_extrapolation_order, init_extrapolation_order, max_extrapolation_order) != (n_min,n_init,n_max)
+    @warn "The range of extrapolation orders and/or the initial order given to the
+      `ExtrapolationMidpointDeuflhard` algorithm are not valid and have been changed:
+      Minimal order: " * lpad(min_extrapolation_order,2," ") * " --> "  * lpad(n_min,2," ") * "
+      Maximal order: " * lpad(max_extrapolation_order,2," ") * " --> "  * lpad(n_max,2," ") * "
+      Initial order: " * lpad(init_extrapolation_order,2," ") * " --> "  * lpad(n_init,2," ")
+  end
+
+  # Warn user if sequence has been changed:
+  if sequence_symbol != :harmonic && sequence_symbol != :romberg && sequence_symbol != :bulirsch
+    @warn "The `sequence_symbol` given to the `ExtrapolationMidpointDeuflhard` algorithm
+       is not valid: it must match `:harmonic`, `:romberg` or `:bulirsch`.
+       Thus it has been changed
+      :$(sequence_symbol) --> :harmonic"
+    sequence_symbol = :harmonic
+  end
+
+  # Initialize algorithm
+  ExtrapolationMidpointDeuflhard(n_min,n_init,n_max,sequence_symbol)
+end
+
+struct ExtrapolationMidpointHairerWanner <: OrdinaryDiffEqExtrapolationVarOrderVarStepAlgorithm
+  n_min::Int # Minimal extrapolation order
+  n_init::Int # Initial extrapolation order
+  n_max::Int # Maximal extrapolation order
+  sequence_symbol::Symbol # Name of the subdividing sequence
+end
+function ExtrapolationMidpointHairerWanner(;min_extrapolation_order=2,init_extrapolation_order=5, max_extrapolation_order=10, sequence_symbol = :harmonic)
+  # Enforce 2 <=  min_extrapolation_order
+  # and min_extrapolation_order + 1 <= init_extrapolation_order <= max_extrapolation_order - 1:
+  n_min = max(2, min_extrapolation_order)
+  n_init = max(n_min, init_extrapolation_order)
+  n_max = max(n_init + 1, max_extrapolation_order)
+
+  # Warn user if orders have been changed
+  if (min_extrapolation_order, init_extrapolation_order, max_extrapolation_order) != (n_min,n_init,n_max)
+    @warn "The range of extrapolation orders and/or the initial order given to the
+      `ExtrapolationMidpointHairerWanner` algorithm are not valid and have been changed:
+      Minimal order: " * lpad(min_extrapolation_order,2," ") * " --> "  * lpad(n_min,2," ") * "
+      Maximal order: " * lpad(max_extrapolation_order,2," ") * " --> "  * lpad(n_max,2," ") * "
+      Initial order: " * lpad(init_extrapolation_order,2," ") * " --> "  * lpad(n_init,2," ")
+  end
+
+  # Warn user if sequence has been changed:
+  if sequence_symbol != :harmonic && sequence_symbol != :romberg && sequence_symbol != :bulirsch
+    @warn "The `sequence_symbol` given to the `ExtrapolationMidpointHairerWanner` algorithm
+       is not valid: it must match `:harmonic`, `:romberg` or `:bulirsch`.
+       Thus it has been changed
+      :$(sequence_symbol) --> :harmonic"
+    sequence_symbol = :harmonic
+  end
+
+  # Initialize algorithm
+  ExtrapolationMidpointHairerWanner(n_min,n_init,n_max,sequence_symbol)
+end
+
 struct RK46NL <: OrdinaryDiffEqAlgorithm end
 struct Heun <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct Ralston <: OrdinaryDiffEqAdaptiveAlgorithm end

src/caches/extrapolation_caches.jl

@@ -47,3 +47,213 @@ function alg_cache(alg::AitkenNeville,u,rate_prototype,uEltypeNoUnits,uBottomElt
   step_no = zero(Int)
   AitkenNevilleConstantCache(dtpropose,T,cur_order,work,A,step_no)
 end
+
+
+struct extrapolation_coefficients
+  # This structure is used by the caches of the algorithms
+  # ExtrapolationMidpointDeuflhard() and  ExtrapolationMidpointHairerWanner().
+  # It contains the constant coefficients used to extrapolate the internal discretisations
+  # in their perfom_step! function and some additional constant data.
+
+  subdividing_sequence::Array{BigInt,1}  # subdividing_sequence[n] is used for the (n -1)th internal discretisation
+
+  # Weights and Scaling factors for extrapolation operators
+  extrapolation_weights::Array{Rational{BigInt},2}
+  extrapolation_scalars::Array{Rational{BigInt},1}
+
+  # Weights and scaling factors for internal extrapolation operators (used for error estimate)
+  extrapolation_weights_2::Array{Rational{BigInt},2}
+  extrapolation_scalars_2::Array{Rational{BigInt},1}
+end
+
+function create_extrapolation_coefficients(alg::algType) where {algType <: Union{ExtrapolationMidpointDeuflhard, ExtrapolationMidpointHairerWanner}}
+  # Compute and return extrapolation_coefficients
+
+  @unpack n_min, n_init, n_max, sequence_symbol = alg
+
+  # Initialize subdividing_sequence:
+  if sequence_symbol == :harmonic
+      subdividing_sequence = [BigInt(n+1) for n = 0:n_max]
+  elseif sequence_symbol == :romberg
+      subdividing_sequence = [BigInt(2)^n for n = 0:n_max]
+  else # sequence_symbol == :bulirsch
+      subdividing_sequence = [n==0 ? BigInt(1) : (isodd(n) ? BigInt(2)^Int64(n/2 + 0.5) : 3BigInt(2^Int64(n/2 - 1))) for n = 0:n_max]
+  end
+
+
+  # Compute nodes corresponding to subdividing_sequence
+  nodes = BigInt(1) .// subdividing_sequence .^ 2
+
+  # Compute barycentric weights for internal extrapolation operators
+  extrapolation_weights_2 = zeros(Rational{BigInt}, n_max, n_max)
+  extrapolation_weights_2[1,:] = ones(Rational{BigInt}, 1, n_max)
+  for n = 2:n_max
+      distance = nodes[2:n] .- nodes[n+1]
+      extrapolation_weights_2[1:(n-1), n] = extrapolation_weights_2[1:n-1, n-1] .// distance
+      extrapolation_weights_2[n, n] = 1 // prod(-distance)
+  end
+
+  # Compute barycentric weights for extrapolation operators
+  extrapolation_weights = zeros(Rational{BigInt}, n_max+1, n_max+1)
+  for n = 1:n_max
+      extrapolation_weights[n+1, (n+1) : (n_max+1)] = extrapolation_weights_2[n, n:n_max] // (nodes[n+1] - nodes[1])
+      extrapolation_weights[1, n] = 1 // prod(nodes[1] .- nodes[2:n])
+  end
+  extrapolation_weights[1, n_max+1] = 1 // prod(nodes[1] .- nodes[2:n_max+1])
+
+  # Rescale barycentric weights to obtain weights of 1. Barycentric Formula
+  for m = 1:(n_max+1)
+      extrapolation_weights[1:m, m] = - extrapolation_weights[1:m, m] .// nodes[1:m]
+      if 2 <= m
+          extrapolation_weights_2[1:m-1, m-1] = -extrapolation_weights_2[1:m-1, m-1] .// nodes[2:m]
+      end
+  end
+
+  # Compute scaling factors for internal extrapolation operators
+  extrapolation_scalars_2 = ones(Rational{BigInt}, n_max)
+  extrapolation_scalars_2[1] = -nodes[2]
+  for n = 1:(n_max-1)
+      extrapolation_scalars_2[n+1] = -extrapolation_scalars_2[n] * nodes[n+2]
+  end
+
+  # Compute scaling factors for extrapolation operators
+  extrapolation_scalars = -nodes[1] * [BigInt(1); extrapolation_scalars_2]
+
+  # Initialize and return extrapolation_coefficients
+  extrapolation_coefficients(subdividing_sequence,
+      extrapolation_weights, extrapolation_scalars,
+      extrapolation_weights_2, extrapolation_scalars_2)
+end
+
+@cache mutable struct ExtrapolationMidpointDeuflhardConstantCache{QType, extrapolation_coefficients} <: OrdinaryDiffEqConstantCache
+  # Values that are mutated
+  Q::Vector{QType} # Storage for stepsize scaling factors. Q[n] contains information for extrapolation order (n + alg.n_min - 1)
+  n_curr::Int64 # Storage for the current extrapolation order
+  n_old::Int64 # Storage for the extrapolation order n_curr before perfom_step! changes the latter
+
+  # Constant values
+  coefficients::extrapolation_coefficients
+  stage_number::Vector{Int64} # stage_number[n] contains information for extrapolation order (n + alg.n_min - 1)
+end
+
+function alg_cache(alg::ExtrapolationMidpointDeuflhard,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
+    # Initialize cache's members
+    QType = tTypeNoUnits <: Integer ? typeof(qmin_default(alg)) : tTypeNoUnits # Cf. DiffEqBase.__init in solve.jl
+
+    Q = fill(zero(QType),alg.n_max - alg.n_min + 1)
+    n_curr = alg.n_init
+    n_old = alg.n_init
+
+    coefficients = create_extrapolation_coefficients(alg)
+    stage_number = [2sum(Int64.(coefficients.subdividing_sequence[1:n+1])) - n for n = alg.n_min:alg.n_max]
+
+    # Initialize cache
+    ExtrapolationMidpointDeuflhardConstantCache(Q, n_curr, n_old, coefficients, stage_number)
+end
+
+
+
+@cache mutable struct ExtrapolationMidpointDeuflhardCache{uType,uNoUnitsType,rateType,QType,extrapolation_coefficients} <: OrdinaryDiffEqMutableCache
+  # Values that are mutated
+  utilde::uType
+  u_temp1::uType
+  u_temp2::uType
+  tmp::uType # for get_tmp_cache()
+  T::Array{uType,1}  # Storage for the internal discretisations obtained by the explicit midpoint rule
+  res::uNoUnitsType # Storage for the scaled residual of u and utilde
+
+  fsalfirst::rateType
+  k::rateType
+
+  # Constant values
+  Q::Vector{QType} # Storage for stepsize scaling factors. Q[n] contains information for extrapolation order (n + alg.n_min - 1)
+  n_curr::Int64 # Storage for the current extrapolation order
+  n_old::Int64 # Storage for the extrapolation order n_curr before perfom_step! changes the latter
+  coefficients::extrapolation_coefficients
+  stage_number::Vector{Int64} # Stage_number[n] contains information for extrapolation order (n + alg.n_min - 1)
+end
+
+function alg_cache(alg::ExtrapolationMidpointDeuflhard,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
+  # Initialize cache's members
+  utilde = zero(u)
+  u_temp1 = zero(u)
+  u_temp2 = zero(u)
+  tmp = zero(u)
+  T = fill(zero(u), alg.n_max + 1)
+  res = uEltypeNoUnits.(zero(u))
+
+  fsalfirst = zero(rate_prototype)
+  k = zero(rate_prototype)
+
+  cc =  alg_cache(alg::ExtrapolationMidpointDeuflhard,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,Val{false})
+  # Initialize cache
+  ExtrapolationMidpointDeuflhardCache(utilde, u_temp1, u_temp2, tmp, T, res, fsalfirst, k,cc.Q, cc.n_curr, cc.n_old, cc.coefficients,cc.stage_number)
+end
+
+@cache mutable struct ExtrapolationMidpointHairerWannerConstantCache{QType,extrapolation_coefficients} <: OrdinaryDiffEqConstantCache
+  # Values that are mutated
+  Q::Vector{QType} # Storage for stepsize scaling factors. Q[n] contains information for extrapolation order (n - 1)
+  n_curr::Int64 # Storage for the current extrapolation order
+  n_old::Int64 # Storage for the extrapolation order n_curr before perfom_step! changes the latter
+
+  # Constant values
+  coefficients::extrapolation_coefficients
+  stage_number::Vector{Int64} # stage_number[n] contains information for extrapolation order (n - 1)
+  sigma::Rational{Int64} # Parameter for order selection
+end
+
+function alg_cache(alg::ExtrapolationMidpointHairerWanner,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
+  # Initialize cache's members
+  QType = tTypeNoUnits <: Integer ? typeof(qmin_default(alg)) : tTypeNoUnits # Cf. DiffEqBase.__init in solve.jl
+
+  Q = fill(zero(QType),alg.n_max + 1)
+  n_curr = alg.n_init
+  n_old = alg.n_init
+
+  coefficients = create_extrapolation_coefficients(alg)
+  stage_number = [2sum(Int64.(coefficients.subdividing_sequence[1:n+1])) - n for n = 0:alg.n_max]
+  sigma = 9//10
+
+  # Initialize the constant cache
+  ExtrapolationMidpointHairerWannerConstantCache(Q, n_curr, n_old, coefficients, stage_number, sigma)
+end
+
+@cache mutable struct ExtrapolationMidpointHairerWannerCache{uType,uNoUnitsType,rateType,QType,extrapolation_coefficients} <: OrdinaryDiffEqMutableCache
+  # Values that are mutated
+  utilde::uType
+  u_temp1::uType
+  u_temp2::uType
+  tmp::uType # for get_tmp_cache()
+  T::Array{uType,1}  # Storage for the internal discretisations obtained by the explicit midpoint rule
+  res::uNoUnitsType # Storage for the scaled residual of u and utilde
+
+  fsalfirst::rateType
+  k::rateType
+
+  # Constant values
+  Q::Vector{QType} # Storage for stepsize scaling factors. Q[n] contains information for extrapolation order (n - 1)
+  n_curr::Int64 # Storage for the current extrapolation order
+  n_old::Int64 # Storage for the extrapolation order n_curr before perfom_step! changes the latter
+  coefficients::extrapolation_coefficients
+  stage_number::Vector{Int64} # stage_number[n] contains information for extrapolation order (n - 1)
+  sigma::Rational{Int64} # Parameter for order selection
+end
+
+
+function alg_cache(alg::ExtrapolationMidpointHairerWanner,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
+  # Initialize cache's members
+  utilde = zero(u)
+  u_temp1 = zero(u)
+  u_temp2 = zero(u)
+  tmp = zero(u)
+  T = fill(zero(u), alg.n_max + 1)
+  res = uEltypeNoUnits.(zero(u))
+  fsalfirst = zero(rate_prototype)
+  k = zero(rate_prototype)
+
+  cc = alg_cache(alg,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,Val{false})
+
+  # Initialize the cache
+  ExtrapolationMidpointHairerWannerCache(utilde, u_temp1, u_temp2, tmp, T, res, fsalfirst, k,
+      cc.Q, cc.n_curr, cc.n_old, cc.coefficients, cc.stage_number, cc.sigma)
+end