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

General Usage
1

Hi

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]

dRdt\R

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)
image

Does anyone know what is wrong?

2

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)
2.885773849638737e32
1 Like
3

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

4

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

julia> svdvals(dRdt)
3-element Vector{Float64}:
 1.0321007400020088e49
 1.432533639792516e33
 3.576512900105512e16

julia> svdvals(dRdt .* 1.0e-46)
3-element Vector{Float64}:
 1032.1007400020085
    1.2316027007572554e-13
    2.9065558295847953e-16
5

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.

7

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.

8

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

9

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.