multiply initial dt by random number to prevent numerical coincendences

lib/OrdinaryDiffEqCore/src/initdt.jl

@@ -9,10 +9,12 @@
     f = prob.f
     p = integrator.p
     oneunit_tType = oneunit(_tType)
+    # a number that's about 1 that we multiply by to prevent numerical coincedences
+    arbitrary_const = oneunit_tType*(93//83)
     dtmax_tdir = tdir * dtmax
 
     dtmin = nextfloat(max(integrator.opts.dtmin, eps(t)))
-    smalldt = max(dtmin, convert(_tType, oneunit_tType * 1 // 10^(6)))
+    smalldt = max(dtmin, convert(_tType, arbitrary_const * 1 // 10^(6)))
 
     if integrator.isdae
         result_dt = tdir * max(smalldt, dtmin)
@@ -132,10 +134,10 @@
 
     if u0 isa Array
         @inbounds @simd ivdep for i in eachindex(u0)
-            tmp[i] = f₀[i] / sk[i] * oneunit_tType
+            tmp[i] = f₀[i] / sk[i] * arbitrary_const
         end
     else
-        @.. broadcast = false tmp = f₀ / sk * oneunit_tType
+        @.. broadcast = false tmp = f₀ / sk * arbitrary_const
     end
 
     d₁ = internalnorm(tmp, t)
@@ -156,20 +158,15 @@
             (d₁ < 1 // 10^(5)), smalldt,
         convert(
             _tType,
-            oneunit_tType * DiffEqBase.value(
+            arbitrary_const * DiffEqBase.value(
                 (d₀ / d₁) /
                     100
             )
         )
     )
-    # if d₀ < 1//10^(5) || d₁ < 1//10^(5)
-    #   dt₀ = smalldt
-    # else
-    #   dt₀ = convert(_tType,oneunit_tType*(d₀/d₁)/100)
-    # end
     dt₀ = min(dt₀, dtmax_tdir)
 
-    if typeof(one(_tType)) <: AbstractFloat && dt₀ < 10eps(_tType) * oneunit(_tType)
+    if typeof(one(_tType)) <: AbstractFloat && dt₀ < 10eps(_tType) * arbitrary_const
         # This catches Andreas' non-singular example
         # should act like it's singular
         result_dt = tdir * max(smalldt, dtmin)
@@ -209,22 +206,22 @@
 
     if u0 isa Array
         @inbounds @simd ivdep for i in eachindex(u0)
-            tmp[i] = (f₁[i] - f₀[i]) / sk[i] * oneunit_tType
+            tmp[i] = (f₁[i] - f₀[i]) / sk[i] * arbitrary_const
         end
     else
-        @.. broadcast = false tmp = (f₁ - f₀) / sk * oneunit_tType
+        @.. tmp = (f₁ - f₀) / sk * arbitrary_const
     end
 
-    d₂ = internalnorm(tmp, t) / dt₀ * oneunit_tType
+    d₂ = internalnorm(tmp, t) / dt₀ * arbitrary_const
     # Hairer has d₂ = sqrt(sum(abs2,tmp))/dt₀, note the lack of norm correction
 
     max_d₁d₂ = max(d₁, d₂)
     if max_d₁d₂ <= 1 // Int64(10)^(15)
-        dt₁ = max(convert(_tType, oneunit_tType * 1 // 10^(6)), dt₀ * 1 // 10^(3))
+        dt₁ = max(convert(_tType, arbitrary_const * 1 // 10^(6)), dt₀ * 1 // 10^(3))
     else
         dt₁ = convert(
             _tType,
-            oneunit_tType *
+            arbitrary_const *
                 DiffEqBase.value(
                 10.0^(
                     -(2 + log10(max_d₁d₂)) /
@@ -277,10 +274,12 @@ end
     f = prob.f
     p = prob.p
     oneunit_tType = oneunit(_tType)
+    # a number that's about 1 that we multiply by to prevent numerical coincedences
+    arbitrary_const = oneunit_tType*(93//83)
     dtmax_tdir = tdir * dtmax
 
     dtmin = nextfloat(max(integrator.opts.dtmin, eps(t)))
-    smalldt = max(dtmin, convert(_tType, oneunit_tType * 1 // 10^(6)))
+    smalldt = max(dtmin, convert(_tType, arbitrary_const * 1 // 10^(6)))
 
     if integrator.isdae
         return tdir * max(smalldt, dtmin)
@@ -303,12 +302,12 @@ end
         throw(TypeNotConstantError(inferredtype, typeof(f₀)))
     end
 
-    d₁ = internalnorm(f₀ ./ sk .* oneunit_tType, t)
+    d₁ = internalnorm(f₀ ./ sk .* arbitrary_const, t)
 
     if d₀ < 1 // 10^(5) || d₁ < 1 // 10^(5)
         dt₀ = smalldt
     else
-        dt₀ = convert(_tType, oneunit_tType * DiffEqBase.value((d₀ / d₁) / 100))
+        dt₀ = convert(_tType, arbitrary_const * DiffEqBase.value((d₀ / d₁) / 100))
     end
     dt₀ = min(dt₀, dtmax_tdir)
     dt₀_tdir = tdir * dt₀
@@ -321,16 +320,16 @@ end
     # Avoids AD issues
     f₀ == f₁ && return tdir * max(dtmin, 100dt₀)
 
-    d₂ = internalnorm((f₁ .- f₀) ./ sk .* oneunit_tType, t) / dt₀ * oneunit_tType
+    d₂ = internalnorm((f₁ .- f₀) ./ sk .* arbitrary_const, t) / dt₀ * arbitrary_const
 
     max_d₁d₂ = max(d₁, d₂)
     if max_d₁d₂ <= 1 // Int64(10)^(15)
-        dt₁ = max(smalldt, dt₀ * 1 // 10^(3))
+        dt₁ = max(smalldt, dt₀ * 1 // 10^3)
     else
         dt₁ = _tType(
-            oneunit_tType *
+            arbitrary_const *
                 DiffEqBase.value(
-                10^(
+                exp10(
                     -(2 + log10(max_d₁d₂)) /
                         get_current_alg_order(
                         integrator.alg,
@@ -354,6 +353,6 @@ end
     _tType = eltype(tType)
     tspan = prob.tspan
     init_dt = abs(tspan[2] - tspan[1])
-    init_dt = isfinite(init_dt) ? init_dt : oneunit(_tType)
-    return convert(_tType, init_dt * 1 // 10^(6))
+    init_dt = isfinite(init_dt) ? init_dt : oneunit(_tType) * 84 // 83 * 1 // 10^6
+    return convert(_tType, init_dt)
 end