Rodas methods do not error control the interpolation of the algebraic variables

[topolarity]

This is in essence the same problem as https://github.com/SciML/OrdinaryDiffEq.jl/issues/2045, except that differential variables are present but evolve far too slowly to provide indirect error control for algebraic variables.

using ModelingToolkit

@variables t
D = Differential(t)

@mtkmodel Foo begin
    @parameters begin
        ω₁
        ω₂
        l
    end
    @variables begin
        x(t) # differential variable
        y(t) # differential variable
        z(t) # algebraic variable
    end
    @equations begin
        D(x) ~  ω₁ * 2pi * y    # solution: x(t) = sin(2πω₁t)
        D(y) ~ -ω₁ * 2pi * x    # solution: y(t) = cos(2πω₁t)
        z^l ~ sin(ω₂ * 2pi * t) # solution: z(t) = sin(2πω₂t) when l == 1.0
    end
end

@mtkbuild foo = Foo();

Solve and plot with a variety of solvers:

using DifferentialEquations: solve, Rosenbrock23, Rodas4, FBDF

prob = ODEProblem(foo,
                  [foo.x => 1.0, foo.y => 0.0, foo.z => 0.0],
                  (0.0, 1.0),
                  [foo.ω₁ => 1.0, foo.ω₂ => 1000.0, foo.l => 1.0])

sol0 = solve(prob, Rosenbrock23())
sol1 = solve(prob, Rodas4())
sol2 = solve(prob, FBDF())

plots = []
for (name, sol) in [
    ("Rosenbrock23", sol0),
    ("Rodas4", sol1),
    ("FBDF", sol2),
]
    z = sol(sol.t, idxs=foo.z)
    plt = plot(z, layout=(2,1), legend=:none, title="$(name) (t=0..1)")
    plot!(plt[2], z[1:findfirst(sol.t .>= 2e-3)], xlim=(0,2.1e-3), legend=:none, title="$(name) (t=0..2e-3)")
    push!(plots, plt)
end
plots

results in:

In summary:

• The FBDF() solution is correct.
• The Rodas4() solution is correct at each sol.t but is heavily under-sampled
• The Rosenbrock23() is wildly incorrect, even at each sol.t

[topolarity]

Worth mentioning also that the dynamics don't have to be very much faster for this to be a noticeable issue:

Notice that the Rosenbrock23() solution still exceeds (-1, 1) significantly. Of course, interpolating this algebraic solution is a whole separate can of worms.

[gstein3m]

I think, we have two problems:

1) As stated in the book of Hairer/Wanner, for stiffly accurate methods like Rodas4, Rodas5 and variants, the numerical solution of a simple algebraic condition like g(t,z)=0 results in a simplified newton iteration. For that reason the method Rodas4 and its embedded scheme both yield nearly exact results for z(t) = sin(2pi fre t) and the stepsize control fails. In that case we should set dtmax < 0.5/freq, a good choice could be 0.1/freq.

2) Why the stepsize control of Rosenbrock23 fails, is not clear to me. Here a large difference between the methods of order 2 and 3 should occur, since the method of order 3 ist not A-stable. In a simple own implementation with reltol=abstol=1.0e-4, I indeed need 2297 steps with a max. error of 6.4e-5 for a frequency = 10. When I apply Rosenbrock23 as below, only 7 steps are needed with an error = 32.9. Thus the implementation should be checked.

using OrdinaryDiffEq, Plots

function f!(du, u, fre, t)
    du[1] = u[1] - sin(fre*2*pi*t)
end

u0 = [0.0]; tspan = (0.0, 1.0); M = zeros(1,1); fre = 10.0
fode = ODEFunction(f!, mass_matrix = M); prob = ODEProblem(fode, u0, tspan, fre)
sol = solve(prob, Rosenbrock23())
println(sol.destats)
err_max = maximum(abs.(sol[1,:] .- sin.(fre*2*pi*sol.t)))
println("err_max:",err_max)

[gstein3m]

Looking at rosenbrock_perform_step.jl, I see:

line 119/120

        @.. broadcast=false vectmp=ifelse(cache.algebraic_vars,
            dto6 * (veck₁ - 2 * veck₂ + veck₃))

and line 241/242

          vectmp[i] = ifelse(cache.algebraic_vars[i], false,
                dto6 * (veck₁[i] - 2 * veck₂[i] + veck₃[i]))

but line 389

utilde = dto6 * (k₁ - 2 * k₂ + k₃)

Do I understand correctly that algebraic variables are excluded from the error estimation? If so, the behaviour of Rosenbrock23 is clear.

I can only assume that this is because R(\infty)>1 applies to the stability function of the method of order p=3 and therefore it is not convergent for DAEs.

Maybe we should try to improve the method, I will start working on this topic.

[ChrisRackauckas]

Indeed, that behavior comes from:

• https://github.com/SciML/OrdinaryDiffEq.jl/issues/1606
• https://github.com/SciML/OrdinaryDiffEq.jl/pull/1644

And you're right. The underlying problem is that the 3rd order method used in the error estimation is not stable, and therefore it's not convergent in the algebraic variables and you get a wonky error estimate. This means that the only realistic way to handle them given the mathematical structure of the method is to exclude the algebraic variables from the error estimate, since otherwise you just always get dt<dtmin failures from the error estimator not actually converging to zero as dt -> 0. This means it's not an implementation issue with Rosenbrock23 but rather an algorithmic one. It would be nice to have a corrected 3rd order embedded form that is stiffly accurate, but I presume doing that and keeping it a RosW method at the same time is difficult.

