Stiff system (of >100 equations) - Rodas5 - multiple runs - instability or not

Specific Domains Modelling & Simulations
,
1

Dear all,

I have a stiff system of >100 equations (which I’m not allowed to post anywhere at the moment). I tried

solve(my_prob, Rodas5(autodiff=false),reltol=1e-8,abstol=1e-8,saveat=1.0) # let's call it  method1

and

solve(my_prob, Rodas5(autodiff=false),reltol=1e-8,abstol=1e-8,saveat=1.0, dt=1e-5, adaptive=true )  # let's call it  method 2

where

my_prob= ODEProblem(my_system!,u0,tspan_m,param)

In the very beginning I set

import Random;
Random.seed!(1234);

However, when I run simulation multiple times with the same initial conditions and same parameter sets, I get different output in a sense that sometimes I get the following errors:

  1. For method 1:
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq /Users/polina/.julia/packages/OrdinaryDiffEq/LQQYm/src/solve.jl:456
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ DiffEqBase /Users/polina/.julia/packages/DiffEqBase/TjqaN/src/integrator_interface.jl:323

and

┌ Warning: Instability detected. Aborting
└ @ DiffEqBase /Users/polina/.julia/packages/DiffEqBase/TjqaN/src/integrator_interface.jl:349
  1. For method 2:
┌ Warning: Instability detected. Aborting
└ @ DiffEqBase /Users/polina/.julia/packages/DiffEqBase/TjqaN/src/integrator_interface.jl:349

Here is an example of how I found that:

du = similar(u0);
for i in 1:20
    param = append!(param_sets[10,1:11],fixedForNowParam[1]) 
    tspan_m = (0.0,1200.0)
    my_prob = ODEProblem(my_system!,u0,tspan_m,param) 
    sol_m = solve(my_prob, Rodas5(autodiff=false),reltol=1e-8,abstol=1e-8,saveat=1.0, dt=1e-5, adaptive=true ) 
    print(i)

end

I understand that my question may be not descriptive enough without the “my_system” defined here, but maybe someone encountered similar problem and knows where can I find info to understand why it happens. I am wondering may it be related to adaptive time step or some other parameters of Rodas5? Several other methods that I tried weren’t really efficient on my system.
Thank you very much for the help and insights!

Status ~/.julia/environments/v1.4/Project.toml
[6e4b80f9] BenchmarkTools v0.5.0
[31c91b34] DiffEqBenchmarks v0.1.0 #master (https://github.com/JuliaDiffEq/DiffEqBenchmarks.jl)
[f3b72e0c] DiffEqDevTools v2.21.0
[aae7a2af] DiffEqFlux v1.12.0
[1130ab10] DiffEqParamEstim v1.13.0
[41bf760c] DiffEqSensitivity v6.10.0
[ef61062a] DiffEqUncertainty v1.4.1
[0c46a032] DifferentialEquations v6.14.0
[31c24e10] Distributions v0.22.6
[587475ba] Flux v0.10.4
[7073ff75] IJulia v1.21.1
[7f56f5a3] LSODA v0.6.1
[c030b06c] ODE v2.8.0
[54ca160b] ODEInterface v0.4.6
[09606e27] ODEInterfaceDiffEq v3.7.0
[429524aa] Optim v0.20.1
[1dea7af3] OrdinaryDiffEq v5.39.1
[65888b18] ParameterizedFunctions v5.0.3
[91a5bcdd] Plots v0.29.9
[90137ffa] StaticArrays v0.12.1
[4c63d2b9] StatsFuns v0.9.5
[c3572dad] Sundials v3.9.0
[44d3d7a6] Weave v0.9.4
[0518478a] deSolveDiffEq v0.1.0

2

It’s hard to debug this without knowing the model, but all signs are pointing to model mispecification. You tried other methods for stiff equations and everything had the same behavior?

3

Thank you @ChrisRackauckas !

I tried the following methods:
1)

solve(my_prob, Rosenbrock23(),reltol=1e-8,abstol=1e-8,saveat=1.0)
  • it gives similar behaviour with instabilities in some runs and adequate results in other runs (u0 and params are the same)
solve(my_prob, Rosenbrock23(autodiff=false),reltol=1e-8,abstol=1e-8,saveat=1.0)
  • either good runs or Domain error or LinearAlgebra.SingularException(163) [in all cases it takes significantly longer time to run than Rodas5]
