Rodas5P gives inaccurate results when integrating triangle wave

[topolarity]

Describe the bug 🐞

The example here is essentially x' ~ f(t) (in semi-explicit DAE form) where f(t) is a continuous, piecewise-linear function with a strongly discontinuous derivative, for which Rodas5P gives poor results.

This is arguably a duplicate of https://github.com/SciML/OrdinaryDiffEq.jl/issues/2054, except:

• the system here has no "fast" / "slow" dynamics
• the error is very bad for both algebraic and differential variables
• the pure ODE form exhibits the same problem

Anyway, here's the MWE:

Minimal Reproducible Example 👇

julia> using OrdinaryDiffEq, LinearAlgebra
julia> triangle(t, T) = 4abs(0.5 - mod(t/T - 0.25, 1)) - 1.0
julia> f(u, p, t) = [u[2], u[2] - triangle(t, 1.0)]
julia> mass_matrix = Diagonal([1.0,0.0])
julia> u0 = [10.0, 0.0]
julia> prob = ODEProblem(ODEFunction(f; mass_matrix), u0, (0., 2π), 1)
julia> sol1 = solve(prob, Rodas5P());
julia> sol2 = solve(prob, ROS34PRw());

julia> t′ = 0:1e-3:2π;
julia> plot(sol1; idxs=2, label="Rodas5P u[2]")
julia> plot!(sol2; idxs=2, label="ROS34PRw u[2]")
julia> plot!(t′, triangle.(t′, 1.0), label="analytic u[2]", color=:gray, linealpha=0.25)

julia> function parabolic_wave(t, T)
           q = (T / 2) * (mod(2t/T-0.5, 1) - 0.5)^2
           return (mod(t/T - 0.25, 1) > 0.5) ? q : (T/4 - q)
       end
julia> plot(sol1; idxs=1, label="Rodas5P u[1]")
julia> plot!(sol2; idxs=1, label="ROS34PRw u[1]")
julia> plot!(t′, u0[1] .+ parabolic_wave.(t′, 1.0), label="analytic u[1]", color=:gray, linealpha=0.5)

The error gets worse as u0 increases, since u[1] is less affected by reltol.

Rodas5P takes only 21 timesteps but seems over-optimistic about the quality of its solution:

For reference, the new ROS34PRw solver does much better (with 103 timesteps):

A few other notes:

• Adding callbacks at the nonlinearity points is a good workaround when available (Rodas5P does quite well in that case)
• The algebraic part of this problem is non-essential (the ODE version x' ~ triangle(t) performs about the same)
• You can run into the same problem solving an infinitely-smooth ODE if the u0 offset is made large enough
  • e.g. x' ~ sin(t), u(0) = 100.0
  • In contrast to this MWE, Rodas5P arguably achieves its tolerances for this system (for the differential variables), but it is still greatly out-performed by ROS34PRw

Environment:

• Output of using Pkg; Pkg.status()

(tmp) pkg> status
Status `~/repos/tmp/tmp/Project.toml`
  [459566f4] DiffEqCallbacks v2.36.1
  [1dea7af3] OrdinaryDiffEq v6.70.0
  [c3572dad] Sundials v4.23.2

[gstein3m]

The theory of a higher order method assumes a sufficiently smooth solution. Here the solution is not continuously differentiable, which can lead to problems. Therefore, lower order methods may provide equally good or better solutions. However, a sufficiently accurate solution can be generated by adjusting reltol, e.g. reltol = 1.0e-7. An error check of the interpolation would also be helpful, see issue #2054 and new method Rodas3P.

[ChrisRackauckas]

This is just a repeat.

[topolarity]

I'm fine with that conclusion - but we should perhaps broaden the scope of #2054 to make it less about the algebraic variables.

The same problem happens here for the fully-explicit ODE x' ~ triangle(t)

[ChrisRackauckas]

It's going to happen on any exactly integrated equation, so any equation that's a polynomial in t less than the power of the order will also apply for the same reason.

[topolarity]

Right, anything that appears to have smooth derivatives at low-orders but whose global behavior is not well-determined by those derivatives (e.g. piecewise functions and also infinitely smooth but non-analytic functions like f(t) = t <= 0 ? 0.0 : exp(-1/t))

Is there a name for the kind of error control that seems to prevent this problem in ROS34PRw?