expanded noise Euler-Heun

Author: ChrisRackauckas

src/alg_utils.jl

@@ -26,3 +26,4 @@ alg_compatible(prob,alg::SRI) = is_diagonal_noise(prob)
 alg_compatible(prob,alg::SRIW1) = is_diagonal_noise(prob)
 alg_compatible(prob,alg::SRA) = is_diagonal_noise(prob)
 alg_compatible(prob,alg::SRA1) = is_diagonal_noise(prob)
+alg_compatible(prob,alg::RKMil) = is_diagonal_noise(prob)

src/caches.jl

@@ -40,14 +40,15 @@ function alg_cache(alg::EM,prob,u,ΔW,ΔZ,rate_prototype,noise_rate_prototype,uE
 end
 
 immutable EulerHeunConstantCache <: StochasticDiffEqConstantCache end
-immutable EulerHeunCache{uType,rateType} <: StochasticDiffEqMutableCache
+immutable EulerHeunCache{uType,rateType,rateNoiseType,rateNoiseCollectionType} <: StochasticDiffEqMutableCache
   u::uType
   uprev::uType
   tmp::uType
-  rtmp1::rateType
-  rtmp2::rateType
-  rtmp3::rateType
-  rtmp4::rateType
+  ftmp1::rateType
+  ftmp2::rateType
+  nrtmp::rateNoiseCollectionType
+  gtmp1::rateNoiseType
+  gtmp2::rateNoiseType
 end
 
 u_cache(c::EulerHeunCache) = ()
@@ -56,9 +57,10 @@ du_cache(c::EulerHeunCache) = (c.rtmp1,c.rtmp2,c.rtmp3,c.rtmp4)
 alg_cache(alg::EulerHeun,prob,u,ΔW,ΔZ,rate_prototype,noise_rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{false}}) = EulerHeunConstantCache()
 
 function alg_cache(alg::EulerHeun,prob,u,ΔW,ΔZ,rate_prototype,noise_rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{true}})
-  tmp = similar(u); rtmp1 = zeros(rate_prototype); rtmp2 = zeros(rate_prototype)
-  rtmp3 = zeros(rate_prototype); rtmp4 = zeros(rate_prototype)
-  EulerHeunCache(u,uprev,tmp,rtmp1,rtmp2,rtmp3,rtmp4)
+  tmp = similar(u); ftmp1 = zeros(rate_prototype); ftmp2 = zeros(rate_prototype)
+  nrtmp = zeros(rate_prototype)
+  gtmp1 = zeros(noise_rate_prototype); gtmp2 = zeros(noise_rate_prototype)
+  EulerHeunCache(u,uprev,tmp,ftmp1,ftmp2,nrtmp,gtmp1,gtmp2)
 end
 
 immutable RandomEMConstantCache <: StochasticDiffEqConstantCache end

src/integrators/low_order.jl

@@ -34,19 +34,41 @@ end
 end
 
 @inline function perform_step!(integrator,cache::EulerHeunCache,f=integrator.f)
-  @unpack rtmp1,rtmp2,rtmp3,rtmp4,tmp = cache
+  @unpack ftmp1,ftmp2,gtmp1,gtmp2,tmp,nrtmp = cache
   @unpack t,dt,uprev,u,ΔW = integrator
-  integrator.f(t,uprev,rtmp1)
-  integrator.g(t,uprev,rtmp2)
+  integrator.f(t,uprev,ftmp1)
+  integrator.g(t,uprev,gtmp1)
+
+  if is_diagonal_noise(integrator.sol.prob)
+    for i in eachindex(u)
+      nrtmp[i]=gtmp1[i]*ΔW[i]
+    end
+  else
+    A_mul_B!(nrtmp,gtmp1,ΔW)
+  end
+
   for i in eachindex(u)
-    tmp[i] = @muladd uprev[i] + rtmp1[i]*dt + rtmp2[i]*ΔW[i]
+    tmp[i] = @muladd uprev[i] + ftmp1[i]*dt + nrtmp[i]
   end
-  integrator.f(t,tmp,rtmp3)
-  integrator.g(t,tmp,rtmp4)
+
+  integrator.f(t,tmp,ftmp2)
+  integrator.g(t,tmp,gtmp2)
+
+  if is_diagonal_noise(integrator.sol.prob)
+    for i in eachindex(u)
+      ΔWo2 = (1/2)*ΔW[i]
+      nrtmp[i]=ΔWo2*(gtmp1[i]+gtmp2[i])
+    end
+  else
+    for i in eachindex(gtmp1)
+      gtmp1[i] = (1/2)*(gtmp1[i]+gtmp2[i])
+    end
+    A_mul_B!(nrtmp,gtmp1,ΔW)
+  end
+
   dto2 = dt*(1/2)
   for i in eachindex(u)
-    ΔWo2 = (1/2)*ΔW[i]
-    u[i] = @muladd uprev[i] + dto2*(rtmp1[i]+rtmp3[i]) + ΔWo2*(rtmp2[i]+rtmp4[i])
+    u[i] = @muladd uprev[i] + dto2*(ftmp1[i]+ftmp2[i]) + nrtmp[i]
   end
   @pack integrator = t,dt,u
 end

test/stratonovich_convergence_tests.jl

@@ -1,4 +1,4 @@
-@everywhere using StochasticDiffEq, DiffEqProblemLibrary, DiffEqDevTools, Base.Test
+using StochasticDiffEq, DiffEqProblemLibrary, DiffEqDevTools, Base.Test
 srand(100)
 dts = 1./2.^(10:-1:2) #14->7 good plot
 
@@ -19,3 +19,6 @@ sim  = test_convergence(dts,prob,EulerHeun(),numMonte=Int(5e1))
 sim2 = test_convergence(dts,prob,RKMil(interpretation=:Stratonovich),numMonte=Int(5e2))
 @test abs(sim2.𝒪est[:l∞]-1) < 0.1
 =#
+
+srand(200)
+sol = solve(prob,EulerHeun(),dt=1/4)