add extensive scalar noise test set for different solvers

Author: frankschae

Project.toml

@@ -60,11 +60,12 @@ DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 DiffEqOperators = "9fdde737-9c7f-55bf-ade8-46b3f136cc48"
 DiffEqProblemLibrary = "a077e3f3-b75c-5d7f-a0c6-6bc4c8ec64a9"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DiffEqCallbacks", "DiffEqDevTools", "DiffEqOperators", "DiffEqProblemLibrary", "LinearSolve", "Pkg", "SafeTestsets", "SparseArrays", "Statistics", "Test"]
+test = ["DiffEqCallbacks", "DiffEqDevTools", "DiffEqOperators", "DiffEqProblemLibrary", "LinearSolve", "ModelingToolkit", "Pkg", "SafeTestsets", "SparseArrays", "Statistics", "Test"]

src/perform_step/SROCK_perform_step.jl

@@ -3,7 +3,7 @@
 
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
-  cache.mdeg = Int(floor(sqrt(2*dt*integrator.opts.internalnorm(integrator.eigen_est,t))+1)) # this is the spectral radius estimate to choose optimal stage
+  cache.mdeg = Int(floor(sqrt(2*abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t))+1)) # this is the spectral radius estimate to choose optimal stage
   choose_deg!(integrator,cache)
 
   mdeg = cache.mdeg
@@ -61,14 +61,14 @@
     elseif (i == mdeg) && alg_interpretation(integrator.alg) == :Ito
       if typeof(W.dW) <: Number
         gₘ₋₂ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
-        uᵢ₋₂ = uᵢ₋₁ + sqrt(dt)*gₘ₋₂
+        uᵢ₋₂ = uᵢ₋₁ + sqrt(abs(dt))*gₘ₋₂
         gₘ₋₁ = integrator.g(uᵢ₋₂,p,tᵢ₋₁)
-        u += gₘ₋₂*W.dW + 1/(2*sqrt(dt))*(gₘ₋₁ - gₘ₋₂)*(W.dW^2 - dt)
+        u += gₘ₋₂*W.dW + 1/(2*sqrt(abs(dt)))*(gₘ₋₁ - gₘ₋₂)*(W.dW^2 - abs(dt))
       elseif is_diagonal_noise(integrator.sol.prob)
         gₘ₋₂ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
-        uᵢ₋₂ .= uᵢ₋₁ .+ sqrt(dt) .* gₘ₋₂
+        uᵢ₋₂ .= uᵢ₋₁ .+ sqrt(abs(dt)) .* gₘ₋₂
         gₘ₋₁ = integrator.g(uᵢ₋₂,p,tᵢ₋₁)
-        u .+= gₘ₋₂ .* W.dW .+ (1/(2*sqrt(dt))) .* (gₘ₋₁ .- gₘ₋₂) .* (W.dW .^ 2 .- dt)
+        u .+= gₘ₋₂ .* W.dW .+ (1/(2*sqrt(abs(dt)))) .* (gₘ₋₁ .- gₘ₋₂) .* (W.dW .^ 2 .- abs(dt))
       else
         gₘ₋₂ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
         u += gₘ₋₂*W.dW
@@ -94,7 +94,7 @@ end
   ccache = cache.constantcache
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
-  ccache.mdeg = Int(floor(sqrt(2*dt*integrator.opts.internalnorm(integrator.eigen_est,t))+1))   # this is the spectral radius estimate to choose optimal stage
+  ccache.mdeg = Int(floor(sqrt(2*abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t))+1))   # this is the spectral radius estimate to choose optimal stage
   choose_deg!(integrator,cache)
 
   mdeg = ccache.mdeg
@@ -155,9 +155,9 @@ end
     elseif (i == mdeg) && alg_interpretation(integrator.alg) == :Ito
       if typeof(W.dW) <: Number || is_diagonal_noise(integrator.sol.prob)
         integrator.g(gₘ₋₂,uᵢ₋₁,p,tᵢ₋₁)
-        @.. uᵢ₋₂ = uᵢ₋₁ + sqrt(dt)*gₘ₋₂
+        @.. uᵢ₋₂ = uᵢ₋₁ + sqrt(abs(dt))*gₘ₋₂
         integrator.g(gₘ₋₁,uᵢ₋₂,p,tᵢ₋₁)
-        @.. u += gₘ₋₂*W.dW + 1/(2*sqrt(dt))*(gₘ₋₁ - gₘ₋₂)*(W.dW^2 - dt)
+        @.. u += gₘ₋₂*W.dW + 1/(2*sqrt(abs(dt)))*(gₘ₋₁ - gₘ₋₂)*(W.dW^2 - abs(dt))
       else
         integrator.g(gₘ₋₂,uᵢ₋₁,p,tᵢ₋₁)
         mul!(uᵢ₋₁,gₘ₋₂,W.dW)
