Hi, I need an indefinite integral interpolation, to be constructed once and evaluated many times. Basically, I want to solve the outer ODE \dot f(t) = A(t) b(t), where the inner integral A(t)=\int_0^t a(t')dt' is to be solved numerically. Because Integrals.jl only gives one value for a fixed t=T, I thought of using ODE solvers:
using OrdinaryDiffEqs
integral_rhs(_u,_p,t)=a(t)
p = ODEProblem(integral_rhs, 0., (0., T))
A = solve(p, Tsit5())
Then A(t) can be evaluated efficiently by interpolation on the interval (0,T).
This works fine, but only for smooth enough a. In my use case with discontinuous a, julia hangs, and I think it’s the ODE solver algorithm that doesn’t like that.
Can I just use another solver? Or is there a better way to do this that I did not think about?
If you know where the discontinuity is, you should tell the solver where it is; otherwise the problem is much harder (and you should probably specify a larger error tolerance). This is true for any numerical method — integration methods are usually designed for smooth functions.
In general, you almost never want to evaluate definite integrals by solving ODEs — specialized integration routines are much more efficient (and simpler to use), because they exploit the fact that you have “random access” to your integrand values at any point in the domain (rather than having to build it up from “left to right” as in a general ODE). To be able to evaluate the integral at arbitrary points in a given interval, see How do I do a fast cumulative integration - #6 by stevengj
Ok, that’s helpful. I should be able to pass the locations of the discontinuities, as these are known. I was already thinking to use the d_discontinuities keyword argument for the ODE solver for that.
But per your point 2., it would still be more efficient to construct a (spline) interpolation first using one of ApproxFun.jl, BSplineKit.jl or some such, which would need to accept the list of discontinuities, and then integrate that via functionality built-in to those packages. Sounds doable.
Actually, I care mostly about the evaluation performance, it’s not a big deal if the initial integration is not fast, as this is a one-time upfront cost.
Splines have fixed polynomial degree and are low-order accurate, vs. Chebyshev polynomials (e.g. in ApproxFun or FastChebInterp.jl) which are exponentially accurate (in the degree) for smooth function. (Or piecewise Chebyshev if you have a discontinuity.)
Unless you have zillions of discontinuities and a high error tolerance, so that on each small segment you are happy with a low-degree polynomial (e.g. a cubic spline)? Where does your integrand come from? What is your error tolerance?
The computational cost of evaluating the approximant is closely related to the number of points you use for the integration. So a faster-converging interpolation will also sometimes be faster to evaluate. (Though not always. e.g. a low-order spline might have O(\log n) cost in the number of knots, or even O(1) if it is a regular grid of knots, whereas evaluating a polynomial interpolant with n points has O(n) cost.) It depends on what your function looks like and your error tolerance.
The integrand is a well-behaved piecewise smooth 1d function that models a time-dependent reaction rate. Not necessarily a low-order polynomial on each piece but not far off. There is a handful of discontinuities at most. It should be numerically an easy problem. Error tolerance maybe 1e-5 rel error to have some margin for later ODE? Not sure. Would you recommend Chebychev in this situation?
Browsing the docs for FastChebInterp.jl and ApproxFun.jl I have not found a convenient way to construct an interpolation function with discontinuities from piecewise defined input though. This may be by design, because then you can’t define a global Derivative operator?
So right now it looks like I would have to handle the discontinuities in the integrand a manually by separately interpolating each piece of it and then wrapping the pieces in an output interpolating function, and manually handling the piecewise integration of that. More tedious than I had hoped for.
I’ll start with passing discontinuities to the ODE approach I was using. If slowness continues to be an issue I may need to do the work.
after some digging i think the following would give me a foot in the door with ApproxFun.jl
using ApproxFun
function approx_piecewise(f::Function, points::Vector{T}) where T<:Real
segs = (Segment(a, b) for (a, b) in zip(points, points[2:end]))
Fun(f, union(segs...))
end
here, points is the list of points of discontinuity, flanked by the left and right boundaries of the support of the function f to be interpolated