Merge pull request #559 from Lightup1/baoab-dt

fix dt scaling when scale_noise==false

Author: ChrisRackauckas

src/caches/dynamical_caches.jl

@@ -20,7 +20,7 @@ end
 function alg_cache(alg::BAOAB,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}
   k = zero(rate_prototype.x[1])
   c1 = exp(-alg.gamma*dt)
-  c2 = sqrt(1 - alg.scale_noise*c1^2) # if scale_noise == false, c2 = 1
+  c2 = alg.scale_noise ? sqrt((1 - c1^2)/abs(dt)) : 1 # if scale_noise == false, c2 = 1
   BAOABConstantCache(k, uEltypeNoUnits(1//2), uEltypeNoUnits(c1), uEltypeNoUnits(c2))
 end
 
@@ -34,7 +34,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 - alg.scale_noise*c1^2) # if scale_noise == false, c2 = 1
+  c2 = alg.scale_noise ? sqrt((1 - c1^2)/abs(dt)) : 1 # if scale_noise == false, c2 = 1
 
   tmp = zero(u)
 

src/perform_step/dynamical.jl

@@ -42,7 +42,7 @@ end
   u2 = u1 + half*dt*du2
 
   # O
-  noise = integrator.g(u2,p,t+dt*half).*W.dW / sqdt
+  noise = integrator.g(u2,p,t+dt*half).*W.dW
   du3 = c1*du2 + c2*noise
 
   # A
@@ -69,7 +69,7 @@ end
 
   # O
   integrator.g(gtmp,utmp,p,t+dt*half)
-  @.. noise = gtmp*W.dW / sqdt
+  @.. noise = gtmp*W.dW
   @.. dutmp = c1*dutmp + c2*noise
 
   # A

test/sde/sde_dynamical.jl

@@ -45,3 +45,43 @@ end
     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
 end
+
+@testset "Scalar u, scale_noise=false" begin
+    u0 = 0
+    v0 = 1
+
+    ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
+    prob1 = DynamicalSDEProblem(ff_harmonic,v0,u0,(0.0,5.0))
+
+    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.
+    sim1  = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ,scale_noise=false),(1/2)^10;trajectories=Int(2e2),use_noise_grid=false)
+    display(sim1.𝒪est)
+    @test abs(sim1.𝒪est[:weak_final]-1) < 0.3
+end
+
+@testset "Vector u, scale_noise=false" begin
+
+    u0 = zeros(2)
+    v0 = ones(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)
+
+    ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
+    prob1 = DynamicalSDEProblem(ff_harmonic,v0,u0,(0.0,5.0))
+    sol1 = solve(prob1,BAOAB(gamma=γ,scale_noise=false);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=γ,scale_noise=false);dt=1/10)
+
+    @test sol1[:] ≈ sol2[:]
+
+    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.
+    sim1  = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ,scale_noise=false),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false)
+    @test abs(sim1.𝒪est[:weak_final]-1) < 0.3
+end