Merge pull request #501 from rmsrosa/random_tamed_em

add RandomTamedEM() solver for RODEs

Author: ChrisRackauckas

src/StochasticDiffEq.jl

@@ -165,9 +165,7 @@ using DocStringExtensions
   export StochasticDiffEqRODEAlgorithm, StochasticDiffEqRODEAdaptiveAlgorithm,
          StochasticDiffEqRODECompositeAlgorithm
 
-  export RandomEM
-  
-  export RandomHeun
+  export RandomEM, RandomTamedEM, RandomHeun
 
   export IteratedIntegralApprox, IICommutative, IILevyArea
 

src/alg_utils.jl

@@ -60,6 +60,7 @@ alg_order(alg::EulerHeun) = 1 // 2
 alg_order(alg::LambaEulerHeun) = 1 // 2
 alg_order(alg::RandomEM) = 1 // 2
 alg_order(alg::RandomHeun) = 1 // 2
+alg_order(alg::RandomTamedEM) = 1 // 2
 alg_order(alg::SimplifiedEM) = 1 // 2
 alg_order(alg::RKMil) = 1 // 1
 alg_order(alg::RKMilCommute) = 1 // 1

src/algorithms.jl

@@ -850,6 +850,8 @@ struct RandomEM <: StochasticDiffEqRODEAlgorithm end
 
 struct RandomHeun <: StochasticDiffEqRODEAlgorithm end
 
