On the below problem, I had two questions
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.
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) k1=0.01 k2=0.01 Dt=2000; du = Dt - u*k1; method= "abs" if method=="max" du = u*k1- k2* (max(0,u))^0.8; elseif method=="abs" du = u*k1- k2* abs(u)^0.8; end return du end println("Starting....") #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() #cb= PositiveDomain() solvers= ["sundials", "kencarp4offauto" , "kencarp4auto"] solver_choose= solvers 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); end p= 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[])); else sol = solve(prob, abstol=tol, reltol=tol,solver, saveat=0.05); end 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), length(MATraw)))) println("Done")