Use keyword arguments to set the tolerance

[mweastwood]

The title is fairly self-explanatory.

[stevengj]

Please use the same names abstol and reltol as in Base.quadgk, since the concept is similar here.

[mweastwood]

I renamed the keywords, but I made a judgment call in the oderkf function. It appears to me that the tolerance that is being used is relative, hence my choice of reltol (previously the variable was simply tol).

[acroy]

@mweastwood : Basically I think this is a step in the right direction. Note that your changes have quite some overlap with @tlycken's API suggestions in #7. Maybe the two PR's could be coordinated somehow?

[tomasaschan]

@mweastwood @acroy @stevengj I think abstol and reltol are good here. However, atol and rtol are used in other (less related) places in Base, so if it matters to be consistent this might be a bigger question than just ODE.jl and Base.quadk. I'm fine with either, and I do think quadk and ODE.jl benefit from similar interfaces, so this is probably a good idea.

Regarding #7, I doubt that it's a good idea to make all changes suggested in that thread at once. I started the discussion there as a pull request mainly to have something to base the discussion on, but it's probably wiser to do it all in small steps. Whatever other changes I made in those commits are going to have to take this one into account anyway, so as far as I care we might as well merge this now and take care of other ideas separately.

[acroy]

@tlycken : I also think starting with this PR would be the right thing. But I would say that your commits in #7 already implement many of the points we want to have eventually (duck typing, norm, coherent arguments). IMO it would be a pity not to use it.

[tomasaschan]

@acroy: Sure. But it's easy to rebase that branch onto something like this if necessary - and since some of the things already implemented in #7 are still under discussion (e.g. function arguments) it might take some time before that can be merged. At least some of the changes in that branch are easy enough to implement to fix that it might be better to do them one at a time instead of all at once - baby steps =)

[acroy]

@tlycken : Sorry for being impatient - let's take it one feature at a time. I guess we shouldn't highjack this PR any more, so let's continue the discussion over at #7.

Regarding this PR: To support abstol in oderkf would it be sufficient to add abstol to tau here? And IMO we shouldn't use kwargs... in ode45_dp, ode45_fb and ode45_ck.

[tomasaschan]

@acroy It depends. In quadgk the condition for accepting an error is E < abstol && E < reltol*f where E is the error and f is the quantity we're looking at, while in other places it's E < abstol + reltol * f. I don't know which I like better.

However, whichever condition we want to use, I think this is the only place where we need to implement it.

[acroy]

@tlycken The easiest solution probably is setting tau=abstol + reltol * f, which automatically takes abstol into account when estimating the new step size. This also has the advantage (?, at least side-effect) that abstol=0 or reltol=0 is possible.

[tomasaschan]

@acroy Sounds good to me. Maybe @stevengj could chime in on whether this makes sense, and why quadgk does it differently?

[acroy]

@stevengj Correct me if I am wrong, but testing reltol and abstol separately (as in quadgk) is more strict than using the sum. I think both criteria are equally valid. Hairer's DOPRI5 for example uses "abstol + reltol * f". We should just make sure that we document the meaning of reltol and abstol accordingly.

[stevengj]

I think of abstol and reltol as separate convergence criteria (potentially just two of many criteria), not as values that should be summed. Hence I would prefer to test them separately, with convergence when either criterion is satisfied.

The condition for convergence in quadgk is the opposite of what @tlycken wrote: E <= abstol || E <= reltol*abs(f), not && (the negation of a conjunction is a disjunction). This is equivalent to E <= max(abstol, reltol*abs(f)).

@acroy: max is slightly stricter than summing, but in practice the two should be very similar, since by summing you can at most double the E at which convergence is achieved. But I don't find doubling the requested tolerance to be a very intuitive behavior, which is another reason why I prefer separate tests.

Note both max and + allow either abstol or reltol to be 0 in order to effectively ignore the corresponding convergence criterion. Note also that you want <=, not <, because E==0 should always be converged.

[tomasaschan]

@stevengj @acroy Yeah, sorry about the misinformation from quadgk there. I must have read it too quickly :P To me, max(abstol, reltol*f) seems the way to go (and maybe we should change it in isapprox as well). Documentation is, of course, key.

@mweastwood Would you mind adding a commit to this PR which fixes this?

[acroy]

@stevengj : Thanks for the clarification! As you say, in practice it probably won't matter much, but I agree that treating abstol and reltol as separate criteria is more intuitive.

@mweastwood : Could you also replace the kwargs... in ode45_dp, ode45_fb and ode45_ck by reltol and abstol? Maybe I am overlooking something, but otherwise the user might provide an unsupported keyword which is silently ignored.