+struct RandomTamedEM <: StochasticDiffEqRODEAlgorithm end
+
 const SplitSDEAlgorithms = Union{IIF1M,IIF2M,IIF1Mil,SKenCarp,SplitEM}
 
 @doc raw"""

src/caches/basic_method_caches.jl

@@ -72,6 +72,22 @@ function alg_cache(alg::RandomEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prot
   RandomEMCache(u,uprev,tmp,rtmp)
 end
 
+struct RandomTamedEMConstantCache <: StochasticDiffEqConstantCache end
+
+@cache struct RandomTamedEMCache{uType,rateType} <: StochasticDiffEqMutableCache
+  u::uType
+  uprev::uType
+  tmp::uType
+  rtmp::rateType
+end
+
+alg_cache(alg::RandomTamedEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = RandomTamedEMConstantCache()
+
+function alg_cache(alg::RandomTamedEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  tmp = zero(u); rtmp = zero(rate_prototype)
+  RandomTamedEMCache(u,uprev,tmp,rtmp)
+end
+
 struct RandomHeunConstantCache <: StochasticDiffEqConstantCache end
 @cache struct RandomHeunCache{uType,rateType,randType} <: StochasticDiffEqMutableCache
   u::uType

src/perform_step/low_order.jl

@@ -121,6 +121,20 @@ end
   @.. u = uprev + dt * rtmp
 end
 
+@muladd function perform_step!(integrator,cache::RandomTamedEMConstantCache)
+    @unpack t,dt,uprev,u,W,p,f = integrator
+    ftmp = integrator.f(uprev,p,t,W.curW)
+    u = uprev .+ dt .* ftmp ./ (1 .+ dt .* norm(ftmp))
+    integrator.u = u
+end
+  
+@muladd function perform_step!(integrator,cache::RandomTamedEMCache)
+    @unpack rtmp = cache
+    @unpack t,dt,uprev,u,W,p,f = integrator
+    integrator.f(rtmp,uprev,p,t,W.curW)
+    @.. u = uprev + dt * rtmp / (1 + dt * norm(rtmp))
+end
+
 @muladd function perform_step!(integrator,cache::RandomHeunConstantCache)
     @unpack t,dt,uprev,u,W,p,f = integrator
     ftmp = integrator.f(uprev,p,t,W.curW)

test/rode_linear_tests.jl

@@ -8,8 +8,9 @@ prob = RODEProblem(f,u0,tspan)
 sol = solve(prob,RandomEM(),dt=1/100, save_noise=true)
 prob = RODEProblem(f,u0,tspan, noise=NoiseWrapper(sol.W))
 sol2 = solve(prob,RandomHeun(),dt=1/100)
-@test abs(sol[end] - sol2[end]) < 0.1
-
+@test abs(sol[end] - sol2[end]) < 0.1 * abs(sol[end])
+sol3 = solve(prob,RandomTamedEM(),dt=1/100)
+@test abs(sol[end] - sol3[end]) < 0.1 * abs(sol[end])
 
 f(du,u,p,t,W) = (du.=1.01u.+0.87u.*W)
 u0 = ones(4)
@@ -18,36 +19,44 @@ sol = solve(prob,RandomEM(),dt=1/100, save_noise=true)
 prob = RODEProblem(f,u0,tspan, noise=NoiseWrapper(sol.W))
 sol2 = solve(prob,RandomHeun(),dt=1/100)
 @test sum(abs,sol[end]-sol2[end]) < 0.1 * sum(abs, sol[end])
+sol3 = solve(prob,RandomEM(),dt=1/100)
+@test sum(abs,sol[end]-sol3[end]) < 0.1 * sum(abs, sol[end])
 
 f(u,p,t,W) = 2u*sin(W)
 u0 = 1.00
-tspan = (0.0,5.0)
+tspan = (0.0,1.0)
 prob = RODEProblem{false}(f,u0,tspan)
 sol = solve(prob,RandomEM(),dt=1/100, save_noise=true)
 prob = RODEProblem{false}(f,u0,tspan, noise=NoiseWrapper(sol.W))
 sol2 = solve(prob,RandomHeun(),dt=1/100)
-@test abs(sol[end]-sol2[end]) < 0.1
+@test abs(sol[end]-sol2[end]) < 0.1 * abs(sol[end])
+sol3 = solve(prob,RandomTamedEM(),dt=1/100)
+@test abs(sol[end]-sol3[end]) < 0.1 * abs(sol[end])
 
 function f(du,u,p,t,W)
-  du[1] = 2u[1]*sin(W[1] - W[2])
-  du[2] = -2u[2]*cos(W[1] + W[2])
+  du[1] = 0.2u[1]*sin(W[1] - W[2])
+  du[2] = -0.2u[2]*cos(W[1] + W[2])
 end
 u0 = [1.00;1.00]
-tspan = (0.0,4.0)
+tspan = (0.0,1.0)
 prob = RODEProblem(f,u0,tspan)
 sol = solve(prob,RandomEM(),dt=1/100, save_noise=true)
 prob = RODEProblem(f,u0,tspan, noise=NoiseWrapper(sol.W))
 sol2 = solve(prob,RandomHeun(),dt=1/100)
 @test sum(abs,sol[end]-sol2[end]) < 0.1 * sum(abs, sol[end])
+sol3 = solve(prob,RandomTamedEM(),dt=1/100)
+@test sum(abs,sol[end]-sol3[end]) < 0.1 * sum(abs, sol[end])
 
 function f(du,u,p,t,W)
-  du[1] = -2W[3]*u[1]*sin(W[1] - W[2])
-  du[2] = -2u[2]*cos(W[1] + W[2])
+  du[1] = -0.2W[3]*u[1]*sin(W[1] - W[2])
+  du[2] = -0.2u[2]*cos(W[1] + W[2])
 end
 u0 = [1.00;1.00]
-tspan = (0.0,5.0)
+tspan = (0.0,1.0)
 prob = RODEProblem(f,u0,tspan,rand_prototype=zeros(3))
 sol = solve(prob,RandomEM(),dt=1/100, save_noise=true)
 prob = RODEProblem(f,u0,tspan, noise=NoiseWrapper(sol.W))
 sol2 = solve(prob,RandomHeun(),dt=1/100)
 @test sum(abs,sol[end]-sol2[end]) < 0.1 * sum(abs, sol[end])
+sol3 = solve(prob,RandomTamedEM(),dt=1/100)
+@test sum(abs,sol[end]-sol3[end]) < 0.1 * sum(abs, sol3[end])