Robertson problem

[devmotion]

I started playing around a bit more with stiff DDEs, and think www.unige.ch/~hairer/radar5-v2.1.tar and accompanying publications contain some nice examples (which should be added to DiffEqProblemLibrary as well, I guess).

However, my first tests show that there are some problems with stiff problems. The first problem is a modified Robertson problem, taken from "The numerical integration of stiff delay differential equations" by Guglielmi; I only changed the integration interval from [0.0, 1e9] to [0.0, 1e5] since otherwise maxiters has to be increased (maybe the number of iterations could be reduced in general by changing some defaults?):

using DelayDiffEq

function f_dde_rober1(du, u, h, p, t)
    # compute repeating terms
    x = 4e-2 * u[1]
    y = 1e4 * u[2] * u[3]
    z = 3e7 * (1 - h(p, t - 1; idxs = 2)) * u[2]^2

    # compute derivatives
    du[1] = y - x
    du[2] = x - y - z
    du[3] = z

    nothing
end
# simplified history function, hence `initial_order` has to be set
h_dde_rober(p, t; idxs = 2) = 0.0

prob_dde_rober1 = DDEProblem(f_dde_rober1, [1.0, 0.0, 0.0], h_dde_rober,
                             (0.0, 1e5);
                             constant_lags = [1])

Using Tsit5 integrations is stopped at around t = 627 since the maximal number of iterations is reached, and the solution is wrong as the following plot shows:

sol = solve(prob_dde_rober1, MethodOfSteps(Tsit5()); initial_order = 1)

using Plots
plotly()
plot(sol)

Using Rosenbrock23 a solution is computed on the whole integration interval:

sol = solve(prob_dde_rober1, MethodOfSteps(Rosenbrock23()); initial_order = 1)

plot(sol)

This solution shows the expected dynamics; moreover, it indicates that the stiff solvers work in principle and are quite performant since according to Guglielmi "RADAR5 is the only code able to correctly integrate the above very stiff problem" out of the solvers he considered.

Interestingly, both in Guglielmi and Hairer's "Implementing Radau IIA methods for stiff delay differential equations" and as example problem in the RADAR5 package a different Robertson problem is mentioned (here again on a reduced integration interval):

function f_dde_rober2(du, u, h, p, t)
    # compute repeating terms
    x = 4e-2 * u[1]
    y = 1e4 * h(p, t - 1e-2; idxs = 2) * u[3]
    z = 3e7 * u[2] * u[2]

    # compute derivatives
    du[1] = y - x
    du[2] = x - y - z
    du[3] = z

    nothing
end

prob_dde_rober2 = DDEProblem(f_dde_rober2, [1.0, 0.0, 0.0], h_dde_rober,
                             (0.0, 1e5);
                             constant_lags = [1e-2])

For this problem, however, Rosenbrock23 computes a solution only on [0.0, 6.14] since the time steps become too small:

sol = solve(prob_dde_rober2, MethodOfSteps(Rosenbrock23()); initial_order = 1)

plot(sol)

Similar dynamics can be observed with Rodas5, although in this case integration stops because of instability (i.e. NaN values):

sol = solve(prob_dde_rober2, MethodOfSteps(Rodas5()); initial_order = 1)

plot(sol)

Using absolute and relative tolerances that are mentioned in Guglielmi and Hairer's publication I can provoke the same explosions as before at a slightly later time point:

sol = solve(prob_dde_rober2, MethodOfSteps(Rodas5()); initial_order = 1,
                                           reltol = 1e-6, abstol = 1e-10)

plot(sol)

According to Guglielmi and Hairer, with a delay of 3e-2 the second component of this problem oscillates and then explodes at around t = 16.8.

I managed to compile and run RADAR5 and the contained Robertson example on my computer; RADAR5 successfully computed a solution on the interval [0.0, 1e10]. Even though I have no experience with Fortran I studied the implementation of the problem to check whether the implementation corresponds to the stated DDE; hereby I noticed that Guglielmi and Hairer apparently used slightly different default parameters in their implementation (tau = 1e-3, reltol = 1e-9, abstol = 1e-14, and an integration interval of [0.0, 1e11]) but I could still compute a correct solution after changing these values to the ones stated above.

