catch bad initdt

Author: ChrisRackauckas

src/initdt.jl

@@ -10,6 +10,38 @@
   f₀ = zeros(u0./t)
   f(t,u0,f₀)
 
+
+  #=
+  Try/catch around the linear solving. This will catch singular matrices defined
+  by DAEs and thus we use the tType(1//10^(6)) default from Hairer. Note that
+  this will not always catch singular matrices, an example from Andreas:
+
+  julia> A = fill(rand(), 2, 2)
+  2×2 Array{Float64,2}:
+   0.637947  0.637947
+   0.637947  0.637947
+
+  julia> inv(A)
+  2×2 Array{Float64,2}:
+    9.0072e15  -9.0072e15
+   -9.0072e15   9.0072e15
+
+  The only way to make this more correct is to check
+
+  issingular(A) = rank(A) < min(size(A)...)
+
+  but that would introduce another svdfact in rank (which may not be possible
+  anyways if the mass_matrix is not actually an array). Instead we stick to the
+  user-chosen factorization. Sometimes this will cause `ftmp` to be absurdly
+  large like shown there, but that later gets caught in the quick estimates
+  below which then makes it spit out the default
+
+  dt₀ = tType(1//10^(6))
+
+  so that is a a cheaper way to get to the same place than an svdfact and
+  still works for matrix-free definitions of the mass matrix.
+  =#
+
   if prob.mass_matrix != I
     ftmp = similar(f₀)
     try
@@ -33,6 +65,13 @@
     dt₀ = tType((d₀/d₁)/100)
   end
   dt₀ = min(dt₀,dtmax_tdir)
+
+  if tType <: AbstractFloat && dt₀ < 10eps(tType)
+    # This catches Andreas' non-singular example
+    # should act like it's singular
+    return tType(1//10^(6))
+  end
+
   dt₀_tdir = tdir*dt₀
 
   u₁ = similar(u0) # required by DEDataArray

test/mass_matrix_tests.jl

@@ -89,3 +89,35 @@ sol = solve(prob,   TRBDF2(),adaptive=false,dt=1/10)
 sol2 = solve(prob2, TRBDF2(),adaptive=false,dt=1/10)
 
 =#
+
+# Singular mass matrices
+
+function f!(t, u, du)
+    du[1] = u[2]
+    du[2] = u[2] - 1.
+    return
+end
+
+u0 = [0.,1.]
+tspan = (0.0, 1.0)
+
+M = zeros(2,2)
+M[1,1] = 1.
+
+m_ode_prob = ODEProblem(f!, u0, tspan, mass_matrix=M)
+sol = solve(m_ode_prob, Rosenbrock23())
+
+M = [0.637947  0.637947
+     0.637947  0.637947]
+
+inv(M) # not caught as singular
+
+function f2!(t, u, du)
+    du[1] = u[2]
+    du[2] = u[1]
+    return
+end
+u0 = zeros(2)
+
+m_ode_prob = ODEProblem(f2!, u0, tspan, mass_matrix=M)
+sol = solve(m_ode_prob, Rosenbrock23())