@@ -187,7 +187,7 @@ end
 
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
-  cache.mdeg = Int(floor(sqrt((2*dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
+  cache.mdeg = Int(floor(sqrt((2*abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
   cache.mdeg = max(3,min(cache.mdeg,200))-2
   choose_deg!(integrator,cache)
 
@@ -319,7 +319,7 @@ end
 
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
-  ccache.mdeg = Int(floor(sqrt((2*dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
+  ccache.mdeg = Int(floor(sqrt((2*abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
   ccache.mdeg = max(3,min(ccache.mdeg,200))-2
   choose_deg!(integrator,cache)
 
@@ -461,9 +461,9 @@ end
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
   if integrator.alg.strong_order_1
-    cache.mdeg = Int(floor(sqrt(dt*integrator.opts.internalnorm(integrator.eigen_est,t)/0.19)+1))
+    cache.mdeg = Int(floor(sqrt(abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)/0.19)+1))
   else
-    cache.mdeg = Int(floor(sqrt(dt*integrator.opts.internalnorm(integrator.eigen_est,t)/0.33)+1))
+    cache.mdeg = Int(floor(sqrt(abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)/0.33)+1))
   end
   cache.mdeg = max(3,min(cache.mdeg,200))
   choose_deg!(integrator,cache)
@@ -552,9 +552,9 @@ end
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
   if integrator.alg.strong_order_1
-    ccache.mdeg = Int(floor(sqrt(dt*integrator.opts.internalnorm(integrator.eigen_est,t)/0.19)+1))
+    ccache.mdeg = Int(floor(sqrt(abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)/0.19)+1))
   else
-    ccache.mdeg = Int(floor(sqrt(dt*integrator.opts.internalnorm(integrator.eigen_est,t)/0.33)+1))
+    ccache.mdeg = Int(floor(sqrt(abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)/0.33)+1))
   end
   ccache.mdeg = max(3,min(ccache.mdeg,200))
   choose_deg!(integrator,cache)
@@ -643,7 +643,7 @@ end
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
   η = oftype(t,0.05)
-  mdeg = Int(floor(sqrt((dt*integrator.opts.internalnorm(integrator.eigen_est,t) + 1.5)/(2-η*4/3))+1))
+  mdeg = Int(floor(sqrt((abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t) + 1.5)/(2-η*4/3))+1))
   mdeg = max(3,min(mdeg,200))
 
   ω₀ = 1 + (η/(mdeg^2))
@@ -729,7 +729,7 @@ end
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
   η = oftype(t,0.05)
-  mdeg = Int(floor(sqrt((dt*integrator.opts.internalnorm(integrator.eigen_est,t) + 1.5)/(2-η*4/3))+1))
+  mdeg = Int(floor(sqrt((abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t) + 1.5)/(2-η*4/3))+1))
   mdeg = max(3,min(mdeg,200))
 
   ω₀ = 1 + (η/(mdeg^2))
@@ -819,8 +819,8 @@ end
 
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
-  (integrator.alg.version_num <= 2) && (cache.mdeg = Int(floor(sqrt((dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1)))
-  (integrator.alg.version_num > 2) && (cache.mdeg = Int(floor(sqrt((dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.611)+1)))
+  (integrator.alg.version_num <= 2) && (cache.mdeg = Int(floor(sqrt((abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1)))
+  (integrator.alg.version_num > 2) && (cache.mdeg = Int(floor(sqrt((abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.611)+1)))
 
   cache.mdeg = max(4,min(cache.mdeg,200))-2
   choose_deg!(integrator,cache)
@@ -974,8 +974,8 @@ end
 
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
-  (integrator.alg.version_num <= 2) && (ccache.mdeg = Int(floor(sqrt((dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1)))
-  (integrator.alg.version_num > 2) && (ccache.mdeg = Int(floor(sqrt((dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.611)+1)))
+  (integrator.alg.version_num <= 2) && (ccache.mdeg = Int(floor(sqrt((abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1)))
+  (integrator.alg.version_num > 2) && (ccache.mdeg = Int(floor(sqrt((abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.611)+1)))
   ccache.mdeg = max(4,min(ccache.mdeg,200))-2
   choose_deg!(integrator,cache)
 
@@ -1109,7 +1109,7 @@ end
 
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
-  cache.mdeg = Int(floor(sqrt((2*dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
+  cache.mdeg = Int(floor(sqrt((2*abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
   cache.mdeg = max(6,min(cache.mdeg,200))-2
   choose_deg!(integrator,cache)
 
@@ -1298,7 +1298,7 @@ end
 
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
-  ccache.mdeg = Int(floor(sqrt((2*dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
+  ccache.mdeg = Int(floor(sqrt((2*abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
   ccache.mdeg = max(6,min(ccache.mdeg,200))-2
   choose_deg!(integrator,cache)
 
@@ -1511,7 +1511,7 @@ end
 
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
-  cache.mdeg = Int(floor(sqrt((2*dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
+  cache.mdeg = Int(floor(sqrt((2*abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
   cache.mdeg = max(3,min(cache.mdeg,200))-2
   choose_deg!(integrator,cache)
 
@@ -1597,7 +1597,7 @@ end
 
   alg = unwrap_alg(integrator, true)
   alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
-  ccache.mdeg = Int(floor(sqrt((2*dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
+  ccache.mdeg = Int(floor(sqrt((2*abs(dt)*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
   ccache.mdeg = max(3,min(ccache.mdeg,200))-2
   choose_deg!(integrator,cache)
 

src/perform_step/explicit_3s_mil_methods.jl

@@ -16,7 +16,7 @@
   k  = integrator.g(u,p,t)
   tmp = u + k*integrator.sqdt
   k₁ = integrator.g(tmp,p,t)
-  k₁ = (k₁ - k)/(integrator.sqdt)
+  k₁ = (k₁ - k)/abs(integrator.sqdt)
   u  = u .+ k .* W.dW .+ 0.5 .* (W.dW .^ 2) .* k₁
   integrator.u = u
 end
@@ -40,7 +40,7 @@ end
   integrator.g(k,u,p,t)
   @.. tmp = u + k*integrator.sqdt
   integrator.g(k₁,tmp,p,t)
-  @.. k₁ = (k₁ - k)/(integrator.sqdt)
+  @.. k₁ = (k₁ - k)/abs(integrator.sqdt)
   @.. u  = u + k*W.dW + (W.dW^2/2)*k₁
   integrator.u = u
 end
@@ -62,7 +62,7 @@ end
   k   = integrator.g(u,p,t)
   tmp = u + k*integrator.sqdt
   k₁  = integrator.g(tmp,p,t)
-  k₁  = (k₁ - k)/integrator.sqdt
+  k₁  = (k₁ - k)/abs(integrator.sqdt)
   u   = u - (dt/2)*k₁
   integrator.u = u
 end
@@ -86,7 +86,7 @@ end
   integrator.g(k,u,p,t)
   @.. tmp = u + k*integrator.sqdt
   integrator.g(k₁,tmp,p,t)
-  @.. k₁  = (k₁ - k)/integrator.sqdt
+  @.. k₁  = (k₁ - k)/abs(integrator.sqdt)
   @.. u   = u - (dt/2)*k₁
 
   integrator.u = u
@@ -109,7 +109,7 @@ end
   k  = integrator.g(u,p,t)
   tmp = u + k*integrator.sqdt
   k₁ = integrator.g(tmp,p,t)
-  k₁ = (k₁ - k)/(integrator.sqdt)
+  k₁ = (k₁ - k)/abs(integrator.sqdt)
   u  = u - (dt/2)*k₁
   integrator.u = u
 end
@@ -133,7 +133,7 @@ end
   integrator.g(k,u,p,t)
   @.. tmp = u + k*integrator.sqdt
   integrator.g(k₁,tmp,p,t)
-  @.. k₁ = (k₁ - k)/(integrator.sqdt)
+  @.. k₁ = (k₁ - k)/abs(integrator.sqdt)
   @.. u  = u - (dt/2)*k₁
   integrator.u = u
 end
@@ -155,7 +155,7 @@ end
   k  = integrator.g(u,p,t)
   tmp = u + k*integrator.sqdt
   k₁ = integrator.g(tmp,p,t)
-  k₁ = (k₁ - k)/(integrator.sqdt)
+  k₁ = (k₁ - k)/abs(integrator.sqdt)
   u  = u .+ k .* W.dW .+ 0.5 .* (W.dW .^ 2) .* k₁
  integrator.u = u
 end
@@ -179,7 +179,7 @@ end
   integrator.g(k,u,p,t)
   @.. tmp = u + k*integrator.sqdt
   integrator.g(k₁,tmp,p,t)
-  @.. k₁ = (k₁ - k)/(integrator.sqdt)
+  @.. k₁ = (k₁ - k)/abs(integrator.sqdt)
   @.. u  = u + k*W.dW + (W.dW^2/2)*k₁
   integrator.u = u
 end
@@ -197,7 +197,7 @@ end
   k  = integrator.g(u,p,t)
   tmp = u + k*integrator.sqdt
   k₁ = integrator.g(tmp,p,t)
-  k₁ = (k₁ - k)/(integrator.sqdt)
+  k₁ = (k₁ - k)/abs(integrator.sqdt)
   u  = u .+ k .* W.dW .+ 0.5 .* (W.dW .^ 2) .* k₁
 
   #stage 3
@@ -221,7 +221,7 @@ end
   integrator.g(k,u,p,t)
   @.. tmp = u + k*integrator.sqdt
   integrator.g(k₁,tmp,p,t)
-  @.. k₁ = (k₁ - k)/(integrator.sqdt)
+  @.. k₁ = (k₁ - k)/abs(integrator.sqdt)
   @.. u  = u + k*W.dW + (W.dW^2/2)*k₁
 
   #stage 3
@@ -243,7 +243,7 @@ end
   k   = integrator.g(u,p,t)
   tmp = u + k*integrator.sqdt
   k₁  = integrator.g(tmp,p,t)
-  k₁  = (k₁ - k)/integrator.sqdt
+  k₁  = (k₁ - k)/abs(integrator.sqdt)
   u   = u - (dt/2)*k₁
 
   #stage 3
@@ -267,7 +267,7 @@ end
   integrator.g(k,u,p,t)
   @.. tmp = u + k*integrator.sqdt
   integrator.g(k₁,tmp,p,t)
-  @.. k₁  = (k₁ - k)/integrator.sqdt
+  @.. k₁  = (k₁ - k)/abs(integrator.sqdt)
   @.. u   = u - (dt/2)*k₁
 
   #stage 3

src/perform_step/low_order.jl

@@ -240,9 +240,9 @@ end
   mil_correction = zero(u)
   if alg_interpretation(integrator.alg) == :Ito
     if typeof(dW) <: Number || is_diagonal_noise(integrator.sol.prob)
-      J = J .- 1//2 .* dt
+      J = J .- 1//2 .* abs(dt)
     else
-      J -= 1//2 .* UniformScaling(dt)
+      J -= 1//2 .* UniformScaling(abs(dt))
     end
   end
 
@@ -308,9 +308,9 @@ end
   @.. mil_correction = zero(u)
   if alg_interpretation(integrator.alg) == :Ito
     if typeof(dW) <: Number || is_diagonal_noise(integrator.sol.prob)
-      @.. J -= 1//2*dt
+      @.. J -= 1 // 2 * abs(dt)
     else
-      J -= 1//2 .* UniformScaling(dt)
+      J -= 1//2 .* UniformScaling(abs(dt))
     end
   end
 
@@ -362,9 +362,9 @@ end
 
   if alg_interpretation(integrator.alg) == :Ito
     if typeof(dW) <: Number || is_diagonal_noise(integrator.sol.prob)
-      J = J .- 1//2 .* dt
+      J = J .- 1//2 .* abs(dt)
     else
-      J -= 1//2 .* UniformScaling(dt)
+      J -= 1//2 .* UniformScaling(abs(dt))
     end
   end
 
@@ -432,9 +432,9 @@ end
 
   if alg_interpretation(integrator.alg) == :Ito
     if typeof(dW) <: Number || is_diagonal_noise(integrator.sol.prob)
-      @.. J -= 1//2*dt
+      @.. J -= 1 // 2 * abs(dt)
     else
-      J -= 1//2 .* UniformScaling(dt)
+      J -= 1 // 2 .* UniformScaling(abs(dt))
     end
   end
 

src/perform_step/sdirk.jl

@@ -29,7 +29,7 @@
       K = @.. uprev + dt * ftmp
       utilde =  K + L*integrator.sqdt
       ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
-      mil_correction = ggprime .* (integrator.W.dW.^2 .- dt)./2
+      mil_correction = ggprime .* (integrator.W.dW.^2 .- abs(dt))./2
       gtmp += mil_correction
     elseif alg_interpretation(alg) == :Stratonovich ||
            typeof(cache) <: ImplicitEulerHeunConstantCache
@@ -158,7 +158,7 @@ end
       @.. z = uprev + dt * tmp + integrator.sqdt * gtmp
       integrator.g(gtmp3,z,p,t)
       @.. gtmp3 = (gtmp3-gtmp)/(integrator.sqdt) # ggprime approximation
-      @.. gtmp2 += gtmp3*(dW.^2 - dt)/2
+      @.. gtmp2 += gtmp3*(dW.^2 - abs(dt))/2
     elseif alg_interpretation(alg) == :Stratonovich
       @.. z = uprev + integrator.sqdt * gtmp
       integrator.g(gtmp3,z,p,t)