I would be interested if there was a second order Rosenbrock that was adaptive and more stable though, but non RosW. Rosenbrock23 is quite a workhorse method in the ODE case in that space but indeed the instability of its error estimator sometimes makes it wonky (even for ODEs), but Shampine must have made a compromise there for a reason.

[ChrisRackauckas]

We should probably throw a warning in general on Rosenbrock23 for DAEs given that behavior then? It might be a sensible thing to do. https://github.com/SciML/OrdinaryDiffEq.jl/pull/2050 we could change that to just always warn on DAEs. It's still a nice method for a lot of situations though.

As stated in the book of Hairer/Wanner, for stiffly accurate methods like Rodas4, Rodas5 and variants, the numerical solution of a simple algebraic condition like g(t,z)=0 results in a simplified newton iteration. For that reason the method Rodas4 and its embedded scheme both yield nearly exact results for z(t) = sin(2pi fre t) and the stepsize control fails. In that case we should set dtmax < 0.5/freq, a good choice could be 0.1/freq.

If you're looking for the higher value topic, I think this one is it. I don't know of a single Rosenbrock that does residual adaptivity. A good example to look at is:

https://www.sciencedirect.com/science/article/abs/pii/S0168927404001187 (1-s2.0-S0168927404001187-main.pdf)

In code terms this is:

https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/src/perform_step/fixed_timestep_perform_step.jl#L295-L321

If we can just derive the maxima and minima points on the interpolation error and do something similar, this would be a good option to enable on DAEs. It does add extra f evaluations, but if you're constructing the Jacobian already then 3 or 4 extra f evaluations isn't the worst thing in the world. We can just alias it to another algorithm so that it's clear it has very different stepping behavior, i.e. Rodas5PR or something for residual error estimator, but internally it's just a switch in that part of the code.

The derivation I think isn't too hard? It's just finding the critical points of the interpolation and hardcoding the error estimate to evaluate at each critical point. This would also make Rosenbrocks in DDEs better, which is a nice extra touch as well.

[gstein3m]

"We should probably throw a warning in general on Rosenbrock23 for DAEs given that behavior then? "

• the user should at least critically examine the result

Regarding an alternative to Rosenbrock23: It is possible to derive a stiffly accurate A-stable Rodas23W with

• s=3, order p=2 for ODEs and DAEs of index-1 and 2, W p=2 for ODEs, W p=2 even for DAEs of index-1 when Jacobian error is O(h)
• embedded method: also stiffly accurate and A-stable, s=4 with no further f evaluation, p=3 for ODEs and index-1 DAEs
• dense output of order p=3, for ODEs and index-1 DAEs with s=5 but no further f evaluation I have to do further accuracy considerations and benchmarks. In particular I am interested in a strategy, to switch between the methods of order p=2 (reusing Jacobian and lu-decomposition for same stepsize as before ) and p=3 (when Jacobian is new).

Regarding residual adaptivity for rosenbrock methods: I will study this ...

[ChrisRackauckas]

It is possible to derive a stiffly accurate A-stable Rodas23W with

Wow, on paper that sounds strictly better than Rosenbrock23.

[gstein3m]

When implementing the new method, I came across a bug in Rodas3. It produces big errors for the example above. The stage 4 is missing in the implementation.

@muladd function perform_step!(integrator, cache::Rosenbrock34Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator @unpack du, du1, du2, fsalfirst, fsallast, k1, k2, k3, k4, dT, J, W, uf, tf, linsolve_tmp, jac_config, atmp, weight = cache @unpack a21, a31, a32, C21, C31, C32, C41, C42, C43, b1, b2, b3, b4, btilde1, btilde2, btilde3, btilde4, gamma, c2, c3, d1, d2, d3, d4 = cache.tab

must be

@muladd function perform_step!(integrator, cache::Rosenbrock34Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator @unpack du, du1, du2, fsalfirst, fsallast, k1, k2, k3, k4, dT, J, W, uf, tf, linsolve_tmp, jac_config, atmp, weight = cache @unpack a21, a31, a32, a41, a42, a43, C21, C31, C32, C41, C42, C43, b1, b2, b3, b4, btilde1, btilde2, btilde3, btilde4, gamma, c2, c3, d1, d2, d3, d4 = cache.tab

and in line 844

@.. broadcast=false u=uprev + a41 * k1 + a42 * k2 + a43 * k3 f(du, u, p, t + dt) #-- c4 = 1 integrator.stats.nf += 1

must included. This is running well. I can't say whether any changes to other functions are necessary

[ChrisRackauckas]

Interesting. Thanks for handling that. It's surprising that dropping that stage passed the 3rd order convergence tests, but I guess it just fell onto a different 3rd or method which is rare but possible.

[ChrisRackauckas]

https://github.com/SciML/OrdinaryDiffEq.jl/pull/2079 handles Rosenbrock23, so I'm going to update this title to just be about the remaining issue with interpolations.

[gstein3m]

In methods Rodas4, 42, 4P, 4P2, 5, 5P we could include an additional error test for ||M y' - f(t,y)|| / abstol < 1. Here y is the interpolation and we could check at time t = t_n + h/2. This test would cost one additional function evaluation. This gave good results in a preliminary program. However, it would be nice if there were an additional option to switch this test on or off. I could implement the test in rosenbrock_perform_step, but I don't know how to integrate an additional option.

[gstein3m]

See recent Preprint and tests in https://github.com/hbrs-cse/RosenbrockMethods

A pull request concerning Rodas5Pr and Rodas5Pe will come soon.