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[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;
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[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);
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[[2]]));
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[1]), length(MATraw))))
println("Done")