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

Specific Domains Numerics
,
1

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:

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

2

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> 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.
3

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

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 :

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
4

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

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…

5

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

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…

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]
6

@ChrisRackauckas do you have any insights here?

7

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

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

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

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.

1 Like