diff --git a/src/caches/dynamical_caches.jl b/src/caches/dynamical_caches.jl index 22dae0cec..5980526fc 100644 --- a/src/caches/dynamical_caches.jl +++ b/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) diff --git a/src/perform_step/dynamical.jl b/src/perform_step/dynamical.jl index 074cbc7dd..68c2d1fa5 100644 --- a/src/perform_step/dynamical.jl +++ b/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 diff --git a/test/sde/sde_dynamical.jl b/test/sde/sde_dynamical.jl index a0a259810..8e4475f04 100644 --- a/test/sde/sde_dynamical.jl +++ b/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