# Solver throwing LinearAlgebra.SingularException(1) error for my ODE System

**URL:** <https://discourse.julialang.org/t/solver-throwing-linearalgebra-singularexception-1-error-for-my-ode-system/88281>\
**Category:** Numerics\
**Tags:** question, diffeq\
**Created:** [October 5, 2022, 10:09am UTC](https://discourse.julialang.org/t/solver-throwing-linearalgebra-singularexception-1-error-for-my-ode-system/88281 "2022-10-05T10:09:10Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![hbooth](https://avatars.discourse-cdn.com/v4/letter/h/e19b73/32.png) [@hbooth](https://discourse.julialang.org/u/hbooth)\
**Post date:** [October 5, 2022, 10:09am UTC](https://discourse.julialang.org/t/solver-throwing-linearalgebra-singularexception-1-error-for-my-ode-system/88281/1 "2022-10-05T10:09:10Z")

</div>

Hi,

When I am trying to solve my ODE system I get a LinearAlgebra.SingularException(1) error and I am not sure why or what I can do about it. Here is my system:

```julia
using DifferentialEquations

const Nc = 100
const L = 1.
const θ = 5.
const c0 = 1.
const λm = 0.4

const tissue = range(0,L,length = Nc)

m(x) = c0*exp(-x/λm)

σ(I) = 1/(1+exp(θ-θ*I)) 

function ode_system!(dg,g,p,t)

    w, λg, = p

    @inbounds for j in 1:Nc
        x = tissue[j]
        dg[1,j] = σ(w[1,1]*g[1,j]+ w[1,2]*g[2,j] + w[1,3]*g[3,j] + w[1,4]*m(x)) - λg[1]*g[1,j]
        dg[2,j] = σ(w[2,1]*g[1,j]+ w[2,2]*g[2,j] + w[2,3]*g[3,j] + w[2,4]*m(x)) - λg[2]*g[2,j]
        dg[3,j] = σ(w[3,1]*g[1,j]+ w[3,2]*g[2,j] + w[3,3]*g[3,j] + w[3,4]*m(x)) - λg[3]*g[3,j]
    end
    
end

w_ex = [-3.43445 1.66673 -0.562781 1.20563;
                0.0 0.626765 -0.782956 0.0;
                -7.24169 -6.14532 0.0 0.0]

λg_ex = 0.05 .* ones(3)

p = (w_ex,λg_ex)

g0 = 0.01 .* ones(3,Nc)
    
prob = ODEProblem(ode_system!,g0,(0,Inf),p)

sol = solve(prob,AutoTsit5(Rosenbrock23()),isoutofdomain=(u,p,t) -> any(x -> x < 0, u), callback = TerminateSteadyState(1e-5,1e-3),maxiters = 1e3, verbose = false, save_everystep = false)

```

Any help would be appreciated.

Thanks

---

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [October 5, 2022, 1:36pm UTC](https://discourse.julialang.org/t/solver-throwing-linearalgebra-singularexception-1-error-for-my-ode-system/88281/2 "2022-10-05T13:36:00Z")

</div>

I cannot reproduce SingularException. But if you try a few other solvers, they all error because `dt<dtmin` as well. This indicates that your system likely incorrect.

```julia
julia> sol = solve(prob,Rodas5(),isoutofdomain=(u,p,t) -> any(x -> x < 0, u), reltol=1e-5, abstol=1e-7, callback = TerminateSteadyState(1e-5,1e-3), maxiters = 10^3);
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/src/julia/SciMLBase/src/integrator_interface.jl:523

julia> sol = solve(prob, CVODE_BDF(), isoutofdomain=(u,p,t) -> any(x -> x < 0, u), reltol=1e-5, abstol=1e-7, callback = TerminateSteadyState(1e-5,1e-3), maxiters = 10^3);
┌ Warning: The isoutofdomain argument is ignored by CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0).
└ @ SciMLBase ~/src/julia/SciMLBase/src/utils.jl:262
┌ Warning: https://diffeq.sciml.ai/dev/basics/compatibility_chart/
└ @ SciMLBase ~/src/julia/SciMLBase/src/utils.jl:274

[CVODES ERROR] CVode
  At t = 0 and h = inf, the corrector convergence test failed repeatedly or with |h| = hmin.

```

---

<div class="post-metadata">

**Author:** ![hbooth](https://avatars.discourse-cdn.com/v4/letter/h/e19b73/32.png) [@hbooth](https://discourse.julialang.org/u/hbooth)\
**Post date:** [October 5, 2022, 1:41pm UTC](https://discourse.julialang.org/t/solver-throwing-linearalgebra-singularexception-1-error-for-my-ode-system/88281/3 "2022-10-05T13:41:56Z")

</div>

Hi. I tried it with the TRBDF solver and it terminates successfully, so I am not sure it is the system…

```julia
sol = solve(prob,TRBDF2(), callback = TerminateSteadyState(1e-5,1e-3),maxiters = 1e3, verbose = false, save_everystep = false)

retcode: Terminated
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
   0.0
 274.93965955950193
u: 2-element Vector{Matrix{Float64}}:
 [0.01 0.01 … 0.01 0.01; 0.01 0.01 … 0.01 0.01; 0.01 0.01 … 0.01 0.01]
 [9.614177032652817 9.60550761535028 … 9.283026231664332 9.282193242446969; 19.681963743795468 19.681929558547157 … 19.654414208139823 19.654198715172864; 1.0325728578241648e-8 1.0336950034137993e-8 … 1.96255649162726e-8 1.9717726818553073e-8]

```

I am not sure why you don’t manage to reproduce the singular exception…this are my packages :

```julia
pkg> status
Status `~/.julia/environments/v1.8/Project.toml`
  [537997a7] AbstractPlotting v0.18.3
  [6e4b80f9] BenchmarkTools v1.3.1
  [3da002f7] ColorTypes v0.11.4
  [2b5f629d] DiffEqBase v6.105.0
  [459566f4] DiffEqCallbacks v2.24.1
  [aae7a2af] DiffEqFlux v1.52.0
⌃ [071ae1c0] DiffEqGPU v1.19.0
⌃ [0c46a032] DifferentialEquations v7.4.0
  [b4f34e82] Distances v0.10.7
  [31c24e10] Distributions v0.25.75
  [5b8099bc] DomainSets v0.5.13
⌃ [634d3b9d] DrWatson v2.9.1
  [aaaaaaaa] DynamicAxisWarping v0.4.11
  [d8961e24] FindPeaks1D v0.1.7 `https://github.com/ymtoo/FindPeaks1D.jl.git#master`
⌃ [e9467ef8] GLMakie v0.4.5
  [e91730f6] Hungarian v0.6.0
  [7073ff75] IJulia v1.23.3
  [40713840] IncompleteLU v0.2.0
  [94925ecb] MethodOfLines v0.5.0
⌃ [961ee093] ModelingToolkit v8.23.0
  [7e02d93a] OptimalTransport v0.3.19
⌃ [91a5bcdd] Plots v1.32.0
  [92933f4c] ProgressMeter v1.7.2
  [90137ffa] StaticArrays v1.5.9
  [2913bbd2] StatsBase v0.33.21
  [f3b207a7] StatsPlots v0.15.3
  [c3572dad] Sundials v4.10.1
  [0c5d862f] Symbolics v4.10.4
Info Packages marked with ⌃ have new versions available

```

---

<div class="post-metadata">

**Author:** ![hbooth](https://avatars.discourse-cdn.com/v4/letter/h/e19b73/32.png) [@hbooth](https://discourse.julialang.org/u/hbooth)\
**Post date:** [October 5, 2022, 2:08pm UTC](https://discourse.julialang.org/t/solver-throwing-linearalgebra-singularexception-1-error-for-my-ode-system/88281/4 "2022-10-05T14:08:12Z")

</div>

I should add, if I drop the relative tolerance on the `TerminateSteadyState` call the original solver I tried `AutoTsit5(Rosenbrock23())` does terminate:

```julia
solve(prob,AutoTsit5(Rosenbrock23()), callback = TerminateSteadyState(1e-5,1e-2),maxiters = 1e3, verbose = false, save_everystep = false)

retcode: Terminated
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
  0.0
 94.09830008270396
u: 2-element Vector{Matrix{Float64}}:
 [0.01 0.01 … 0.01 0.01; 0.01 0.01 … 0.01 0.01; 0.01 0.01 … 0.01 0.01]
 [0.4211883433342539 0.4135235406632426 … 0.15502069856838743 0.15446894770089728; 0.2897922503571214 0.2897882442706354 … 0.2867584062496022 0.2867356485393743; 9.202138757164838e-5 9.212647292613664e-5 … 0.00018143431424400678 0.00018217531267239543]

```

This makes me think the error is something to do with the tolerances - the third component du goes extremely small - maybe its something to do with the relative condition in the default terminating function used in the callback - `allDerivPass` ?

Its not so much of a problem if the problem is unstable, but I need the solver to reliably return a retcode rather throw an error, as I am performing a parameter search on the matrices `w_ex`…

---

<div class="post-metadata">

**Author:** ![hbooth](https://avatars.discourse-cdn.com/v4/letter/h/e19b73/32.png) [@hbooth](https://discourse.julialang.org/u/hbooth)\
**Post date:** [October 6, 2022, 11:18am UTC](https://discourse.julialang.org/t/solver-throwing-linearalgebra-singularexception-1-error-for-my-ode-system/88281/5 "2022-10-06T11:18:33Z")

</div>

Another peculiarity I have found in this example relates to the use of `isoutofdomain` - suppose we change the matrix to

`w_ex = [0.0 0.0 0.0 6.65803; -0.101676 0.0 1.00214 0.0; 0.0 -0.0664579 0.599044 0.0]`

Then when solving for the steady state without domain control I get

```julia
solve(prob,AutoTsit5(Rosenbrock23()), callback = TerminateSteadyState(1e-5,1e-3),maxiters = 1e3,save_everystep = false,verbose = false)

retcode: Terminated
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
    0.0
 5437.216599828782
u: 2-element Vector{Matrix{Float64}}:
 [0.01 0.01 … 0.01 0.01; 0.01 0.01 … 0.01 0.01; 0.01 0.01 … 0.01 0.01]
 [19.99999999998965 19.999999999976268 … 1.9995154232199426 1.8772194358947054; 20.0 20.0 … 0.2353371973781851 0.2410830963626622; 20.0 20.0 … 0.31649911264762154 0.3089632849374852]

```

however with domain control…

```julia
solve(prob,AutoTsit5(Rosenbrock23()),isoutofdomain=(u,p,t) -> any(x -> x < 0, u), callback = TerminateSteadyState(1e-5,1e-3),maxiters = 1e3,save_everystep = false,verbose = false)

retcode: MaxIters
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
     0.0
 24128.188732221657
u: 2-element Vector{Matrix{Float64}}:
 [0.01 0.01 … 0.01 0.01; 0.01 0.01 … 0.01 0.01; 0.01 0.01 … 0.01 0.01]
 [19.99999999998965 19.999999999976268 … 1.9995154232199426 1.8772194358947054; 20.0 20.0 … 0.23533719737818506 0.24108309636266156; 20.0 20.0 … 0.31649911264762154 0.30896328493748465]

```

---

<div class="post-metadata">

**Author:** ![hbooth](https://avatars.discourse-cdn.com/v4/letter/h/e19b73/32.png) [@hbooth](https://discourse.julialang.org/u/hbooth)\
**Post date:** [October 13, 2022, 8:43am UTC](https://discourse.julialang.org/t/solver-throwing-linearalgebra-singularexception-1-error-for-my-ode-system/88281/6 "2022-10-13T08:43:10Z")

</div>

@ChrisRackauckas do you have any insights here?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 13, 2022, 9:19am UTC](https://discourse.julialang.org/t/solver-throwing-linearalgebra-singularexception-1-error-for-my-ode-system/88281/7 "2022-10-13T09:19:33Z")

</div>

That means that the true solution is going negative. You’ve tried lower tolerances?

> [@hbooth](#):
>
> Its not so much of a problem if the problem is unstable, but I need the solver to reliably return a retcode rather throw an error, as I am performing a parameter search on the matrices `w_ex`…

Yes, it’s in my list of things to do rather soon.

> <https://github.com/SciML/LinearSolve.jl/pull/207>
>
> @ChrisRackauckas acceptable? I'll add a check option to UMFPACK for 1.9, and fig…ure out how I want to do this in KLU. But this should suffice for now, we can always set \`check = false\` default later once those checks are integrated upstream.

But for that, I am fixing up the retcodes a bit before @rayegun goes to add that all to LinearSolve.

> <https://github.com/SciML/SciMLBase.jl/pull/274>
>
> I think I found a way that is relatively safe. In all of our code throughout Sci…ML, I noticed that the only thing that was truly done was \`sol.retcode == :Success\`, etc. in which case a relatively small amount of functions put over a retcode enum could be non-breaking to all functional uses of the previous retcode. If this is the case, we may finally have a good way to upgrade the one piece I find most hideous in the interface.
> 
> Fixes https://github.com/SciML/DifferentialEquations.jl/issues/867

That’s my today. It’ll also help the solvers in that if they warn exit they can reject the step instead of erroring, and then adaptivity could do its trick to try other steps.