I do not know why DelayDiffEq does not compute a correct solution for the second problem. I'm also surprised that Rodas5 with default settings returns NaN values; I assumed this could (and should) not happen since then the time step should be decreased... I also played around with the number of fixed point iterations; setting force_dtmin = true in the last example resulted in NaN values as well. Are there any settings I could tune? Do you have any ideas why RADAR5 can solve this problem whereas I could not compute a solution with DelayDiffEq?

[ChrisRackauckas]

Yes, I noticed this awhile ago but forgot to document it.

There's a few things that can be happening. First of all, we don't actually know that the stability of the Rosenbrock methods are like when applied to DDEs. I think there's some concepts like Q-stability and nothing is fully Q-stable, so it would be an interesting analysis to see what kind of stability we actually get. It sounds like it's quite good for most "stiff DDEs" (from my tests as well), but you can break it if you try hard enough.

Another thing to note is that our extrapolation method could reduce stability. Functional iteration is not used in stiff ODE solvers because of this property. Anderson acceleration can help, but Newton's method is much more stable. With our generic method of steps we cannot take into account a full Newton on every step including the extrapolation, so that is one issue we do have that RADAR5 doesn't. I'm not sure how much this matters though, and the fact that RADAR5 is solving such a large nonlinear system is a huge cost burden so there is a tradeoff.

Another thing we might want to look into is the extrapolation over discontinuities. When extrapolating over a discontinuity, the calculation becomes incorrect because of that discontinuity in the interval effecting the interpolation result. I think one of the papers we've looked at before suggested doing something like:

1. Solve a step
2. If you find a discontinuity in the interval via event handling, shorten the step
3. If not in the interval, take the step. If in the interval, shorten a bit more.
4. Keep shortening till you hone in on that discontinuity to be near the edge of the step. Then take that step.

The reason here is that then less of the estimate uses the values past the discontinuity so it gives a more accurate answer. I'm not sure if that could be affecting us here, but it would be good to open an issue to investigate it.

These last two points can be eliminated with constrained stepping though. If constrained stepping still is unstable, then it is likely just a property of the Rosenbrock integrator. But the interesting thing is:

Using Tsit5 integrations is stopped at around t = 627 since the maximal number of iterations is reached, and the solution is wrong as the following plot shows:

No it's not. That looks correct. You can even check it with:

sol = solve(prob_dde_rober1, MethodOfSteps(Rosenbrock23()); abstol=1e-8,reltol=1e-8,initial_order = 1)
sol2 = solve(prob_dde_rober1, MethodOfSteps(Tsit5()); abstol=1e-8,reltol=1e-8,initial_order = 1)
sol2[end] - sol(sol2.t[end])

So it's just exiting early. The Runge-Kutta methods at extremely high accuracy (with bigfloats) should be able to calculate the solution correctly and not have any stability issues just by using a small dt (and no numerical instability issues because of bigfloats). They might take forever to run, but they will be correct.
So I am curious how Hairer and Gugliemi know that the value should be t = 16.8 since I can't seem to get an RK method to say it's above 9. We should run some really detailed solves:

prob_dde_rober2_big = DDEProblem(f_dde_rober2, big.([1.0, 0.0, 0.0]), h_dde_rober,
                             big.((0.0, 1e5));
                             constant_lags = [1e-2])
sol = solve(prob_dde_rober2_big, MethodOfSteps(Tsit5()); abstol=1e-20,reltol=1e-20,initial_order = 1)
Float64(sol.t[end])

sol = solve(prob_dde_rober2_big,
            MethodOfSteps(Rosenbrock23(linsolve=LinSolveFactorize(qrfact!)));
            abstol=1e-12,reltol=1e-12,initial_order = 1)
Float64(sol.t[end])

and get a definitive answer and double check that t=16.8 value. Also, I think you're using a slightly different delay than them?

As a side note, I know that maybe a 6 months from now @YingboMa will want to do fully implicit RK methods in OrdinaryDiffEq, in which case we should make sure we get a full RADAR5 implementation to see for ourselves how well it can do.

[ChrisRackauckas]

We should also check with autodiff=false.

[devmotion]

No it's not. That looks correct. You can even check it with:

Ah, of course. Don't know why I didn't notice that Tsit5 just exits early.

[devmotion]

Another thing we might want to look into is the extrapolation over discontinuities. When extrapolating over a discontinuity, the calculation becomes incorrect because of that discontinuity in the interval effecting the interpolation result. I think one of the papers we've looked at before suggested doing something like:

Interestingly, RADAR5 only estimates discontinuity points if

