Here is an extreme case I posed. Here, I end up with a domain error, because I start at 0.1 ( any higher ) one parameter will give a domain error, I chose du[2] to be 0 or close to it with a negative value. I get a domain error because of floating rounding of or precision as I describe it. I wanted to understand best practices of how to solve this, Chris suggested me earlier and it was very helpful to add abs ( ) , I have seen also practices of using max(0, (…)) , (…) being the negative value due to rounding i.e param_cal= (1-u[2]/0.1)^0.8. I have not tried this , but if param_cal was= (0.1- u[2])^0.8 , would it have also given similar error IC of u[2] is 0.1, or is it the division that creates this issue.
Hope you find the message clear, thank you again!
Christian
using DifferentialEquations
using Sundials
function ode_functionminie(du,u, p,t)
k1=0.01
k2=0.01
Dt=2000;
println(u[2])
param_cal= (1-u[2]/0.1)^0.8
du[1] = Dt - u[1]*k1;
du[2] = -1e-9;]
return du
end
println("Starting....")
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, 0.1]
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);
elseif solver_choose=="kencarp4offauto" || solver_choose=="kencarp4auto"
sol = solve(prob, abstol=tol, reltol=tol,solver, dtmax=0.1);
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")