 # 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 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!

Christian

``````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")

``````

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

1 Like

I think I understand , thank you , Chris!