New to Julia and to DifferentialEquations.jl, so please be forgiving.
I am trying to solve a linear ODE with constant coefficients whose solution I know for sure is stable.
I would like to use one of the radau solvers (RadauIIA3, RadauIIA5). Relative tolerance has to be $0$. (I need to work with the unscaled local error estimates.)
After experiencing issues with failures caused by too small a timestep at $t = 0$ and trying to ensure the issue was not with my specification of the problem, I have turned to $$\frac{du}{dt} = - \lambda\cdot u, u(0)=1,$$ which solution is $$u(t)=e ^ {- \lambda\cdot t}.$$
If attempted to solve above problem with absolute tolerance $0$, then both solvers fail because the timestep at $t=0$ is too small. However, by setting the absolute tolerance to the smallest positve Float64 (i. e., nextfloat(0.0)), both solvers succeed.
The code at the end of this message will allow to exercise the failure by uncommenting the line #alg = RadauIIA3().
Any help would be greatly appreciated.
Regards,
Víctor Suñé
PS Minimal working example
using DifferentialEquations
function f(u,p,t)
return - p[1] * u
end
lambda = 1 / 100000.0
tspan = (0.0, 1 / lambda) # Time span
rTol = 0.0 # Relative tolerance for the ODE solver
rTol = nextfloat(0.0) # Uncomment if a radau solver is to be used
aTol = 1E-9 # Absolute one
alg = nothing
alg = RadauIIA3()
odef = ODEFunction{false}(f)
problem = ODEProblem(odef, 1.0, tspan, (lambda))
res = solve(problem, alg, reltol = rTol, abstol = aTol)
Hello,
New to
Julia
and toDifferentialEquations.jl
, so please be forgiving.I am trying to solve a linear ODE with constant coefficients whose solution I know for sure is stable.
I would like to use one of the radau solvers (RadauIIA3, RadauIIA5). Relative tolerance has to be $0$. (I need to work with the unscaled local error estimates.)
After experiencing issues with failures caused by too small a timestep at $t = 0$ and trying to ensure the issue was not with my specification of the problem, I have turned to $$\frac{du}{dt} = - \lambda\cdot u, u(0)=1,$$ which solution is $$u(t)=e ^ {- \lambda\cdot t}.$$
If attempted to solve above problem with absolute tolerance $0$, then both solvers fail because the timestep at $t=0$ is too small. However, by setting the absolute tolerance to the smallest positve
Float64
(i. e.,nextfloat(0.0)
), both solvers succeed.The code at the end of this message will allow to exercise the failure by uncommenting the line
#alg = RadauIIA3()
.Any help would be greatly appreciated.
Regards, Víctor Suñé
PS Minimal working example
using DifferentialEquations
function f(u,p,t) return - p[1] * u end lambda = 1 / 100000.0
tspan = (0.0, 1 / lambda) # Time span rTol = 0.0 # Relative tolerance for the ODE solver
rTol = nextfloat(0.0) # Uncomment if a radau solver is to be used
aTol = 1E-9 # Absolute one alg = nothing
alg = RadauIIA3()
odef = ODEFunction{false}(f) problem = ODEProblem(odef, 1.0, tspan, (lambda)) res = solve(problem, alg, reltol = rTol, abstol = aTol)