• the Newton process does not converge
• the estimated error is not under a given tolerance
• the error increases (from a step to the subsequent) over a prescribed threshold

[ChrisRackauckas]

On unconstrained meshes good results are obtained only when the integrator is equipped by a smothness indicator for the control of the continuous solution

Interesting. To fix the issue, they:

1. Modify the interpolation to only use values before
2. Incorporate the state-dependent delay terms into the Jacobian to get better convergence and discontinuity checks if applicable
3. Use a residual control error estimator

[devmotion]

So I am curious how Hairer and Gugliemi know that the value should be t = 16.8 since I can't seem to get an RK method to say it's above 9.

As I mentioned above, Hairer and Guglielmi computed this value for a delay of 3e-2, and not for problem prob_dde_rober2; they even provide a plot of the problem with delay 1e-2 and one of the problem with delay 3e-2. However, I can only get to a time point of around 10.29 with Tsit5 and around 10.31 with Rosenbrock23 before dt becomes too small:

function f_dde_rober3(du, u, h, p, t)
    # compute repeating terms
    x = 4e-2 * u[1]
    y = 1e4 * h(p, t - 3e-2; idxs = 2) * u[3]
    z = 3e7 * u[2] * u[2]
       
    # compute derivatives
    du[1] = y - x
    du[2] = x - y - z
    du[3] = z
       
    nothing
end

prob_dde_rober3_big = DDEProblem(f_dde_rober3, big.([1.0, 0.0, 0.0]), h_dde_rober,
                                 big.((0.0, 1e5));
                                 constant_lags = [3e-2])

sol1 = solve(prob_dde_rober3_big, MethodOfSteps(Tsit5()); initial_order = 1,
             abstol = 1e-20, reltol = 1e-20)
Float64(sol1.t[end])

sol2 = solve(prob_dde_rober3_big,
             MethodOfSteps(Rosenbrock23(linsolve=LinSolveFactorize(qrfact!)));
             initial_order = 1, abstol = 1e-12, reltol = 1e-12)
Float64(sol2.t[end])

[devmotion]

These last two points can be eliminated with constrained stepping though. If constrained stepping still is unstable, then it is likely just a property of the Rosenbrock integrator.

I could reproduce the instability for constrained stepping as well:

sol = solve(prob_dde_rober2, MethodOfSteps(Rosenbrock23(); constrained=true);
            initial_order = 1)

plot(sol)

sol = solve(prob_dde_rober2, MethodOfSteps(Rodas5(); constrained=true); initial_order = 1)

plot(sol)

[devmotion]

We should also check with autodiff=false.

The algorithms are unstable with finite differencing as well:

 sol = solve(prob_dde_rober2, MethodOfSteps(Rosenbrock23(; autodiff = false));
             initial_order = 1)

sol = solve(prob_dde_rober2, MethodOfSteps(Rodas5(; autodiff = false)); initial_order = 1)

And finite differencing with constrained stepping:

sol = solve(prob_dde_rober2, MethodOfSteps(Rosenbrock23(; autodiff = false); constrained=true);
            initial_order = 1)

sol = solve(prob_dde_rober2, MethodOfSteps(Rodas5(; autodiff = false); constrained=true);
            initial_order = 1)

[ChrisRackauckas]

I'm running a really refined RK right now to double check some things. Can you run an SDIRK like TRBDF2 and see if that gives the same results as Rosenbrock?

[devmotion]

Can you run an SDIRK like TRBDF2 and see if that gives the same results as Rosenbrock?

Unfortunately, using TRBDF2 does not help:

 sol = solve(prob_dde_rober2, MethodOfSteps(TRBDF2());
             initial_order = 1)

 sol = solve(prob_dde_rober2, MethodOfSteps(TRBDF2(); constrained = true));
             initial_order = 1)

 sol = solve(prob_dde_rober2, MethodOfSteps(TRBDF2(; autodiff = false));
             initial_order = 1)

 sol = solve(prob_dde_rober2, MethodOfSteps(TRBDF2(; autodiff = false); constrained = true);
             initial_order = 1)

[ChrisRackauckas]

In theory it should be stable.

http://qktg.shnu.edu.cn/zrb/shsfqkszrb/ch/reader/create_pdf.aspx?file_no=201402001&year_id=2014&quarter_id=2&falg=1
https://www.sciencedirect.com/science/article/pii/S0096300303002893

[ChrisRackauckas]

