Best practices in dealing avoiding with precision/floating round of problems when solving ODE's? Thank you

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!


using DifferentialEquations
using Sundials

function  ode_functionminie(du,u, p,t)


param_cal= (1-u[2]/0.1)^0.8

du[1] = Dt - u[1]*k1; 
du[2] = -1e-9;]

return du



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, 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);
    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))))


Either way it’s the same thing if you want to reject the step.

1 Like

I think I understand , thank you , Chris!