Fast quadrature of solutions via interpolation

[ChrisRackauckas]

We have the interpolation, which means we can easily integrate the interpolation on each interval. Analogous to the derivatives,

sol(t,Val{-1})

can give the integral of the interpolation function. Since the interpolation is continuous, by the Fundamental Theorem of Calculus we can get definite integrals like:

sol(b,Val{-1}) - sol(a,Val{-1})

or indefinite integrals:

sol(t,Val{-1}) - sol(t0,Val{-1})

This means that people who want to perform quadrature on the resulting solution could get the results in O(1) using just two interpolation calls, if we define the coefficients for the integrated interpolation which is easy to do.

[ChrisRackauckas]

To do this, you'd look at the interpolation functions, like:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/dense/interpolants.jl#L151
https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/dense/interpolants.jl#L183

The first is the Val{0} and the second is Val{1}, i.e. the derivative. It's pretty clear how to reverse that via standard calculus.

[ranocha]

Integration is a bit more involved than differentiation because you have to choose an additional constant. Thus, the primitive function in the n-th time step depends on all previous time steps. How should this be handled without loosing too much performance?

[ChrisRackauckas]

Oh that's a very good point. Basically if Val{x} for x < 0, for every previous step it should sum up the Val{x} interpolant on the right endpoint to get the constant. It should then add that constant to the last one. This should still be much cheaper than doing full quadrature though.

[ranocha]

That's probably correct. However, in order to compute definite integrals, a second method might be added in order to avoid the additional overhead of summing up the integrals up to the starting point, since sol(t,Val{-1}) is O(t) and not O(1).

[ChrisRackauckas]

That's a very good point. Or there can be a keyword argument for the starting point which starts at sol.t[1], and modify the summation idea to start at the first interval that this falls in. Care would have to be taken for the first value since that's the [theta,1] integral instead of the [0,theta].