DIfferentialEquations with the driving function defined on a fixed grid

Certainly you will have to do some kind of interpolation (linear, cubic, etc.), e.g. via Interpolations.jl or Dierckx.jl. However, beware that such interpolations have discontinuities in some derivatives at the “knots” of the interpolant (typically knots = the data points or a subset thereof). This will greatly degrade the accuracy/efficiency of the solver algorithms in DifferentialEquations.jl, since they assume smoothness. A workaround is to pass a tstops parameter listing all of the knots, which will allow the solver to integrate over the smooth segments.

If you have a lot of knots (finely spaced data), however, breaking integration at every knot might cause the solver to use more steps than it would for a similar ODE with smooth data.

An alternative is to find an interpolation that is smooth everywhere, but this is hard to do in general. Simply fitting to a high-degree polynomial is likely to suffer from a Runge phenomenon in the common case of equally spaced data points. If your data represents a smooth function and if you can choose the grid of data points, then by choosing a set of Chebyshev points you can use a smooth Chebyshev interpolant; ApproxFun.jl may help with this. Alternatively, if your data are noisy, then you might consider some kind of smooth fit.

(In the special case where you are just computing a definite integral, i.e. a separable ODE, then you can use QuadGK.jl to design a specialized quadrature rule that takes your nonsmooth interpolant into account a priori.)

5 Likes