solve(my_prob, TRBDF2(),reltol=1e-8,abstol=1e-8,saveat=1.0)
  • it gives similar behaviour with instabilities in some runs and adequate results in other runs (u0 and params are the same)
solve(my_prob, TRBDF2(autodiff=false),reltol=1e-8,abstol=1e-8,saveat=1.0)
  • either good runs or LinearAlgebra.SingularException(163)
solve(my_prob, Rodas4P(autodiff=false),reltol=1e-8,abstol=1e-8,saveat=1.0)
  • either good runs or “Warning: Instability detected. Aborting” or “Warning: Automatic dt set the starting dt as NaN, causing instability.”
solve(my_prob, Rodas4P(),reltol=1e-8,abstol=1e-8,saveat=1.0)
  • either good runs or “Warning: Instability detected. Aborting” or "DomainError with -0.002495613279934969: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar. "
solve(my_prob, Kvaerno5(),reltol=1e-8,abstol=1e-8,saveat=1.0)
  • either good runs or “Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.”
solve(my_prob, KenCarp4(),reltol=1e-8,abstol=1e-8,saveat=1.0)
  • either good runs or “Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.” or “DomainError with -3.212142508058485e-5: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar.”

My main question:
I am trying to understand whether the problem is in my system of equations, or solver, or both. What is not clear to me, is that if I have a misspecification in the system, why does the problem occur only in ~13% of the cases (in case of Rodas5, or in ~10% - ~25% in other methods)? If there is misspecification in the system, shouldn’t it fail in every run? Why do I even get quite a large percent of simulations which run well (when all parameters and initial conditions are the same to those simulations which failed)? I am searching for the theoretical explanation for these why… but maybe it’s really system dependent…

Thank you very much!

4

Are you doing something that assumes time is monotonic, like using a buffer to carry values over from one f call to the next? That kind of stuff breaks the assumptions of adaptivity.

1 Like
5

I didn’t expect it to be a case, but it sounds so close to the answer to my questions, especially given that restarting Jupyter kernel improves the expectation of good run. I have some local variables, maybe they can be the cause. I will have a more close look at how to use the best practices of Julia, but I have to use only those specifications that are compatible with diffeqpy.

Thank you very much again! I think this is the hint I was looking for!

6

I solved it, it was such a stupid mistake! Thank you so much for your advice, @ChrisRackauckas !

So apparently as this model has quite a long history by now, it became a bit redundant, and a few derivatives are zeros now, so my mistake was that I just commented them while working on the equations themselves, instead of explicit assignment of these derivatives to zeros. Surprisingly, overall results weren’t affected, apart from that weird behaviour with instabilities at some runs.

In addition to it, I also replaced du = similar(u0); to du = zeros(length(u0)); [just in case because similar(u0) gave a vector with some non-zero numbers and I wasn’t sure how important is it], and included “return nothing”.

Just to be double-sure that now my code does what I expect it to do, if you have a little bit of time and interest to the topic @ChrisRackauckas, could you please have a quick look at my little piece of code below? It’s a toy model which is solvable and it runs well, it doesn’t mean anything, but it features all the aspects that I found doubtable after going through Julia tutorials and docs.

function my_system!(du,u,p,t)
    
    a = 0.18
    b = 1e-5
    c = 0.0001
    d = a + b/c
    
    j1 = a * ( u[1] * u[2]*exp(-b/a) ) 
    j2 = p[1]*((u[1])^0.1 - u[2])^(-1)
    j3 = p[2]*u[3]*log(clamp(u[1],1e-15,u[1])/clamp(u[2],1e-15,u[2])) 
    
    du[1] =  0.5*(j1 + j2)  
    du[2] =  j2 - j3
    du[3] = 0
        
    return nothing 
end
    

ini1 = 0.001;
ini2 = 0.5;
ini3 = 1e-5;

u0 = [ini1,ini2,ini3];

du = zeros(length(u0));
    
param = [0.0003, 1.2]

tspan = (0.0,100.0)

my_prob = ODEProblem(my_system!,u0,tspan,param)  

sol = solve(my_prob, Rodas5(autodiff=false),reltol=1e-8,abstol=1e-8,saveat=1.0, dt=1e-5, adaptive=true )

Thank you very much again!

7

Yeah, that looks good, though you don’t have to do adaptive = true since that’s the default. Also, using autodiff here should be fine, and I’m not sure setting the initial dt is necessary.

3 Likes
8

Great! Thank you very much @ChrisRackauckas !