Issues with radau solvers

[vsunye]

Hello,

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

[ChrisRackauckas]

I am trying to solve a linear ODE with constant coefficients whose solution I know for sure is stable.

If you're solving a linear ODE, then you don't need to use general linear methods. It can be much faster to use things which specialize on the matrix exponential. See https://docs.sciml.ai/DiffEqDocs/stable/solvers/nonautonomous_linear_ode/#Exponential-Methods-for-Linear-and-Affine-Problems which shows how to define an operator and use that for LinearExponential if you have a lazy operator, or https://github.com/SciML/ExponentialUtilities.jl for faster matrix exponentials. Or just use exp(A*t) from Julia base which has decent performance.

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.)

Note that none of those are 0. It's actually impossible to represent e^x for almost any value except 0 to infinite digits of accuracy. A 64-bit float on a computer only can represent numbers to 16 digits of accuracy. Using exp(A*t) is 1ulp which means 1 unit in last place error, so accuracy around 15 digits.

If attempted to solve above problem with absolute tolerance , then both solvers fail because the timestep at
is too small. However, by setting the absolute tolerance to the smallest positve Float64 (i. e., nextfloat(0.0)), both solvers succeed.

The standard numerical methods are not expected to work at 0 tolerance. They would require a dt = 0 in order to do so, which would take infinitely many steps. All methods are only numerical approximations.

If you really need fast and stable calculations of linear ODE solutions, just use the matrix exponential tools.