sol1 = solve(prob_dde_rober3_big, MethodOfSteps(Tsit5()); initial_order = 1,
             abstol = 1e-28, reltol = 1e-28, maxiters = 1e7,
             unstable_check = (dt,u,p,t) -> abs(u[3])>1e2
             )
Float64(maximum(diff(sol1.t))) # 3.750184086763931e-5
Float64(sol1.t[end]) 11.25522208749271

and exited due to stability. So there likely is an instability right around there. Hmm.

[devmotion]

Hmm it's really strange that we can not verify the instability around 16.8.

[ChrisRackauckas]

Maybe it could be due to not using a residual control?

[devmotion]

I ran RADAR5 on this problem with a delay of 0.03 (which should explode according to Hairer and Guglielmi) and could compute a reasonable solution on the whole interval [0.0, 1e11]. So I assume the publication might be incorrect and the problem might not be unstable for this delay. However, using a delay of 0.3 RADAR5 exits at around 15.5...

[ChrisRackauckas]

I wonder if we can save out values of the history function and plot that just to make sure those are following the true trajectory. I think there may be a bug in here. Check OrdinaryDiffEq:

using OrdinaryDiffEq
function f_rober2(du, u, p, t)
    # compute repeating terms
    x = 4e-2 * u[1]
    y = 1e4 * u[2] * u[3]
    z = 3e7 * u[2] * u[2]

    # compute derivatives
    du[1] = y - x
    du[2] = x - y - z
    du[3] = z

    nothing
end
prob_ode_rober2 = ODEProblem(f_rober2, [1.0, 0.0, 0.0], (0.0, 1e5))
sol = solve(prob_ode_rober2,Rosenbrock23())

which works. However,

using DelayDiffEq
function f_dde_rober2(du, u, h, p, t)
    # compute repeating terms
    x = 4e-2 * u[1]
    y = 1e4 * h(p, t-1e-16)[2] * u[3]
    z = 3e7 * u[2] * u[2]

    # compute derivatives
    du[1] = y - x
    du[2] = x - y - z
    du[3] = z

    nothing
end
h_dde_rober(p, t) = zeros(3)
prob_dde_rober2 = DDEProblem(f_dde_rober2, [1.0, 0.0, 0.0], h_dde_rober,
                             (0.0, 1e5);
                             constant_lags = [1e-2])
sol = solve(prob_dde_rober2,MethodOfSteps(Rosenbrock23()))

fails with maxiters. Checking:

function f_dde_rober4(du, u, h, p, t)
    # compute repeating terms
    x = 4e-2 * u[1]
    y = 1e4 * h(p, t-1e-16)[2] * u[3]
    z = 3e7 * u[2] * u[2]

    @show h(p, t-1e-16)[2] - u[2]
    # compute derivatives
    du[1] = y - x
    du[2] = x - y - z
    du[3] = z

    nothing
end
h_dde_rober(p, t) = zeros(3)
prob_dde_rober2 = DDEProblem(f_dde_rober4, [1.0, 0.0, 0.0], h_dde_rober,
                             (0.0, 1e5))
sol = solve(prob_dde_rober2,MethodOfSteps(Rosenbrock23(),constrained=true),maxiters=1000)

gives things like:

(h(p, t - 1.0e-16))[2] - u[2] = 2.2215283979943968e-5
(h(p, t - 1.0e-16))[2] - u[2] = 6.0901529403518086e-5
(h(p, t - 1.0e-16))[2] - u[2] = 0.0
(h(p, t - 1.0e-16))[2] - u[2] = 5.041613513194181e-7
(h(p, t - 1.0e-16))[2] - u[2] = -1.360253148542067e-5
(h(p, t - 1.0e-16))[2] - u[2] = 0.0
(h(p, t - 1.0e-16))[2] - u[2] = -9.044269046495332e-7
(h(p, t - 1.0e-16))[2] - u[2] = -8.82513835403601e-7

