How can I implement positive Domain for sundials CVODE_BDF ? Also, try to compare both isoutofdomain , with julia solvers, one question there, Thanks

On the below problem, I had two questions

  1. When applying kencarp4 isoutofdomain is being applied to both states rather than the second du2 like I want to , so the IC is being populated for both states, but I expect only the second being affect, why is this? I also get a dt < dt min error , some type of conergence error or system unstable, I do not see this if I remove isoutofdomain option.

  2. I read the documentation of positive domain, but is still not clear what I need to pass to u argument, in positive domain, in this case I want to say if u2 is negative start halving the steps until you get + for u2 , but I am assuming they will ignore the abs () or I want it to be never , I want u2 to stay negative essentially, if the initial IC is negative.

I am confused if I will need general domain here instead. Again thank you, I would appreciate your help in the given problem , I have not been able to find an example that matches what I want to do, hope my question and problem is clear, thanks.

My main objective is to understand how to use Positive Domain , beyond what I read in the documentation, and I thought while I am at it , comparing with isoutofdomain was a good idea, I look forward to your response, thank you!!

using DifferentialEquations
using Sundials
using DiffEqCallbacks
function  ode_functionminie(du,u, p,t)


du[1] = Dt - u[1]*k1; 
method= "abs"

if method=="max"
    du[2] = u[1]*k1- k2* (max(0,u[2]))^0.8; 

elseif method=="abs"
    du[2] = u[1]*k1- k2* abs(u[2])^0.8; 


return du



#isoutofdomain=(u,p,t) -> any(x -> (x < 0), u)
#cb=PositiveDomain(u=nothing; save=true, abstol=nothing, scalefactor=nothing)
#u: A prototype of the state vector of the integrator. A copy of it is saved and extrapolated values are written to it. If it is not specified every application of the callback allocates a new copy of the state vector.

#cb= PositiveDomain()

solvers= ["sundials", "kencarp4offauto" , "kencarp4auto"]
solver_choose= solvers[2]

if solver_choose=="lsoda"
    solver= lsoda();

    #sol = solve(prob, abstol=1e-6, reltol=1e-6, lsoda())
elseif solver_choose=="sundials"
    solver= CVODE_BDF();

elseif solver_choose=="kencarp4auto"
    solver = KenCarp4();

elseif solver_choose=="kencarp4offauto"
    solver = KenCarp4(autodiff=false);

uo_e=[0.0, -654]
tol= 1e-6
t_start= 0.0;
sim_time_e= 24.0;

prob = ODEProblem(ode_functionminie,uo_e,(t_start,sim_time_e), p)
if solver_choose=="sundials"
    sol = solve(prob, abstol=tol, reltol=tol,solver, dtmax=0.1, cb);
elseif solver_choose=="kencarp4offauto" || solver_choose=="kencarp4auto"
     sol = solve(prob, abstol=tol, reltol=tol,solver, dtmax=0.1, isoutofdomain=(u,p,t) -> any(x -> (x < 0), u[[2]]));
    sol = solve(prob, abstol=tol, reltol=tol,solver, saveat=0.05);

times= 0.0:1.0:12.0
time_output, output= collect(times), sol(times).u
final_output= Array(sol(times)')
MATraw= output;
final_output2= permutedims(reshape(hcat(MATraw...), (length(MATraw[1]), length(MATraw))))


You can enclose your code in triple backticks ```


Thanks Nilshg!

1 Like

It looks like you’re applying isoutofdomain to keep the solution positive when the true solution actually goes negative. DifferentialEquations.jl solves the differential equation, so if it goes negative, we will error saying that what was specified was impossible (here, dt goes to dtmin because it keeps rejecting). Check your model.

1 Like

Chris, thanks for looking into it, that make sense, thanks for the clarification, help!