From ac1d8fd251e9ecb3933a4b1c85325a1dc07bae63 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 5 Jul 2017 16:36:41 +0200 Subject: [PATCH] Initial integration of Anderson acceleration --- src/DelayDiffEq.jl | 5 +- src/algorithms.jl | 20 ++++---- src/integrator_interface.jl | 99 +++++++++++++++++++++++++------------ src/integrator_type.jl | 27 ++++++---- src/iteration_function.jl | 41 +++++++++++++++ src/solve.jl | 37 ++++++++------ test/unconstrained.jl | 20 ++++---- test/units.jl | 28 +++++------ 8 files changed, 188 insertions(+), 89 deletions(-) create mode 100644 src/iteration_function.jl diff --git a/src/DelayDiffEq.jl b/src/DelayDiffEq.jl index d92715ec..18f31ce2 100644 --- a/src/DelayDiffEq.jl +++ b/src/DelayDiffEq.jl @@ -9,7 +9,7 @@ using OrdinaryDiffEq, DataStructures, RecursiveArrayTools, Combinatorics using Compat -using Base.Test +using NLsolve import OrdinaryDiffEq: initialize!, perform_step!, loopfooter!, loopheader!, alg_order, handle_tstop!, ODEIntegrator, savevalues!, @@ -22,6 +22,9 @@ import DiffEqBase: solve, solve!, init, resize!, u_cache, user_cache, import OrdinaryDiffEq: Rosenbrock23Cache, Rosenbrock32Cache, ImplicitEulerCache, TrapezoidCache +import NLsolve: AbstractDifferentiableMultivariateFunction + +include("iteration_function.jl") include("integrator_type.jl") include("integrator_interface.jl") include("history_function.jl") diff --git a/src/algorithms.jl b/src/algorithms.jl index 9fe2b1e4..36c8b572 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -3,15 +3,17 @@ immutable MethodOfSteps{algType,AType,RType,NType,constrained} <: AbstractMethodOfStepsAlgorithm{constrained} alg::algType - picardabstol::AType - picardreltol::RType - picardnorm::NType - max_picard_iters::Int + fixedpoint_abstol::AType + fixedpoint_reltol::RType + max_fixedpoint_iters::Int + picardnorm::NType # anderson acceleration with nlsolve always uses infinite norm of residuals + m::Int # controls history size of anderson acceleration (m=0 corresponds to simple fixed-point iteration) end Base.@pure MethodOfSteps(alg;constrained=false, - picardabstol = nothing, - picardreltol = nothing, - picardnorm = nothing, - max_picard_iters = 10) = - MethodOfSteps{typeof(alg),typeof(picardabstol),typeof(picardreltol),typeof(picardnorm),constrained}(alg,picardabstol,picardreltol,picardnorm,max_picard_iters) + fixedpoint_abstol = nothing, + fixedpoint_reltol = nothing, + max_fixedpoint_iters = 10, + picardnorm = nothing, + m = 0) = + MethodOfSteps{typeof(alg),typeof(fixedpoint_abstol),typeof(fixedpoint_reltol),typeof(picardnorm),constrained}(alg,fixedpoint_abstol,fixedpoint_reltol,max_fixedpoint_iters,picardnorm,m) diff --git a/src/integrator_interface.jl b/src/integrator_interface.jl index ee0c2fa2..b1bd9c48 100644 --- a/src/integrator_interface.jl +++ b/src/integrator_interface.jl @@ -1,67 +1,102 @@ function savevalues!(integrator::DDEIntegrator,force_save=false) + # update ODE integrator with values of DDE integrator integrator.integrator.u = integrator.u integrator.integrator.k = integrator.k integrator.integrator.t = integrator.t + + # add steps for interpolation when needed OrdinaryDiffEq.ode_addsteps!(integrator.integrator,integrator.f) + + # update solution of ODE integrator savevalues!(integrator.integrator,force_save) end function postamble!(integrator::DDEIntegrator) + # update ODE integrator with values of DDE integrator integrator.integrator.u = integrator.u integrator.integrator.k = integrator.k integrator.integrator.t = integrator.t + + # clean up solution of ODE integrator OrdinaryDiffEq.postamble!(integrator.integrator) end function perform_step!(integrator::DDEIntegrator) integrator.tprev = integrator.t # this is necessary to extrapolate from the current interval + + # update ODE integrator with values of DDE integrator integrator.integrator.uprev = integrator.uprev integrator.integrator.tprev = integrator.tprev integrator.integrator.fsalfirst = integrator.fsalfirst integrator.integrator.t = integrator.t integrator.integrator.dt = integrator.dt - # if dt>max lag, then it's explicit so use Picard iteration + # if dt>min lag, then it's explicit so use fixed-point iteration if integrator.dt >minimum(integrator.prob.lags) - - # the is done to correct the extrapolation - t_prev_cache = integrator.tprev + # this is done to correct the extrapolation t_cache = integrator.t - uprev_cache = integrator.uprev + tprev_cache = integrator.tprev + if typeof(integrator.uprev_cache) <: AbstractArray + copy!(integrator.uprev_cache,integrator.uprev) + else + integrator.uprev_cache = integrator.uprev + end - numiters = 1 - while true - if typeof(integrator.u) <: AbstractArray - copy!(integrator.u_cache,integrator.u) - else - integrator.u_cache = integrator.u - end - perform_step!(integrator,integrator.cache) + if eltype(integrator.u) <: Real # not possible to use nlsolve with units (yet) + integrator.first_iteration = true # first iteration step needs special updates - if typeof(integrator.resid) <: AbstractArray - integrator.resid .= (integrator.u .- integrator.u_cache)./(integrator.picardabstol .+ max.(abs.(integrator.u),abs.(integrator.u_cache)).*integrator.picardreltol) + # execute Anderson acceleration of fixed-point iteration + if typeof(integrator.u) <: Vector + nlsolve(integrator.iterator,integrator.u; + method=:anderson,m=integrator.m,iterations=integrator.max_fixedpoint_iters,xtol=0,ftol=1) else - integrator.resid = (integrator.u .- integrator.u_cache)./(integrator.picardabstol .+ max.(abs.(integrator.u),abs.(integrator.u_cache)).*integrator.picardreltol) - end - - picardEEst = integrator.picardnorm(integrator.resid) - if picardEEst < 1 || numiters > integrator.max_picard_iters - break + nlsolve(integrator.iterator,vec(collect(integrator.u)); + method=:anderson,m=integrator.m,iterations=integrator.max_fixedpoint_iters,xtol=0,ftol=1) end - if numiters == 1 - integrator.integrator.tprev = integrator.t - integrator.integrator.t = integrator.t+integrator.dt - integrator.integrator.uprev = integrator.u + else # simple Picard iteration for units + for i in 1:integrator.max_fixedpoint_iters + # update cache of u + if typeof(integrator.u) <: AbstractArray + copy!(integrator.u_cache,integrator.u) + else + integrator.u_cache = integrator.u + end + + perform_step!(integrator,integrator.cache) + + # calculate residuals + if typeof(integrator.resid) <: AbstractArray + @. integrator.resid = (integrator.u - integrator.u_cache)/ + @muladd(integrator.fixedpoint_abstol + max(abs(integrator.u),abs(integrator.u_cache))*integrator.fixedpoint_reltol) + else + integrator.resid = @. (integrator.u - integrator.u_cache)/ + (integrator.fixedpoint_abstol + max(abs(integrator.u),abs(integrator.u_cache))*integrator.fixedpoint_reltol) + end + + # special updates of ODE integrator in first iteration step + if i == 1 + integrator.integrator.tprev = integrator.t + integrator.integrator.t = integrator.t+integrator.dt + integrator.integrator.uprev = integrator.u + end + + # update ODE integrator with values of DDE integrator + integrator.integrator.u = integrator.u + integrator.integrator.k = integrator.k + + # stop fixed-point iteration when residuals are small + integrator.picardnorm(integrator.resid) < 1 && break end - integrator.integrator.u = integrator.u - integrator.integrator.k = integrator.k - numiters += 1 end + + # reset to cached values integrator.t = t_cache - integrator.integrator.t = t_cache - integrator.tprev = t_prev_cache - integrator.uprev = uprev_cache - + integrator.tprev = tprev_cache + if typeof(integrator.uprev) <: AbstractArray + copy!(integrator.uprev,integrator.uprev_cache) + else + integrator.uprev = integrator.uprev_cache + end else # no iterations perform_step!(integrator,integrator.cache) end diff --git a/src/integrator_type.jl b/src/integrator_type.jl index 00ec3846..a29a0986 100644 --- a/src/integrator_type.jl +++ b/src/integrator_type.jl @@ -9,11 +9,13 @@ type DDEIntegrator{algType<:OrdinaryDiffEqAlgorithm,uType,tType,absType,relType, uprev::uType tprev::tType u_cache::uType - picardabstol::absType - picardreltol::relType + uprev_cache::uType + fixedpoint_abstol::absType + fixedpoint_reltol::relType resid::residType # This would have to resize for resizing DDE to work picardnorm::NType - max_picard_iters::Int + max_fixedpoint_iters::Int + m::Int alg::algType rate_prototype::rateType notsaveat_idxs::Vector{Int} @@ -39,26 +41,33 @@ type DDEIntegrator{algType<:OrdinaryDiffEqAlgorithm,uType,tType,absType,relType, integrator::IType fsalfirst::rateType fsallast::rateType + first_iteration::Bool + iterator::IterationFunction{DDEIntegrator{algType,uType,tType,absType,relType,residType,tTypeNoUnits,tdirType,ksEltype,SolType,rateType,F,ProgressType,CacheType,IType,ProbType,NType,O}} (::Type{DDEIntegrator{algType,uType,tType,absType,relType, residType,tTypeNoUnits,tdirType,ksEltype,SolType,rateType,F, ProgressType,CacheType,IType,ProbType,NType,O}}){ algType,uType,tType,absType, relType,residType,tTypeNoUnits,tdirType,ksEltype,SolType,rateType,F,ProgressType, - CacheType,IType,ProbType,NType,O}(sol,prob,u,k,t,dt,f,uprev,tprev,u_cache, - picardabstol,picardreltol,resid,picardnorm,max_picard_iters, + CacheType,IType,ProbType,NType,O}(sol,prob,u,k,t,dt,f,uprev,tprev,u_cache,uprev_cache, + fixedpoint_abstol,fixedpoint_reltol,resid,picardnorm,max_fixedpoint_iters,m, alg,rate_prototype,notsaveat_idxs,dtcache,dtchangeable,dtpropose,tdir, EEst,qold,q11, iter,saveiter,saveiter_dense,prog,cache, kshortsize,just_hit_tstop,accept_step,isout,reeval_fsal,u_modified, - integrator,opts) = new{algType,uType,tType,absType,relType,residType,tTypeNoUnits,tdirType, + opts,integrator) = begin + dde_int = new{algType,uType,tType,absType,relType,residType,tTypeNoUnits,tdirType, ksEltype,SolType,rateType,F,ProgressType,CacheType,IType,ProbType,NType,O}( - sol,prob,u,k,t,dt,f,uprev,tprev,u_cache, - picardabstol,picardreltol,resid,picardnorm,max_picard_iters, + sol,prob,u,k,t,dt,f,uprev,tprev,u_cache,uprev_cache, + fixedpoint_abstol,fixedpoint_reltol,resid,picardnorm,max_fixedpoint_iters,m, alg,rate_prototype,notsaveat_idxs,dtcache,dtchangeable,dtpropose,tdir, EEst,qold,q11, iter,saveiter,saveiter_dense,prog,cache, - kshortsize,just_hit_tstop,accept_step,isout,reeval_fsal,u_modified,integrator,opts) # Leave off fsalfirst and last + kshortsize,just_hit_tstop,accept_step,isout,reeval_fsal,u_modified,opts,integrator) # Leave off fsalfirst, fsallast, first_iteration and iterator + dde_int.iterator = IterationFunction{DDEIntegrator{algType,uType,tType,absType,relType,residType,tTypeNoUnits,tdirType, + ksEltype,SolType,rateType,F,ProgressType,CacheType,IType,ProbType,NType,O}}(dde_int) + dde_int + end end function (integrator::DDEIntegrator)(t,deriv::Type=Val{0};idxs=nothing) diff --git a/src/iteration_function.jl b/src/iteration_function.jl new file mode 100644 index 00000000..9aa56009 --- /dev/null +++ b/src/iteration_function.jl @@ -0,0 +1,41 @@ +# nlsolve optimizes subtypes of AbstractDifferentiableMultivariateFunction +mutable struct IterationFunction{I<:AbstractDDEIntegrator} <: AbstractDifferentiableMultivariateFunction + integrator::I + f!::IterationFunction{I} # Anderson acceleration needs field f! + + function IterationFunction{I}(integrator::I) where I <:AbstractDDEIntegrator + f = new(integrator) + f.f! = f + end +end + +function (f::IterationFunction)(x,fvec) + integrator = f.integrator + + # update u + if typeof(integrator.u) <: AbstractArray + copy!(integrator.u,x) + else + integrator.u = x[1] # x is always a vector + end + + perform_step!(integrator,integrator.cache) + + # update output vector of residuals + # in contrast to picard iteration we can not use integrator.resid + # as output vector for different fixed-point iterations + @. fvec = (integrator.u - x)/ + @muladd(integrator.fixedpoint_abstol + max(abs(integrator.u),abs(x))*integrator.fixedpoint_reltol) + + # special updates of ODE integrator after the first iteration necessary + if integrator.first_iteration + integrator.integrator.tprev = integrator.t + integrator.integrator.t = integrator.t+integrator.dt + integrator.integrator.uprev = integrator.u + integrator.first_iteration = false + end + + # update ODE integrator with values of DDE integrator + integrator.integrator.u = integrator.u + integrator.integrator.k = integrator.k +end diff --git a/src/solve.jl b/src/solve.jl index 5de46acf..a7b02e59 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -73,10 +73,10 @@ function init{uType,tType,isinplace,algType<:AbstractMethodOfStepsAlgorithm,lTyp end integrator.dt = dt - if typeof(alg.picardabstol) <: Void - picardabstol_internal = map(eltype(uType),integrator.opts.abstol) + if typeof(alg.fixedpoint_abstol) <: Void + fixedpoint_abstol_internal = map(eltype(uType),integrator.opts.abstol) else - picardabstol_internal = map(eltype(uType),alg.picardabstol) + fixedpoint_abstol_internal = map(eltype(uType),alg.fixedpoint_abstol) end if typeof(alg.picardnorm) <: Void picardnorm = integrator.opts.internalnorm @@ -85,23 +85,32 @@ function init{uType,tType,isinplace,algType<:AbstractMethodOfStepsAlgorithm,lTyp uEltypeNoUnits = typeof(recursive_one(integrator.u)) - if typeof(alg.picardreltol) <: Void - picardreltol_internal = map(uEltypeNoUnits,integrator.opts.reltol) + if typeof(alg.fixedpoint_reltol) <: Void + fixedpoint_reltol_internal = map(uEltypeNoUnits,integrator.opts.reltol) else - picardreltol_internal = map(uEltypeNoUnits,alg.picardreltol) + fixedpoint_reltol_internal = map(uEltypeNoUnits,alg.fixedpoint_reltol) end if typeof(integrator.u) <: AbstractArray - resid = similar(integrator.u,uEltypeNoUnits) u_cache = similar(integrator.u) + uprev_cache = similar(integrator.u) else - resid = one(uEltypeNoUnits) u_cache = oneunit(eltype(uType)) + uprev_cache = oneunit(eltype(uType)) + end + + # for real numbers nlsolve is used for Anderson acceleration of fixed-point iteration which creates vectors of residuals + if eltype(integrator.u) <: Real + resid = nothing + elseif typeof(integrator.u) <: AbstractArray + resid = similar(integrator.u,uEltypeNoUnits) + else + resid = one(uEltypeNoUnits) end dde_int = DDEIntegrator{typeof(integrator.alg), uType,tType, - typeof(picardabstol_internal), - typeof(picardreltol_internal), + typeof(fixedpoint_abstol_internal), + typeof(fixedpoint_reltol_internal), typeof(resid), tTypeNoUnits,typeof(integrator.tdir), typeof(integrator.k),typeof(sol), @@ -112,9 +121,9 @@ function init{uType,tType,isinplace,algType<:AbstractMethodOfStepsAlgorithm,lTyp typeof(picardnorm), typeof(integrator.opts)}( sol,prob,integrator.u,integrator.k,integrator.t,integrator.dt, - dde_f2,integrator.uprev,integrator.tprev,u_cache, - picardabstol_internal,picardreltol_internal, - resid,picardnorm,alg.max_picard_iters, + dde_f2,integrator.uprev,integrator.tprev,u_cache,uprev_cache, + fixedpoint_abstol_internal,fixedpoint_reltol_internal, + resid,picardnorm,alg.max_fixedpoint_iters,alg.m, integrator.alg,integrator.rate_prototype,integrator.notsaveat_idxs,integrator.dtcache, integrator.dtchangeable,integrator.dtpropose,integrator.tdir, integrator.EEst,integrator.qold,integrator.q11, @@ -122,7 +131,7 @@ function init{uType,tType,isinplace,algType<:AbstractMethodOfStepsAlgorithm,lTyp integrator.saveiter_dense,integrator.prog,integrator.cache, integrator.kshortsize,integrator.just_hit_tstop,integrator.accept_step, integrator.isout, - integrator.reeval_fsal,integrator.u_modified,integrator.opts,integrator) # Leave off fsalfirst and last + integrator.reeval_fsal,integrator.u_modified,integrator.opts,integrator) # Leave off fsalfirst, fasllast, first_iteration, and iterator initialize!(dde_int) initialize!(integrator.opts.callback,integrator.t,integrator.u,dde_int) diff --git a/test/unconstrained.jl b/test/unconstrained.jl index 1af9d718..6d50a1ca 100644 --- a/test/unconstrained.jl +++ b/test/unconstrained.jl @@ -14,16 +14,16 @@ h = (t) -> 0.0 prob = ConstantLagDDEProblem(f,h,1.0,lags,(0.0,100.0)) -alg1 = MethodOfSteps(Tsit5(),constrained=false,max_picard_iters=100,picardabstol=1e-12,picardreltol=1e-12) +alg1 = MethodOfSteps(Tsit5(),constrained=false,max_fixedpoint_iters=100,fixedpoint_abstol=1e-12,fixedpoint_reltol=1e-12) sol1 = solve(prob,alg1) -alg2 = MethodOfSteps(DP8(),constrained=false,max_picard_iters=10,picardabstol=1e-8,picardreltol=1e-10) +alg2 = MethodOfSteps(DP8(),constrained=false,max_fixedpoint_iters=10,fixedpoint_abstol=1e-8,fixedpoint_reltol=1e-10) sol2 = solve(prob,alg2) -alg3 = MethodOfSteps(Tsit5(),constrained=true,max_picard_iters=10,picardabstol=1e-8,picardreltol=1e-10) +alg3 = MethodOfSteps(Tsit5(),constrained=true,max_fixedpoint_iters=10,fixedpoint_abstol=1e-8,fixedpoint_reltol=1e-10) sol3 = solve(prob,alg3) -alg4 = MethodOfSteps(DP5(),constrained=false,max_picard_iters=100,picardabstol=1e-12,picardreltol=1e-12) +alg4 = MethodOfSteps(DP5(),constrained=false,max_fixedpoint_iters=100,fixedpoint_abstol=1e-12,fixedpoint_reltol=1e-12) sol4 = solve(prob,alg4) @test abs(sol1[end] - sol2[end]) < 1e-3 @@ -38,16 +38,16 @@ h = (t) -> 0.0 prob = ConstantLagDDEProblem(f,h,1.0,lags,(0.0,10.0);iip=DiffEqBase.isinplace(f,4)) -alg1 = MethodOfSteps(Tsit5(),constrained=false,max_picard_iters=100,picardabstol=1e-12,picardreltol=1e-12) +alg1 = MethodOfSteps(Tsit5(),constrained=false,max_fixedpoint_iters=100,fixedpoint_abstol=1e-12,fixedpoint_reltol=1e-12) sol1 = solve(prob,alg1) -alg2 = MethodOfSteps(DP8(),constrained=false,max_picard_iters=10,picardabstol=1e-8,picardreltol=1e-10) +alg2 = MethodOfSteps(DP8(),constrained=false,max_fixedpoint_iters=10,fixedpoint_abstol=1e-8,fixedpoint_reltol=1e-10) sol2 = solve(prob,alg2) -alg3 = MethodOfSteps(Tsit5(),constrained=true,max_picard_iters=10,picardabstol=1e-8,picardreltol=1e-10) +alg3 = MethodOfSteps(Tsit5(),constrained=true,max_fixedpoint_iters=10,fixedpoint_abstol=1e-8,fixedpoint_reltol=1e-10) sol3 = solve(prob,alg3) -alg4 = MethodOfSteps(DP5(),constrained=false,max_picard_iters=100,picardabstol=1e-12,picardreltol=1e-12) +alg4 = MethodOfSteps(DP5(),constrained=false,max_fixedpoint_iters=100,fixedpoint_abstol=1e-12,fixedpoint_reltol=1e-12) sol4 = solve(prob,alg4) @test abs(sol1[end] - sol2[end]) < 1e-3 @@ -72,7 +72,7 @@ end prob = ConstantLagDDEProblem(f,h,[1.0],lags,(0.0,100.0)) -alg1 = MethodOfSteps(Tsit5(),constrained=false,max_picard_iters=100,picardabstol=1e-12,picardreltol=1e-12) +alg1 = MethodOfSteps(Tsit5(),constrained=false,max_fixedpoint_iters=100,fixedpoint_abstol=1e-12,fixedpoint_reltol=1e-12) @time sol1 = solve(prob,alg1) f = function (t,u,h,du) @@ -94,5 +94,5 @@ end prob = ConstantLagDDEProblem(f,h,[1.0],lags,(0.0,100.0)) -alg1 = MethodOfSteps(Tsit5(),constrained=false,max_picard_iters=100,picardabstol=1e-12,picardreltol=1e-12) +alg1 = MethodOfSteps(Tsit5(),constrained=false,max_fixedpoint_iters=100,fixedpoint_abstol=1e-12,fixedpoint_reltol=1e-12) @time sol1 = solve(prob,alg1) diff --git a/test/units.jl b/test/units.jl index 65c9e0ec..d38fde7a 100644 --- a/test/units.jl +++ b/test/units.jl @@ -11,37 +11,37 @@ h = (t) -> 0.0u"N" prob = ConstantLagDDEProblem(f,h,1.0u"N",lags,(0.0u"s",100.0u"s")) # Unconstrained algorithm without explicit absolut or relative tolerance -alg = MethodOfSteps(Tsit5(),constrained=false,max_picard_iters=100) +alg = MethodOfSteps(Tsit5(),constrained=false,max_fixedpoint_iters=100) solve(prob,alg) # Unconstrained algorithm without units -alg = MethodOfSteps(Tsit5(),constrained=false,max_picard_iters=100,picardabstol=1e-12,picardreltol=1e-12) +alg = MethodOfSteps(Tsit5(),constrained=false,max_fixedpoint_iters=100,fixedpoint_abstol=1e-12,fixedpoint_reltol=1e-12) sol = solve(prob,alg) # Unconstrained algorithm with correct units -alg2 = MethodOfSteps(Tsit5(),constrained=false,max_picard_iters=100,picardabstol=1e-9u"mN",picardreltol=1e-12) +alg2 = MethodOfSteps(Tsit5(),constrained=false,max_fixedpoint_iters=100,fixedpoint_abstol=1e-9u"mN",fixedpoint_reltol=1e-12) sol2 = solve(prob,alg2) @test sol.t == sol2.t && sol.u == sol2.u # Unconstrained algorithm with incorrect units for absolute tolerance -alg = MethodOfSteps(Tsit5(),constrained=false,max_picard_iters=100,picardabstol=1e-12u"s",picardreltol=1e-12) +alg = MethodOfSteps(Tsit5(),constrained=false,max_fixedpoint_iters=100,fixedpoint_abstol=1e-12u"s",fixedpoint_reltol=1e-12) @test_throws Unitful.DimensionError solve(prob,alg) # Unconstrained algorithm with units for relative tolerance -alg = MethodOfSteps(Tsit5(),constrained=false,max_picard_iters=100,picardabstol=1e-12,picardreltol=1e-12u"N") +alg = MethodOfSteps(Tsit5(),constrained=false,max_fixedpoint_iters=100,fixedpoint_abstol=1e-12,fixedpoint_reltol=1e-12u"N") @test_throws Unitful.DimensionError solve(prob,alg) # Constrained algorithm without explicit absolute or relative tolerance -alg = MethodOfSteps(Tsit5(),constrained=true,max_picard_iters=100) +alg = MethodOfSteps(Tsit5(),constrained=true,max_fixedpoint_iters=100) solve(prob,alg) # Constrained algorithm without units -alg = MethodOfSteps(Tsit5(),constrained=true,max_picard_iters=100,picardabstol=1e-12,picardreltol=1e-12) +alg = MethodOfSteps(Tsit5(),constrained=true,max_fixedpoint_iters=100,fixedpoint_abstol=1e-12,fixedpoint_reltol=1e-12) sol = solve(prob,alg) # Constrained algorithm with correct units -alg2 = MethodOfSteps(Tsit5(),constrained=true,max_picard_iters=100,picardabstol=1e-9u"mN",picardreltol=1e-12) +alg2 = MethodOfSteps(Tsit5(),constrained=true,max_fixedpoint_iters=100,fixedpoint_abstol=1e-9u"mN",fixedpoint_reltol=1e-12) sol2 = solve(prob,alg2) @test sol.t == sol2.t && sol.u == sol2.u @@ -55,29 +55,29 @@ h = (t) -> [0.0u"N"] prob = ConstantLagDDEProblem(f,h,[1.0u"N"],lags,(0.0u"s",100.0u"s")) # Unconstrained algorithm without explicit absolute or relative tolerance -alg = MethodOfSteps(Tsit5(),constrained=false,max_picard_iters=100) +alg = MethodOfSteps(Tsit5(),constrained=false,max_fixedpoint_iters=100) solve(prob,alg) # Unconstrained algorithm without units -alg = MethodOfSteps(Tsit5(),constrained=false,max_picard_iters=100,picardabstol=1e-12,picardreltol=1e-12) +alg = MethodOfSteps(Tsit5(),constrained=false,max_fixedpoint_iters=100,fixedpoint_abstol=1e-12,fixedpoint_reltol=1e-12) sol = solve(prob,alg) # Unconstrained algorithm with correct units and both absolute and relative tolerance as vector -alg2 = MethodOfSteps(Tsit5(),constrained=false,max_picard_iters=100,picardabstol=[1e-9u"mN"],picardreltol=[1e-12]) +alg2 = MethodOfSteps(Tsit5(),constrained=false,max_fixedpoint_iters=100,fixedpoint_abstol=[1e-9u"mN"],fixedpoint_reltol=[1e-12]) sol2 = solve(prob,alg2) @test sol.t == sol2.t && sol.u == sol2.u # Constrained algorithm without explicit absolute or relative tolerance -alg = MethodOfSteps(Tsit5(),constrained=true,max_picard_iters=100) +alg = MethodOfSteps(Tsit5(),constrained=true,max_fixedpoint_iters=100) solve(prob,alg) # Constrained algorithm without units -alg = MethodOfSteps(Tsit5(),constrained=true,max_picard_iters=100,picardabstol=1e-12,picardreltol=1e-12) +alg = MethodOfSteps(Tsit5(),constrained=true,max_fixedpoint_iters=100,fixedpoint_abstol=1e-12,fixedpoint_reltol=1e-12) sol = solve(prob,alg) # Constrained algorithm with correct units and both absolute and relative tolerance as vector -alg2 = MethodOfSteps(Tsit5(),constrained=true,max_picard_iters=100,picardabstol=[1e-9u"mN"],picardreltol=[1e-12]) +alg2 = MethodOfSteps(Tsit5(),constrained=true,max_fixedpoint_iters=100,fixedpoint_abstol=[1e-9u"mN"],fixedpoint_reltol=[1e-12]) sol2 = solve(prob,alg2) @test sol.t == sol2.t && sol.u == sol2.u