while on the ODE solution, doing stuff like sol(sol.t[30]-1e-16) - sol[30] is almost exactly zero (and I know this is almost always the case. The numerical issues should cause it to return exactly the value it's close to rather than oscillate).

So it sounds to me like there may be a problem with the history values.

[ChrisRackauckas]

Now he's something that's interesting:

function f_dde_rober4(du, u, h, p, t)
    # compute repeating terms
    x = 4e-2 * u[1]
    y = 1e4 * h(p, t-1e-16)[2] * u[3]
    z = 3e7 * u[2] * u[2]

    @show h(p, t-1e-16)[2] - u[2]
    # compute derivatives
    du[1] = y - x
    du[2] = x - y - z
    du[3] = z

    nothing
end
h_dde_rober(p, t) = zeros(3)
prob_dde_rober2 = DDEProblem(f_dde_rober4, [1.0, 0.0, 0.0], h_dde_rober,
                             (0.0, 1e5))
sol = solve(prob_dde_rober2,MethodOfSteps(Rodas5(autodiff=false)),dt=1e-2)

(h(p, t - 1.0e-16))[2] - u[2] = -6.0554544523933395e-6
(h(p, t - 1.0e-16))[2] - u[2] = 6.0554544523933395e-6
(h(p, t - 1.0e-16))[2] - u[2] = 0.0
(h(p, t - 1.0e-16))[2] - u[2] = 0.0
(h(p, t - 1.0e-16))[2] - u[2] = -0.00015198844887788534
(h(p, t - 1.0e-16))[2] - u[2] = 0.001216570117719098
(h(p, t - 1.0e-16))[2] - u[2] = -0.07309512770559262
(h(p, t - 1.0e-16))[2] - u[2] = 87.69690577013398
(h(p, t - 1.0e-16))[2] - u[2] = 1.0578873888263168e10
(h(p, t - 1.0e-16))[2] - u[2] = 6.379016646395286e24
(h(p, t - 1.0e-16))[2] - u[2] = 2.3194356423743255e54
(h(p, t - 1.0e-16))[2] - u[2] = 3.0664755684963483e113
(h(p, t - 1.0e-16))[2] - u[2] = 0.0
(h(p, t - 1.0e-16))[2] - u[2] = 2.525568517554205e231
(h(p, t - 1.0e-16))[2] - u[2] = 2.5976895522572776e231
(h(p, t - 1.0e-16))[2] - u[2] = 3.348148356188627e231
(h(p, t - 1.0e-16))[2] - u[2] = NaN
(h(p, t - 1.0e-16))[2] - u[2] = NaN
(h(p, t - 1.0e-16))[2] - u[2] = NaN
(h(p, t - 1.0e-16))[2] - u[2] = NaN

it goes out of control?

[ChrisRackauckas]

I think it's highlighting that we're not handling step rejections properly? It seems we advance_ode_integrator!(integrator) on every perform_step!, but don't pull back if the step failed or if integrator.EEst>1.

[devmotion]

That might very well be the case... I already assumed something along those lines when I got NaN values in the first runs above.

[ChrisRackauckas]

That test case was bad because it requires using the isout handling. However, I am narrowing down a bug which is requiring isout when dt is constrained to be smaller than the delay size, and that is for sure an issue.

[ChrisRackauckas]

On the bug branch, it errors if history extrapolation is ever used. I run the problem:

using DelayDiffEq
function f_dde_rober4(du, u, h, p, t)
    # compute repeating terms
    x = 4e-2 * u[1]
    y = 1e4 * h(p, t-1e-4)[2] * u[3]
    z = 3e7 * u[2] * u[2]
    # compute derivatives
    du[1] = y - x
    du[2] = x - y - z
    du[3] = z

    nothing
end
h_dde_rober(p, t) = zeros(3)
prob_dde_rober2 = DDEProblem(f_dde_rober4, [1.0, 0.0, 0.0], h_dde_rober,
                             (0.0, 1e5))
sol = solve(prob_dde_rober2,MethodOfSteps(Rodas5(autodiff=false)),dt=5e-5,dtmax=5e-5)

and after quite awhile, I get that the history function goes beyond f.sol.t[end] by ~1e-11. However, that should never happen since dtmax < 1e-2 which is the delay size, so there's definitely a bug there and I haven't been able to narrow down how that could happen.

"here2" = "here2"
integrator.dt = 5.0e-5
integrator.sol.t[end] = 16.513999999986442
integrator.integrator.sol.t[end] = 16.513999999986442
"out of perform step" = "out of perform step"
"here2" = "here2"
integrator.dt = 5.0e-5
integrator.sol.t[end] = 16.514049999986444
integrator.integrator.sol.t[end] = 16.514049999986444
f.sol.t[end] = 16.514049999986444
t = 16.514050000064042

[ChrisRackauckas]

Nevermind. That "bug" was just because finite differencing went too far.