From 5928190bddd085098352015b4d23a7f503355789 Mon Sep 17 00:00:00 2001 From: James Gardner Date: Fri, 5 Feb 2021 15:22:50 +0000 Subject: [PATCH 1/6] Add BAOAB algorithm --- src/StochasticDiffEq.jl | 4 ++ src/alg_utils.jl | 3 ++ src/algorithms.jl | 10 +++++ src/caches/dynamical_caches.jl | 39 ++++++++++++++++ src/perform_step/dynamical.jl | 81 ++++++++++++++++++++++++++++++++++ src/solve.jl | 8 +++- test/sde/sde_dynamical.jl | 25 +++++++++++ 7 files changed, 168 insertions(+), 2 deletions(-) create mode 100644 src/caches/dynamical_caches.jl create mode 100644 src/perform_step/dynamical.jl create mode 100644 test/sde/sde_dynamical.jl diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl index d769c0752..b8e479732 100644 --- a/src/StochasticDiffEq.jl +++ b/src/StochasticDiffEq.jl @@ -84,6 +84,7 @@ module StochasticDiffEq include("caches/predcorr_caches.jl") include("caches/SROCK_caches.jl") include("caches/tau_caches.jl") + include("caches/dynamical_caches.jl") include("integrators/type.jl") include("dense.jl") include("alg_utils.jl") @@ -110,6 +111,7 @@ module StochasticDiffEq include("perform_step/composite.jl") include("perform_step/SROCK_perform_step.jl") include("perform_step/tau_leaping.jl") + include("perform_step/dynamical.jl") include("tableaus.jl") include("SROCK_tableaus.jl") include("iterated_integrals.jl") @@ -145,6 +147,8 @@ module StochasticDiffEq export TauLeaping, CaoTauLeaping + export BAOAB + export StochasticDiffEqRODEAlgorithm, StochasticDiffEqRODEAdaptiveAlgorithm, StochasticDiffEqRODECompositeAlgorithm diff --git a/src/alg_utils.jl b/src/alg_utils.jl index c99e1920e..1f1f1f5f0 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -95,6 +95,8 @@ alg_order(alg::SMEB) = 1//1 alg_order(alg::TauLeaping) = 1//1 alg_order(alg::CaoTauLeaping) = 1//1 +alg_order(alg::BAOAB) = 1//1 # Not sure what order it is + alg_order(alg::SKenCarp) = 2//1 alg_order(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = maximum(alg_order.(alg.algs)) get_current_alg_order(alg::StochasticDiffEqAlgorithm,cache) = alg_order(alg) @@ -196,6 +198,7 @@ alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::RKMilGeneral) = true alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::IIF1M) = true alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::IIF2M) = true alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = max((alg_compatible(prob,a) for a in alg.algs)...) +alg_compatible(prob::DiffEqBase.AbstractSDEProblem,alg::BAOAB) = is_diagonal_noise(prob) function alg_compatible(prob::JumpProblem,alg::Union{StochasticDiffEqJumpAdaptiveAlgorithm,StochasticDiffEqJumpAlgorithm}) prob.prob isa DiscreteProblem diff --git a/src/algorithms.jl b/src/algorithms.jl index fb2dbfb2f..ae67f9731 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -555,3 +555,13 @@ end struct RandomEM <: StochasticDiffEqRODEAlgorithm end const SplitSDEAlgorithms = Union{IIF1M,IIF2M,IIF1Mil,SKenCarp,SplitEM} + +""" +Leimkuhler B., Matthews C., Robust and efficient configurational molecular sampling via +Langevin dynamics, J. Chem. Phys. 138, 174102 (2013) +DOI:10.1063/1.4802990 +""" +struct BAOAB{T} <: StochasticDiffEqAlgorithm + gamma::T +end +BAOAB(;gamma=1.0) = BAOAB(gamma) diff --git a/src/caches/dynamical_caches.jl b/src/caches/dynamical_caches.jl new file mode 100644 index 000000000..5cfd45300 --- /dev/null +++ b/src/caches/dynamical_caches.jl @@ -0,0 +1,39 @@ + +struct BAOABConstantCache{uType,uEltypeNoUnits} <: StochasticDiffEqConstantCache + k::uType + half::uEltypeNoUnits + c1::uEltypeNoUnits + c2::uEltypeNoUnits +end +@cache struct BAOABCache{uType,uEltypeNoUnits,rateNoiseType} <: StochasticDiffEqMutableCache + utmp::uType + dutmp::uType + k::uType + gtmp::uType + noise::rateNoiseType + half::uEltypeNoUnits + c1::uEltypeNoUnits + c2::uEltypeNoUnits +end + +function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) + k = zero(rate_prototype.x[1]) + c1 = exp(-alg.gamma*dt) + c2 = sqrt(1 - c1^2) + BAOABConstantCache(k, uEltypeNoUnits(1//2), uEltypeNoUnits(c1), uEltypeNoUnits(c2)) +end + +function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}}) + dutmp = zero(u.x[1]) + utmp = zero(u.x[2]) + k = zero(rate_prototype.x[1]) + + gtmp = zero(rate_prototype.x[1]) + noise = zero(rate_prototype.x[1]) + + half = uEltypeNoUnits(1//2) + c1 = exp(-alg.gamma*dt) + c2 = sqrt(1 - c1^2) + + BAOABCache(utmp, dutmp, k, gtmp, noise, half, uEltypeNoUnits(c1), uEltypeNoUnits(c2)) +end diff --git a/src/perform_step/dynamical.jl b/src/perform_step/dynamical.jl new file mode 100644 index 000000000..cc925f12a --- /dev/null +++ b/src/perform_step/dynamical.jl @@ -0,0 +1,81 @@ +function verify_f2(f, p, q, pa, t, integrator, ::BAOABConstantCache) + res = f(p, q, pa, t) + res != p && throwex(integrator) +end +function verify_f2(f, res, p, q, pa, t, integrator, ::BAOABCache) + f(res, p, q, pa, t) + res != p && throwex(integrator) +end +function throwex(integrator) + algn = typeof(integrator.alg) + throw(ArgumentError("Algorithm $algn is not applicable if f2(p, q, t) != p")) +end + +function initialize!(integrator, cache::BAOABConstantCache) + @unpack t,dt,uprev,u,p,W = integrator + du1 = integrator.uprev.x[1] + u1 = integrator.uprev.x[2] + + verify_f2(integrator.f.f2, du1, u1, p, t, integrator, cache) + cache.k .= integrator.f.f1(du1,u1,p,t) +end + +function initialize!(integrator, cache::BAOABCache) + @unpack t,dt,uprev,u,p,W = integrator + du1 = integrator.uprev.x[1] + u1 = integrator.uprev.x[2] + + verify_f2(integrator.f.f2, cache.k, du1, u1, p, t, integrator, cache) + integrator.f.f1(cache.k,du1,u1,p,t) +end + +@muladd function perform_step!(integrator,cache::BAOABConstantCache,f=integrator.f) + @unpack t,dt,uprev,u,p,W = integrator + @unpack k, half, c1, c2 = cache + du1 = uprev.x[1] + u1 = uprev.x[2] + + # B + du2 = du1 + half*dt*k + + # A + u2 = u1 + half*dt*du2 + + # O + noise = integrator.g(u2,p,t+dt*half).*W.dW + du3 = c1*du2 + c2*noise + + # A + u = u2 + half*dt*du3 + + # B + k .= f.f1(du3,u,p,t+dt) + du = du3 + half*dt*k + + integrator.u = ArrayPartition((du, u)) +end + +@muladd function perform_step!(integrator,cache::BAOABCache,f=integrator.f) + @unpack t,dt,uprev,u,p,W = integrator + @unpack utmp, dutmp, k, gtmp, noise, half, c1, c2 = cache + du1 = uprev.x[1] + u1 = uprev.x[2] + + # B + @.. dutmp = du1 + half*dt*k + + # A + @.. utmp = u1 + half*dt*dutmp + + # O + integrator.g(gtmp,utmp,p,t+dt*half) + @.. noise = gtmp*W.dW + @.. dutmp = c1*dutmp + c2*noise + + # A + @.. u.x[2] = utmp + half*dt*dutmp + + # B + f.f1(k,dutmp,u.x[2],p,t+dt) + @.. u.x[1] = dutmp + half*dt*k +end diff --git a/src/solve.jl b/src/solve.jl index 77089137e..c59a78753 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -199,7 +199,9 @@ function DiffEqBase.__init( end rateType = typeof(rate_prototype) ## Can be different if united - if is_diagonal_noise(prob) + if prob.f isa DynamicalSDEFunction + noise_rate_prototype = rate_prototype.x[1] + elseif is_diagonal_noise(prob) noise_rate_prototype = rate_prototype elseif prob isa DiffEqBase.AbstractRODEProblem if prob isa DiffEqBase.AbstractSDEProblem @@ -273,7 +275,9 @@ function DiffEqBase.__init( randType = typeof(rand_prototype) else randElType = uBottomEltypeNoUnits # Strip units and type info - if is_diagonal_noise(prob) + if prob.f isa DynamicalSDEFunction + rand_prototype = copy(noise_rate_prototype) + elseif is_diagonal_noise(prob) if typeof(u) <: SArray rand_prototype = zero(u) # TODO: Array{randElType} for units else diff --git a/test/sde/sde_dynamical.jl b/test/sde/sde_dynamical.jl new file mode 100644 index 000000000..cd9af34d0 --- /dev/null +++ b/test/sde/sde_dynamical.jl @@ -0,0 +1,25 @@ + +using StochasticDiffEq, DiffEqNoiseProcess, Test, DiffEqDevTools +using Plots +using RecursiveArrayTools + +u0 = zeros(2) +v0 = ones(2) +γ = 1 +f1_harmonic(v,u,p,t) = -u +f2_harmonic(v,u,p,t) = v +g(u,p,t) = 0.2 + +ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g) +prob = DynamicalSDEProblem(ff_harmonic,g,v0,u0,(0.0,5.0)) + +sol1 = solve(prob,BAOAB(gamma=γ);dt=1/10,save_noise=true) + +f1_harmonic_iip(dv,v,u,p,t) = dv .= f1_harmonic(v,u,p,t) +f2_harmonic_iip(du,v,u,p,t) = du .= f2_harmonic(v,u,p,t) +g_iip(du,u,p,t) = du .= g(u,p,t) + +prob = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W)) + +sol2 = solve(prob,BAOAB(gamma=γ);dt=1/10) +@test sol1[:] ≈ sol2[:] From 277c3937ad780b09677b431201aedf67cc003929 Mon Sep 17 00:00:00 2001 From: James Gardner Date: Mon, 8 Feb 2021 06:47:21 +0000 Subject: [PATCH 2/6] Clean up and add to runtests --- src/alg_utils.jl | 2 +- src/algorithms.jl | 7 ++++++- src/caches/dynamical_caches.jl | 4 ++-- test/runtests.jl | 1 + test/sde/sde_dynamical.jl | 3 --- 5 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 1f1f1f5f0..27885c9d0 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -95,7 +95,7 @@ alg_order(alg::SMEB) = 1//1 alg_order(alg::TauLeaping) = 1//1 alg_order(alg::CaoTauLeaping) = 1//1 -alg_order(alg::BAOAB) = 1//1 # Not sure what order it is +alg_order(alg::BAOAB) = 2//1 alg_order(alg::SKenCarp) = 2//1 alg_order(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = maximum(alg_order.(alg.algs)) diff --git a/src/algorithms.jl b/src/algorithms.jl index ae67f9731..a6349cc2e 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -556,10 +556,15 @@ struct RandomEM <: StochasticDiffEqRODEAlgorithm end const SplitSDEAlgorithms = Union{IIF1M,IIF2M,IIF1Mil,SKenCarp,SplitEM} -""" +@doc raw""" Leimkuhler B., Matthews C., Robust and efficient configurational molecular sampling via Langevin dynamics, J. Chem. Phys. 138, 174102 (2013) DOI:10.1063/1.4802990 + +```math +du = vdt \\ +dv = f(v,u) dt - \gamma v dt + g(u) dW +``` """ struct BAOAB{T} <: StochasticDiffEqAlgorithm gamma::T diff --git a/src/caches/dynamical_caches.jl b/src/caches/dynamical_caches.jl index 5cfd45300..b3357247a 100644 --- a/src/caches/dynamical_caches.jl +++ b/src/caches/dynamical_caches.jl @@ -19,7 +19,7 @@ end function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) k = zero(rate_prototype.x[1]) c1 = exp(-alg.gamma*dt) - c2 = sqrt(1 - c1^2) + c2 = sqrt(1 - c1^2) / sqrt(2alg.gamma) BAOABConstantCache(k, uEltypeNoUnits(1//2), uEltypeNoUnits(c1), uEltypeNoUnits(c2)) end @@ -33,7 +33,7 @@ function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototy half = uEltypeNoUnits(1//2) c1 = exp(-alg.gamma*dt) - c2 = sqrt(1 - c1^2) + c2 = sqrt(1 - c1^2) / sqrt(2alg.gamma) BAOABCache(utmp, dutmp, k, gtmp, noise, half, uEltypeNoUnits(c1), uEltypeNoUnits(c2)) end diff --git a/test/runtests.jl b/test/runtests.jl index 214b65e02..70babc4a5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -57,6 +57,7 @@ const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR") if !is_APPVEYOR && (GROUP == "All" || GROUP == "AlgConvergence") @time @safetestset "Convergence Tests" begin include("sde/sde_convergence_tests.jl") end + @time @safetestset "Dynamical SDE Tests" begin include("sde/sde_dynamical.jl") end end if !is_APPVEYOR && GROUP == "AlgConvergence2" diff --git a/test/sde/sde_dynamical.jl b/test/sde/sde_dynamical.jl index cd9af34d0..6048ac8c9 100644 --- a/test/sde/sde_dynamical.jl +++ b/test/sde/sde_dynamical.jl @@ -1,7 +1,4 @@ - using StochasticDiffEq, DiffEqNoiseProcess, Test, DiffEqDevTools -using Plots -using RecursiveArrayTools u0 = zeros(2) v0 = ones(2) From 3d0e1b2e5c07c0aa4189656977eeb7fcd26178a9 Mon Sep 17 00:00:00 2001 From: James Gardner Date: Tue, 9 Feb 2021 11:42:24 +0000 Subject: [PATCH 3/6] Add convergence tests --- test/sde/sde_dynamical.jl | 29 +++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/test/sde/sde_dynamical.jl b/test/sde/sde_dynamical.jl index 6048ac8c9..a4c145c3e 100644 --- a/test/sde/sde_dynamical.jl +++ b/test/sde/sde_dynamical.jl @@ -1,22 +1,31 @@ -using StochasticDiffEq, DiffEqNoiseProcess, Test, DiffEqDevTools +using StochasticDiffEq, DiffEqNoiseProcess, Test, DiffEqDevTools, Random +Random.seed!(1) u0 = zeros(2) v0 = ones(2) -γ = 1 +γ = 50 + f1_harmonic(v,u,p,t) = -u f2_harmonic(v,u,p,t) = v -g(u,p,t) = 0.2 - -ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g) -prob = DynamicalSDEProblem(ff_harmonic,g,v0,u0,(0.0,5.0)) - -sol1 = solve(prob,BAOAB(gamma=γ);dt=1/10,save_noise=true) +g(u,p,t) = sqrt(2γ) f1_harmonic_iip(dv,v,u,p,t) = dv .= f1_harmonic(v,u,p,t) f2_harmonic_iip(du,v,u,p,t) = du .= f2_harmonic(v,u,p,t) g_iip(du,u,p,t) = du .= g(u,p,t) -prob = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W)) +ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g) +prob1 = DynamicalSDEProblem(ff_harmonic,g,v0,u0,(0.0,5.0)) +sol1 = solve(prob1,BAOAB(gamma=γ);dt=1/10,save_noise=true) + +prob2 = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W)) +sol2 = solve(prob2,BAOAB(gamma=γ);dt=1/10) -sol2 = solve(prob,BAOAB(gamma=γ);dt=1/10) @test sol1[:] ≈ sol2[:] + +dts = (1/2) .^ (9:-1:3) +# Can't use NoiseGrid as noise is not generated with the correct size in convergence.jl. We require noise with shape of v. +# I don't think there are any analytic solutions for problems of this type +sim1 = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false) +@test abs(sim1.𝒪est[:weak_final]-1) < 0.3 # Gives order of only 1? Should be weak order 2. +sim2 = analyticless_test_convergence(dts,prob2,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false) +@test abs(sim2.𝒪est[:weak_final]-1) < 0.3 # Gives order of only 1? Should be weak order 2. From f4c28a9e00739cb128bfd205ff9d858036e56908 Mon Sep 17 00:00:00 2001 From: James Gardner Date: Wed, 24 Feb 2021 10:14:23 +0000 Subject: [PATCH 4/6] Fix mistake in O step --- src/perform_step/dynamical.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/perform_step/dynamical.jl b/src/perform_step/dynamical.jl index cc925f12a..c9526db37 100644 --- a/src/perform_step/dynamical.jl +++ b/src/perform_step/dynamical.jl @@ -30,7 +30,7 @@ function initialize!(integrator, cache::BAOABCache) end @muladd function perform_step!(integrator,cache::BAOABConstantCache,f=integrator.f) - @unpack t,dt,uprev,u,p,W = integrator + @unpack t,dt,sqdt,uprev,u,p,W = integrator @unpack k, half, c1, c2 = cache du1 = uprev.x[1] u1 = uprev.x[2] @@ -42,7 +42,7 @@ end u2 = u1 + half*dt*du2 # O - noise = integrator.g(u2,p,t+dt*half).*W.dW + noise = integrator.g(u2,p,t+dt*half).*W.dW / sqdt du3 = c1*du2 + c2*noise # A @@ -56,7 +56,7 @@ end end @muladd function perform_step!(integrator,cache::BAOABCache,f=integrator.f) - @unpack t,dt,uprev,u,p,W = integrator + @unpack t,dt,sqdt,uprev,u,p,W = integrator @unpack utmp, dutmp, k, gtmp, noise, half, c1, c2 = cache du1 = uprev.x[1] u1 = uprev.x[2] @@ -69,7 +69,7 @@ end # O integrator.g(gtmp,utmp,p,t+dt*half) - @.. noise = gtmp*W.dW + @.. noise = gtmp*W.dW / sqdt @.. dutmp = c1*dutmp + c2*noise # A From 9f0b29c3a62b9afadbe6d0583972488018f552d9 Mon Sep 17 00:00:00 2001 From: James Gardner Date: Wed, 24 Feb 2021 10:15:00 +0000 Subject: [PATCH 5/6] Modify form of equation --- src/algorithms.jl | 2 +- src/caches/dynamical_caches.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index a6349cc2e..d5498c8cc 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -563,7 +563,7 @@ DOI:10.1063/1.4802990 ```math du = vdt \\ -dv = f(v,u) dt - \gamma v dt + g(u) dW +dv = f(v,u) dt - \gamma v dt + g(u) \sqrt{2\gamma} dW ``` """ struct BAOAB{T} <: StochasticDiffEqAlgorithm diff --git a/src/caches/dynamical_caches.jl b/src/caches/dynamical_caches.jl index b3357247a..5cfd45300 100644 --- a/src/caches/dynamical_caches.jl +++ b/src/caches/dynamical_caches.jl @@ -19,7 +19,7 @@ end function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) k = zero(rate_prototype.x[1]) c1 = exp(-alg.gamma*dt) - c2 = sqrt(1 - c1^2) / sqrt(2alg.gamma) + c2 = sqrt(1 - c1^2) BAOABConstantCache(k, uEltypeNoUnits(1//2), uEltypeNoUnits(c1), uEltypeNoUnits(c2)) end @@ -33,7 +33,7 @@ function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototy half = uEltypeNoUnits(1//2) c1 = exp(-alg.gamma*dt) - c2 = sqrt(1 - c1^2) / sqrt(2alg.gamma) + c2 = sqrt(1 - c1^2) BAOABCache(utmp, dutmp, k, gtmp, noise, half, uEltypeNoUnits(c1), uEltypeNoUnits(c2)) end From 28446aa9e6d8bae80ac52ff847de2b4300c0cb86 Mon Sep 17 00:00:00 2001 From: James Gardner Date: Wed, 24 Feb 2021 10:16:02 +0000 Subject: [PATCH 6/6] Fix order and testing Seems like the algorithm is actually just 1st order. Order estimates both strong and weak are 1, even with small dt and many trajectories. --- src/alg_utils.jl | 2 +- test/sde/sde_dynamical.jl | 12 +++++------- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 27885c9d0..2cd7b260f 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -95,7 +95,7 @@ alg_order(alg::SMEB) = 1//1 alg_order(alg::TauLeaping) = 1//1 alg_order(alg::CaoTauLeaping) = 1//1 -alg_order(alg::BAOAB) = 2//1 +alg_order(alg::BAOAB) = 1//1 alg_order(alg::SKenCarp) = 2//1 alg_order(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = maximum(alg_order.(alg.algs)) diff --git a/test/sde/sde_dynamical.jl b/test/sde/sde_dynamical.jl index a4c145c3e..3c49159b6 100644 --- a/test/sde/sde_dynamical.jl +++ b/test/sde/sde_dynamical.jl @@ -3,11 +3,11 @@ Random.seed!(1) u0 = zeros(2) v0 = ones(2) -γ = 50 +γ = 1 f1_harmonic(v,u,p,t) = -u f2_harmonic(v,u,p,t) = v -g(u,p,t) = sqrt(2γ) +g(u,p,t) = 1 f1_harmonic_iip(dv,v,u,p,t) = dv .= f1_harmonic(v,u,p,t) f2_harmonic_iip(du,v,u,p,t) = du .= f2_harmonic(v,u,p,t) @@ -22,10 +22,8 @@ sol2 = solve(prob2,BAOAB(gamma=γ);dt=1/10) @test sol1[:] ≈ sol2[:] -dts = (1/2) .^ (9:-1:3) +dts = (1/2) .^ (8:-1:4) + # Can't use NoiseGrid as noise is not generated with the correct size in convergence.jl. We require noise with shape of v. -# I don't think there are any analytic solutions for problems of this type sim1 = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false) -@test abs(sim1.𝒪est[:weak_final]-1) < 0.3 # Gives order of only 1? Should be weak order 2. -sim2 = analyticless_test_convergence(dts,prob2,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false) -@test abs(sim2.𝒪est[:weak_final]-1) < 0.3 # Gives order of only 1? Should be weak order 2. +@test abs(sim1.𝒪est[:weak_final]-1) < 0.3