[stevengj]

@acroy, unrecognized keywords won't be silently ignored, because they are eventually passed to oderkf which will throw an error.

[acroy]

@stevengj: Sorry, you are right of course!

[mweastwood]

Following the discussion, I added abstol to ode45, but now the behavior between ode23 and ode45 is inconsistent.

[acroy]

@mweastwood : Thanks for the commit! What do you mean by inconsistent behavior?

[mweastwood]

@acroy : maybe I need to review the ode23 code a little bit more, but it seems to me that ode23 and ode45 are using the relative and absolute tolerances in different ways.

For instance ode23 only checks to see if err <= reltol when making a step, and the functionality associated with the absolute tolerance is incorporated into the threshold variable (whose purpose isn't immediately obvious to me). Whereas the discussion above concluded that we want to make sure the error is less than max(f*reltol,abstol).

[acroy]

@mweastwood : Mhmm, I think err<=rtol in ode23 is basically identical to

 norm(e,Inf) <= rtol*max(max(abs(y), abs(ynew)), threshold) = max(rtol*max(abs(y), abs(ynew)), atol) 

(given the definition of err here) Apart from abs(ynew)) this seems to be consistent, but I guess we could also adopt the notation to match ode45? I have to admit though, that the reason for introducing threshold is completely unclear to me. Maybe this was used in the original implementation?

[acroy]

@mweastwood : I think we should actually replace the convergence stuff in ode23 and match it to what you implemented for ode45. The threshold version does not work for reltol==0.

[mweastwood]

@acroy : Ok, I probably won't get around to this until early next week.

[acroy]

Hang on. Why don't we just use our oderkf to realize the Bogacki–Shampine method (aka ode23)? We just need the right Butcher tableau and make the variable pow a function argument (used to estimate the new step size). Then we would only have one function for ode45 and ode23 to worry about. I just tried this locally and it seems to work nicely. But maybe there is a caveat?

[tomasaschan]

@acroy I haven't looked at it at all, but it sounds to me like a good idea; less code to maintain is better code to maintain. However, we should make sure we carefully compare the current implementation versus using oderkf with a different Butcher tableau and pow, to make sure we don't loose performance if/when we make the switch. If it works well, I'm thinking it might be a feasible idea to generalize oderkf enough to be (almost) any RK-method, and then implement all other RK methods as wrappers that send the right tableau and pow. That would be an easy way to, for instance, add an equivalent of MATLAB's ode113.

BUT it's a big thing to do, and as such out of scope for this PR =)

[acroy]

@tlycken : Indeed, doing just the minimal changes in oderkf (introducing an argument pow and some adjustments of the tolerances to match ode23) gives (slightly) worse performance than ode23. With some optimizations (mainly de-vectorization) the performance becomes comparable if not better (as far as I can tell). Since the Butcher tableaus are fixed, maybe one could even use code-generation magic to further optimize this? If you are interested I can post the modified code as a gist ...

Anyways, this should definitely be a new PR. But if we are going in this direction, we might just use @mweastwood's changes as they are (maybe add some comment/warning that reltol==0 in ode23 gives zero step size).

[tomasaschan]

@pao I think this is good to merge.

[acroy]

@tlycken : I guess it has to be rebased & squashed first?

@mweastwood : In view of #16 it is very likely that ode23 is soon replaced. I would suggest to just add a check reltol!=0 to ode23 and then leave it at that.

[pao]

@tlycken @acroy Correct, merge status is not green; the changeset conflicts right now.

[mweastwood]

I believe I managed to rebase+squash my commits (never done this before, so I might have messed up).

Also added a simple check for reltol = 0.

[acroy]

@mweastwood : This looks good to me. I don't know if you have to squash all commits into just one?

[ivarne]

As long as the merge works fine and gives the correct result, there is no reason to be bullish about rebasing and squashing. We should spend our time on writing code to change history, not to prettify the git log. This is anyway a low traffic repo, so the git log will be readable anyway.

[acroy]

@ivarne: thanks for the comment, I just wasn't sure about the policy. Anyways I would me more than happy to see this merged and move on with tackling the other issues.

[pao]

@acroy @ivarne @tlycken One of you should probably request write privileges to this repository from @ViralBShah or @StefanKarpinski, since you are actively working here and seem to know what you're talking about. (Five intentional pings in one comment!)

[ViralBShah]

@acroy added. @ivarne and @tlycken already have access.

[ivarne]

Stefan added us yesterday.