Benchmarks for DDEs

[ChrisRackauckas]

We should start getting some benchmarks together for DDEs. There probably aren't many good ones with analytical solutions, so the appxsol stuff will have to be used on the simplest examples. Show the difference between constrained and unconstrained, the the difference that order makes in DDEs. Also see the difference between the OwrenZen optimized interpolant methods and the classic optimized RK methods to see if it really makes a difference.

@devmotion do you have any model equations to consider?

[devmotion]

In my master thesis I'm mostly working with the model (and some variants of it) in https://link.springer.com/article/10.1007%2Fs00216-014-8063-6. However, it contains a bunch of parameters, only one variable with delay, one variable that tends to oscillate betweeen positive and negative values although mathematically it has to be positive all the time, and some variables that have to be positive in the original model to ensure a real solution (there are some arbitrary real exponents). So I guess this might not be the best suited model for benchmarks ;)

However, in http://www.radford.edu/%7Ethompson/RP/dklag6.pdf which you posted in the issue about generic state-dependent delay handling the authors benchmark different versions of their program with three different models. The first two of them are state-dependent delay models (so not suited yet), but the third example, especially the integral-free version on p. 327, could be used. They had to run their program twice since after a certain time point D the equations change from an ODE to a DDE system but with DelayDiffEq we can compute the solution in one run.

[ChrisRackauckas]

although mathematically it has to be positive all the time

Is this a numerical issue?

So I guess this might not be the best suited model for benchmarks ;)

That sounds like a perfect benchmark. ROBER and VanDerPol became popular ODE benchmarks precisely because they could make a lot of stiff solvers fail. If this is such an equation, then it would nice to have a reference version as a benchmark and be able to track how well we do, and find out what we can do to improve it, or at least what kinds of options are required to handle it well.

[devmotion]

Is this a numerical issue?

Maybe to some part, but I already could improve the results by setting element-wise absolute tolerances (the variables are on quite different scales). However, only with the additional effort of restricting the solution to positive values, as discussed in SciML/DiffEqCallbacks.jl#12, I could obtain quite good results.

If this is such an equation, then it would nice to have a reference version as a benchmark and be able to track how well we do, and find out what we can do to improve it, or at least what kinds of options are required to handle it well.

I will reevaluate how well it performs since we changed the algorithm meanwhile, so the plots shown in SciML/DiffEqCallbacks.jl#12 have changed probably.

[devmotion]

I just did some plots (without any benchmarking regarding time or memory) of the critical variable in my system. It is very small compared to the other four variables, so without element-wise absolut tolerances and without any other options I usually do not get a satisfying solution:

BS3

Tsit5

DP5

DP8

RK4

OwrenZen3

OwrenZen4

OwrenZen5

But wow, I'm impressed by the OwrenZen5 result 👍

[ChrisRackauckas]

But wow, I'm impressed by the OwrenZen5 result

Yes, it controls the interpolation error the best. We should probably update the docs and suggest it as the Owren methods as the defaults for delay equations. Also try RK4 with its residual error estimator.

[devmotion]

And then I repeated these calculation with element-wise absolute tolerances. This improves all results (well, except of OwrenZen5) but still most of them still show some tiny oscillations.

BS3

Tsit5

DP5

DP8

RK4

OwrenZen3

OwrenZen4

OwrenZen5

At first glance, all OwrenZen methods and RK4 outperform the other algorithms here.

[ChrisRackauckas]

This is with no callback or isoutofdomain?

[devmotion]

Yes, exactly. With the callback I could get similar results for DP8 as for the OwrenZen methods shown above, but some of the other methods still showed these tiny fluctuations.

[ChrisRackauckas]

Part of it with Tsit5 could be numerical issues with the interpolant. It's not using Horner's rule yet.

SciML/OrdinaryDiffEq.jl#130

[devmotion]

I added RK4 to the plots above. With element-wise absolute tolerances it gives also a satisfying result.

[ChrisRackauckas]

Can we make this into a notebook?

[devmotion]

I'll transfer it to a notebook 👍 Which solver would you suggest for computing reference solutions?

In the references you posted I also discovered some other equations for benchmarking, so probably we can add some other notebooks as well.

[ChrisRackauckas]

Which solver would you suggest for computing reference solutions?

Well, with lazy interpolants working Vern9 at very low tolerances should be very accurate. But I think for this one the plots are a very good indicator of qualitative behavior in favor of the OwrenZen methods which should be showcased somewhere.

[ChrisRackauckas]

Cool. Future issues on this will be more specific.