LinearAlgebra.SingularException(2) when solving matrix equation which is solvable in Maple


I’m doing an assignment where I use my own implementation of Implicit Eulers with adaptive step size using error estimation. For this, I have made my own implementation of the Newton solver. In one ODE system, my numerical solver fails to solve the following matrix equation.

dRdt = [2.7347115660463565e46 1.4462465159832052e46 7.06637105328182e46; 5.469423132092713e46 2.8924930319664105e46 1.413274210656364e47; -3.658477011433253e48 -1.9347779478036194e48 -9.453339201714811e48]
R = [1.0026696232160644e47, 2.0053392464321287e47, -1.3413640444362068e49]


This returns the error

ERROR: LoadError: LinearAlgebra.SingularException(2)

I have tried to solve the equation in Maple and it solves it without an issue (see image below)

Does anyone know what is wrong?

The system is ill-conditioned

julia> dRdt = [2.7347115660463565e46 1.4462465159832052e46 7.06637105328182e46
               5.469423132092713e46 2.8924930319664105e46 1.413274210656364e47
              -3.658477011433253e48 -1.9347779478036194e48 -9.453339201714811e48];

julia> cond(dRdt)
1 Like

I would factor out 10^{46} from the system :wink:

I don’t think that will affect the condition number, will it?

julia> svdvals(dRdt)
3-element Vector{Float64}:

julia> svdvals(dRdt .* 1.0e-46)
3-element Vector{Float64}:

I guess I should have looked at the results before I said that it wouldn’t affect the condition number. The condition number for the scaled matrix is on the order of 10^19 instead of 10^33 but it is still computationally singular.

The condition number won’t change in infinite precision, but that’s not what we’re dealing with here. Anyhow, as @dmbates observes, 10^19 is still bad enough to make most solvers complain.

Indeed, my bad. You probably only gain floating-point representation space.
My guess is that Maple computes the solution analytically (relatively easy for 3x3).

Thank you all for your answers they were very helpful to understand what was going on! I ended up tweaking the tolerances of my solver, which avoided that ill-conditioned system and solved the ODEs.