[WIP] adding output at specified points

[berceanu]

I just added the initstep, minstep and maxstep keywords specified in the API.
Output at specified time-points is still wip, as for a tspan of [0.:10.:50.] I currently get output at

tout = [0.0 6.0e-7 2.351e-5 0.0005 0.009 10.013 19.981 30.016 39.983 49.997 50.0]

[ivarne]

Why do you have abs in eps(abs(...))? According to the documentation eps(x) is the difference between x and the next larger number.

[berceanu]

@tlycken the "prob" with Base.isapprox() is I didn't know it had configurable tolerances. just love it when i reinvent the wheel :P sorry about that!

[berceanu]

@ivarne ty, fixed the signature
@tlycken using norm instead is trivial to implement, but could you detail on basing the step estimation on stability properties?

[coveralls]

Coverage increased (+0.04%) when pulling 8684b86 on berceanu:apistep into 913a3a5 on JuliaLang:master.

[berceanu]

any ideas on fixing the points = :specified behaviour? i'm refering to the extra points which should not be included in tout.

[ivarne]

I think points=:specified should do something smart about interpolation. For such an exotic behaviour, accuracy is more important than efficiency, so we could spend some extra time calculating a better guess for the values at the specified points. How about picking the 3 closest values and using a 2. degree polynomial for estimation?

[ivarne]

As we also have the derivatives, we could probably get a higher polynomial from 3 points.

[berceanu]

the polynomial interpolation that you refer to is normally referred to as dense output. i thought this was a separate issue of the api. in fact i just did a pr on api.md about it. the thing is, for 5th order rk, it can be done without further function evaluations, but not for higher orders. perhaps we neee a separate implementations depending on the method order?

[acroy]

@ivarne: this is what dense output is actually about, doing interpolation
in a smart way and compatible with with the order/accuracy of the RK method
;-) As an alternative we could just say that :specified adjusts the time
steps such that the given times are attained. That would be simple, but
most likely inefficient. In addition we could have points=:dense for dense
output.

@TLYKEN, @berceanu: I think the initial step is typically estimated by some
lower order RK method (ode23 does already something like this and dop853
has this as well). Regarding tspan[end]=Inf: obviously this is not
problematic for maxstep. And minstep is anyways somewhat ill defined, but
seems to be used by other implementations; might be better to have a
maximum number of step-size refinements?

On Wednesday, June 25, 2014, Ivar Nesje [email protected] wrote:

As we also have the derivatives, we could probably get a higher polynomial
from 3 points.

—
Reply to this email directly or view it on GitHub
#43 (comment).

[berceanu]

it looks like 3rd order hermite interpolation using the derivatives we already calculated should Just Work for any method order p

[tomasaschan]

Instead of interpolating for points=:specified, can't we just make sure that the solver steps exactly to the specified points? I'm thinking something along the lines of

if t + h > next_specified_point
    h = next_specified_point - t
end

and then only add the selected points to the output.

Re. step sizes: my main concern in the case of tspan[end]=Inf is if we base any step size guess on tspan[end]-tspan[1] we're going to end up with (Inf - tspan[1])/x == Inf. For maxstep this might not break anything, other than making the max step useless, but for either of the other two things will definitely break. Estimating based on a lower-order method might be a good way to do it; otherwise, it should in principle be possible to decide on some relative change in x that we allow, evaluate the system function once to see how large the changes are, and take a step that gives approximately the correct magnitude of the change in x. Perhaps something like this:

dx0 = F(tspan[1], x0)
h = 0.01 * norm(x) / norm(dx0)

With this approach, the initial step with length h will change x only by something on the order of 1%.

[acroy]

@tlycken :

can't we just make sure that the solver steps exactly to the specified points?

This is what I meant by "adjusts the time steps such that the given times are attained." :-) One might argue that the same effect can be achieved by calling odeXX in a loop and using only yout[end].

evaluate the system function once to see how large the changes are ...

That would basically correspond to the lowest order method (the estimate should also depend on the tolerances). I guess dopri5 et al are using 2nd order to get a better estimate, but in any case there is certainly some heuristics involved. minstep is IMO more complicated to do "right".

[tomasaschan]

@acroy: Great minds think alike! (But I guess they're only half great when it takes so long to successfully communicate that...)

That would basically correspond to the lowest order method

Ah, that makes sense - I don't know exactly what I thought "using a lower order method to estimate step size" meant in practice, but it's probably the way to go then. It seems to me that creating good heuristics for the step sizes is complicated enough to warrant its own discussion; I say let's fork that issue out into its own, and focus on the rest of the changes here. We can always improve that part of the code later.

@berceanu You mention problems with points=:specified, but I don't see that option referred to in the code anywhere. Do you have a buggy implementation not yet committed, or do you need any help in figuring out how to attack the problem?

[berceanu]

I propose to move the interpolation/dense output discussion to #40 and restrict this to approaches such as that suggested by @tlycken that make the solver pass exactly through the suggested set of time steps.

[coveralls]

Coverage decreased (-0.01%) when pulling 7548518 on berceanu:apistep into 913a3a5 on JuliaLang:master.

[berceanu]

I just implemented an algorithm based on explicit Euler stepping for initial step heuristics.
Still a bit stuck on the points=:specified issue.

[berceanu]

@tlycken the points=:specified is implicit in the code, because output is only appended if t is in tspan.

if points == :all || approxin(t, tspan; atol=.02)

[tomasaschan]

@berceanu Hm... that's not the behavior I'd expect from points=:specified. If I do tout, yout = ode45(F, x0, tspan; points=:specified) then I'd expect tout .== tspan, exactly, i.e. I'd expect the solver to give me solution points exactly at the points in tspan, and only there. Part of the reason this is important, is that it makes the exact size of tout and yout predictable: they're the same length, exactly, as tspan. This is important in applications where you want the solution on a specific grid.

With your approach, quite surprising things could happen. For example, if I pass tspan=[0,0.5,1]; points=:specified and the solver needs to take quite small steps close to t=0.5, it's quite possible that more than one point matches approxin(t,tspan;atol=.02); accurate control of the choice of solution points is lost.

Instead, if points==:specified, I think we need a check at each time step, to make sure we don't step past a specified point, similar to what we do to ensure we don't step past tspan[end].

[acroy]

I think it would be great to have hinit in the package. Regarding points=:specified I would say discuss/handle this separately. In #47 I used an interpolation scheme to get the results at the specified points and I think this makes most sense if there is such an interpolation (for RK methods this is typically known).

[mauro3]

This has been superseded, as far as I can tell, closing. (Please let me/us know if that was wrong.)