Fixed edge case infinite loop on stiff-ish problems (+very bad luck) #86
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Fixes #85.
This turns out to be a "fun" mix of adaptive time stepping and numerical stability.
The problem was that the solver would enter an infinite loop right at the end of integration. The scenario is that a step was proposed which was (a) rejection-worthy and (b) within the 1e-6 clip-to-end behaviour used for numerical stability of dense interpolation, and (c) the proposed-new-step was also rejection-worthy and within the clip-to-end behaviour. And one that only triggers on at least mildly stiff problems. And of course, the proposed-new-step got its step size increased by the clipping and therefore remained rejection-worthy and hence the infinite loop. Quite the edge case!
The fix is to only trigger clipping after an accepted step. This is probably nearly as good as it gets: if a step is rejected, and the proposed-new-endpoint is clippable, then certainly the previous endpoint was clippable as well. (As it is strictly larger.) But we already know that step got rejected, so we should ignore clippability of the proposed-new-endpoint; else we'd just be repeating our failed step.
Unfortunately, this does mean that any subsequent steps may be below the threshold for numerical stability of dense or
ts
output. Mostly this is occurs when usingSaveAt(ts=<something containing t1>)
.So to handle that case, in the case that we're (a) rejecting a step and (b) the new step is clippable, we override the adaptive step size controller and force-propose a step that's exactly half the length of the distance to the end. This is both (a) a step reduction, so possibly an acceptable step, and (b) maximises floating-point stability of both subsequent half-steps.
Finally, this PR also reduces the tolerance for float64, from 1e-6 down to 1e-10.
Phew. Hopefully that works!