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 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/0.1)^0.8. I have not tried this , but if param_cal was= (0.1- u)^0.8 , would it have also given similar error IC of u is 0.1, or is it the division that creates this issue.
Hope you find the message clear, thank you again!
using DifferentialEquations using Sundials function ode_functionminie(du,u, p,t) k1=0.01 k2=0.01 Dt=2000; println(u) param_cal= (1-u/0.1)^0.8 du = Dt - u*k1; du = -1e-9;] return du end println("Starting....") 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, 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), length(MATraw)))) println("Done")