WIP: Fix different solutions #17
Conversation
src/integrator_interface.jl
Outdated
| # separate k to avoid influence of calculations on interpolation | ||
| # constant caches that contain fsalfirst must be updated for interpolation | ||
| # has to be changed all such caches - move to OrdinaryDiffEq? | ||
| if typeof(integrator.cache) <: Union{BS3ConstantCache} |
There was a problem hiding this comment.
Just forgot to remove Union after experimenting with Tsit5 😄
There was a problem hiding this comment.
To update k[1]. For some caches like BS3ConstantCache it contains fsalfirst - however, at this point of the algorithm its an old, outdated value of fsalfirst, and hence also any interpolation that uses k[1] is wrong. There are no problems with BS3Cache since it only points to fsalfirst and hence the value of k[1] always equals fsalfirst.
Of course, one has to add all affected caches (there are many more besides BS3ConstantCache - all not-inplace caches that contain fsalfirst and use it for interpolation) to this Union. I just wanted to wait til I know whether you would like to apply this update of k here or rather in OrdinaryDiffEq.jl when fsalfirst is updated.
There was a problem hiding this comment.
won't this make k[1] == k[2] in the first round, making the interpolation/extrapolation invalid?
There was a problem hiding this comment.
Yes, it will make k[1] == k[2] in the first round. However, I'm still not completely sure about which values one would like to use for the interpolation - I mean, I know what we want to achieve but the problem is the correct assignment of u, uprev, t, and tprev for integrator.integrator. With the current (and the old) setup we have u == uprev (by the ODE algorithm) and t == tprev (by our assignments) in the first round. In the next rounds u is always updated to the newest estimate, uprev stays fixed, t is fixed to integrator.t + integrator.dt, i.e. the time point after this step, and tprev is fixed to integrator.t, i.e. the current time point.
There was a problem hiding this comment.
it should never have k[1] == k[2]. That would make it a constant extrapolation instead of the Hermite. In the first round it should stay in the previous setup, unchanged. After the first iteration, tprev = integrator.t, t = integrator.t + integrator.dt, uprev should be u(tprev), and u should be u(t)
There was a problem hiding this comment.
That's exactly one of the reasons why I was confused about this part. So the problem in our case is rather that for in-place algorithms with k[1] = integrator.fsalfirst and k[2] = integrator.fsallast (like BS3Cache) we always automatically have k[1] == k[2] at this point.... Maybe it's best to always separate k by integrator.integrator.k = recursivecopy(integrator.k) right after initialization and then use recursivecopy! in all other places? But I guess I have to try it to see whether it actually works...
src/integrator_interface.jl
Outdated
| # update ODE integrator | ||
| integrator.integrator.uprev = integrator.uprev | ||
| integrator.integrator.tprev = integrator.tprev | ||
| integrator.integrator.fsalfirst = integrator.fsalfirst | ||
| integrator.integrator.tprev = integrator.t # to extrapolate from current interval |
There was a problem hiding this comment.
This can't update yet: it doesn't have valid values to extrapolate with.
There was a problem hiding this comment.
Ah OK, then it is also incorrect on the master branch, since it just shortens lines 41 and 45 of the current implementation. The value at time integrator.t is integrator.uprev = integrator.u at this point - this is also not useful for the extrapolation... It should rather be integrator.uprev = u(integrator.tprev) (not updated tprev), shouldn't it?
There was a problem hiding this comment.
Oh wait, was integrator.integrator.tprev never updated in the previous round? I think I'm going to have to run through it all again and print values out because it's definitely tricky.
src/integrator_interface.jl
Outdated
| # integrator.EEst = integrator.integrator.EEst | ||
| #end | ||
| # update ODE integrator | ||
| integrator.integrator.u = integrator.u |
There was a problem hiding this comment.
To update the ODE integrator. Otherwise integrator.integrator.u still contains u(integrator.t), when peforming a step without iterations, or u(integrator.t+integrator.dt) after the second to last iteration, when performing fixed-point iterations. For any interpolations that occur after this step it should rather be the current, updated value u(integrator.t+integrator.dt).
|
Since both the current implementation and the first draft of this PR are not completely correct, I try to summarize both what we want to achieve and my suggestions/ideas, even though it is quite lengthy. I want to provide a rough overview of our desired algorithm using an example, hopefully as an assistance and a reminder to myself. Maybe it can also speed up the process of fixing the discovered problems. I surely made some mistakes, so please correct anything that is wrong. ExampleAssume we already calculated a solution u(t) up to time point t1, and the last step of our computation we moved from time point t0 to t1 such that [t0, t1] is the last calculated interval. Let dt be the proposed next time step, i.e. we would like to advance to time point t1 + dt in the next step. Moreover, let k[a, b](x, y) be the interpolation data of time interval [a, b] with function values x and y at a and b, respectively. Then our algorithm works in the following fashion:
Now, if there are still stopping time points left but there is no stopping time point between t1 and
When we reached the end point, the solution is cleaned by
Considerations
|
|
I tried to implement exactly what I wrote down above. It seems to work quite well, I could improve some of the error bounds in the tests significantly. Partly this might be due to the improved extrapolation and due to the fact that the error estimate now combines both the error estimate of the integrator and the error estimate of the fixed-point iteration, which might lead to a higher rate of unaccepted time steps and thus sometimes smaller time steps. It is even possible to use methods with lazy interpolants (#16) since additional steps are added to the interpolation before values are saved or a new step is calculated. I tested the problem Two other problems I faced:
Moreover, with this changed implementation I get the best error estimates for my toy example so far; before (#15 (comment)) the error estimates were julia> sol.errors
Dict{Symbol,Float64} with 3 entries:
:l∞ => 3.33289e-5
:final => 1.86351e-5
:l2 => 1.40652e-5
julia> sol2.errors
Dict{Symbol,Float64} with 3 entries:
:l∞ => 3.31446e-5
:final => 1.86005e-5
:l2 => 1.401e-5now they are julia> sol.errors
Dict{Symbol,Float64} with 3 entries:
:l∞ => 2.98268e-5
:final => 1.7883e-5
:l2 => 1.30246e-5
julia> sol2.errors
Dict{Symbol,Float64} with 3 entries:
:l∞ => 2.98268e-5
:final => 1.7883e-5
:l2 => 1.30246e-5The tests pass locally but need SciML/OrdinaryDiffEq.jl#80 which is merged but not released yet. |
|
Tests pass on Windows 64bit, Linux, and Mac - but the stricter error bounds lead to failures on Windows 32bit. Is there a workaround without sacrificing the better bounds completely? |
|
I still think that the outline above summarises what we want to do - but does the extrapolation really work correctly? I found some papers regarding our approach, e.g. https://link.springer.com/article/10.1007/BF02072017 and http://www.sciencedirect.com/science/article/pii/S0168927497000226, and because of the obvious formula for First of all, every extra- and intrapolation depends on both the last time point of the So my question is: don't we have to update But maybe I'm also just confused right now... It's really tricky to get it right it seems. |
We just get a I'd point to Hairer II or the papers on LSODE/Sundials as arguments for this as being a good starting scheme for implicit methods. |
|
Ah OK, good to know. I am just worried there might still be some failures in the algorithm - in particular since suddenly the iterations were exploding for some special cases. And I was thinking whether one can stop the fixed point iterations as soon as the error estimate increases between two subsequent iterations. |
We should switch over to Anderson first, and this could probably be something in NLsolve.jl and we then interpret a retcode. |
|
Where exactly are we on this PR? Ready for review again? Has some merge conflicts after the other PR |
|
Finally I fixed the conflicts, so it is ready for review again. Took me awhile since I ran into PainterQubits/Unitful.jl#95 and I figured out that in fact the interpolation was done incorrectly - Moreover, due to this change lazy interpolants do not work anymore, so for now I completely removed the code which was supposed to add additional steps and handle these algorithms. |
|
Again HDF5 build errors on Windows 32bit but it seems these issues are common and known: JuliaPackaging/WinRPM.jl#120, JuliaPackaging/WinRPM.jl#118, JuliaPackaging/WinRPM.jl#107, and JuliaPackaging/WinRPM.jl#117 |
src/integrator_utils.jl
Outdated
| # update time interval of ODE integrator | ||
| integrator.integrator.t = integrator.t | ||
| integrator.integrator.tprev = integrator.tprev | ||
| integrator.integrator.dt = integrator.integrator.t - integrator.integrator.tprev |
There was a problem hiding this comment.
why not integrator.dt? This could be open to floating point issues.
src/solve.jl
Outdated
| @@ -75,7 +75,7 @@ function init(prob::AbstractDDEProblem{uType,tType,lType,isinplace}, alg::algTyp | |||
| ode_prob, | |||
| OrdinaryDiffEq.alg_order(alg))) | |||
| end | |||
| integrator.dt = dt | |||
| integrator.dt = zero(dt) | |||
There was a problem hiding this comment.
won't this override the init dt choice from the user?
There was a problem hiding this comment.
No, this just assures that the ODE integrator always satisfies tprev + dt == t. The DDE integrator is still initialized with the init dt choice from the user (or the calculated init dt).
There was a problem hiding this comment.
Setting dt = t - tprev does not assure that the ODE integrator satisfies tprev + dt == t
There was a problem hiding this comment.
I'm not sure I understand you. Of course, for all other steps (besides the initial step) lines https://github.com/devmotion/DelayDiffEq.jl/blob/8008e8f2d0837c944e7a7093a948a4fb094f0bb5/src/integrator_utils.jl#L239-L242 and https://github.com/devmotion/DelayDiffEq.jl/blob/8008e8f2d0837c944e7a7093a948a4fb094f0bb5/src/integrator_interface.jl#L70-L72 guarantee that t, tprev, and dtof the ODE integrator integrator.integrator belong to the same time interval which is important for correct interpolation. So setting dt to zero initially (only for integrator.integrator, the DDE integrator still uses the supplied or calculated dt!) since tprev = t initially is not needed for any calculations but ensures this consistent relationship of these three fields throughout the whole algorithm.
Of course, because of floating point issues the subtraction is a bad idea, and hence I corrected this in the last commit.
There was a problem hiding this comment.
Does this actually have a practical impact, or just makes things nicer because then tprev + dt == t?
There was a problem hiding this comment.
However, for all other steps it makes a huge difference
Not this line though, just the PR in general?
So this specific line could be removed but I thought it would be preferable to keep dt, t, and tprev synchronized all the time.
Nah, it makes sense. Just was curious about the practical effect. In theory it would do a constant extrapolation if dt=0, I just was wondering if that ever came into play.
There was a problem hiding this comment.
Not this line though, just the PR in general?
Yes just the PR in general.
Nah, it makes sense. Just was curious about the practical effect. In theory it would do a constant extrapolation if dt=0, I just was wondering if that ever came into play.
But it is not implemented yet, is it? As far as I see, all calls of current_interpolant and current_interpolant! (and also of current_extrapolant and current_extrapolant!) will calculate a Θ that is not a number for dt=0 - which will always lead to interpolations/extrapolations that are not a number, I guess.
There was a problem hiding this comment.
But it is not implemented yet, is it? As far as I see, all calls of current_interpolant and current_interpolant! (and also of current_extrapolant and current_extrapolant!) will calculate a Θ that is not a number for dt=0 - which will always lead to interpolations/extrapolations that are not a number, I guess.
But almost all of the interpolations have u0 + dt* ..., so it'll "accidentally" work out as a constant extrapolation. I think just the SSPRK methods don't. But I don't think we're actually hitting this case so it probably doesn't come up. Probably worthwhile to address in another PR if it does, with the SSPRK methods/interpolants as the test case which doesn't have the accidental workaround.
There was a problem hiding this comment.
But almost all of the interpolations have u0 + dt* ..., so it'll "accidentally" work out as a constant extrapolation.
I guess it won't, and it didn't in my tests with BS3 and Tsit5. The problem is that the expression in brackets after u0 + dt* is not a number in these cases; an example is line https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/dense/interpolants.jl#L123.
There was a problem hiding this comment.
oh, 0*NaN = NaN? I never knew that. Okay, then I guess that doesn't happen.
Codecov Report
@@ Coverage Diff @@
## master #17 +/- ##
==========================================
+ Coverage 82.08% 82.31% +0.22%
==========================================
Files 9 9
Lines 268 294 +26
==========================================
+ Hits 220 242 +22
- Misses 48 52 +4
Continue to review full report at Codecov.
|
I finally found the missing piece, and could fix #15 for BS3 algorithms. The fix can easily be extended to all other affected algorithms, however, it might be preferable to add it to
OrdinaryDiffEq.jl.Besides the already discussed separation of
uandkinperform_steps!it is necessary to update the value offsalfirstink(if existent): for in-place functions this happens automatically via pointers, whenfsalfirstis updated in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/integrator_utils.jl#L280; for not in-place functionsfsalfirstis updated in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/integrator_utils.jl#L282, but this change is propagated tokonly at the end of the next execution ofperform_step!in https://github.com/devmotion/OrdinaryDiffEq.jl/blob/master/src/integrators/low_order_rk_integrators.jl#L43. Hence the interpolation, which depends onk, differs. This should apply to all not in-place algorithms which savefsalfirstinkand use its entry inkfor interpolation; I could provoke the problem with Tsit5 and apply the same fix. Only the fix for BS3 is included here so far, since I was not sure whether a fix in OrdinaryDiffEq.jl would be more appropriate.The tests pass locally but need SciML/OrdinaryDiffEq.jl#80.