How can I numerically pre-compute an indefinite integral with a parameter?

Suppose I have a function $f(t)$, and I want to compute its indefinite integral $$F(t)=\int_0^tf(\tau)\mathrm d\tau.$$ Moreover, suppose that, for any of a number of reasons, I require this integral to be done numerically, I need it to be done in a single swoop for all $t$ in a range $[0,T]$ of interest, and I need to be able to evaluate $F(t)$ for any $t$ in that range.

One standard trick (which is perhaps not as well known as it should be) is that one of the better ways to do this is to reinterpret the indefinite integral as an initial-value problem, $$F'(t)=f(t)\quad\text{under }F(0)=0$$ and then simply give it to the differential-equations solver NDSolve directly, e.g.

NDSolve[{F'[t]==f[t], F==0}, F, {t, 0, T}]


and this will return $F$ as an InterpolatingFunction object which can be evaluated anywhere in $[0,T]$.

I would like to extend this to the case where the indefinite integral comes with a parameter, say, of the form $$F(t,t')=\int_{0}^tf(\tau,t')\mathrm d\tau.$$ Is there an equivalently clean way to obtain $F$ as an indefinite integral?

One such way is to invoke a similar solver, ParametricNDSolve, in the form

ParametricNDSolve[{F'[t]==f[t,tt], F==0}, F, {t, 0, T}, {tt}]


which will return a solution of the form {f-> ParametricFunction}, where the ParametricFunction is evaluated as ParametricFunction[tt] (for a numeric tt) to give an InterpolatingFunction object, so e.g. ParametricFunction will return a number.

I find this solution unsatisfactory because it is very slow when one needs to repeat the process over and over for different $t'$s. While the ParametricFunction object it returns is nice (including, for example, the ability to compute derivatives like $\frac{\partial F}{\partial t'}$), it is essentially based on a single ray on the $(t,t')$ plane, and there's no way for the different rays to talk to each other, with the downside that Mathematica ends up calculating many things which are very similar in a highly sub-optimal way.

I would like, then, clean solutions where I can run the whole calculation once, and then simply query for the values of $F(t,t')$ without incurring in additional computation.

• A representative example of an integrand and parameter range would be helpful. – Daniel Lichtblau Feb 5 '16 at 18:36
• @DanielLichtblau I would rather keep this a bit open to make it more generally useful, but if something specific is required the cosine example in my answer will do fine, and even better if it can avoid these domain issues. – Emilio Pisanty Feb 5 '16 at 18:51
• I could be misunderstanding, but this does not look like a PDE scenario. Rather it seems to be a parametrized ODE that gets treated numerically in a parameter range by independent single ODEs. As such I would (naively) expect better reliability numerically than from a PDE. (If I am correct then the method you show might be generally quite good.) – Daniel Lichtblau Feb 5 '16 at 19:28
• @DanielLichtblau Indeed, it's just an ODE at heart. However, I'm struggling to 'glue' together all the different single ODEs to get a fast calculation of the whole thing. In particular, if $f(t,t')$ doesn't change much for any $t$ for all $t'$ between $t_1'$ and $t_2'$ and I've already run the ODE for $F(t,t_1')$ and $F(t,t_2')$, I would like a scheme that doesn't need to re-run the entire ODE every time I ask for $F(t,t')$ for some $t'$ in between the two. – Emilio Pisanty Feb 5 '16 at 19:40
• The PDE approach does sometimes work but it feels really clunky to me. (And, indeed, since it's not really a PDE, it does cause NDSolve to stop and say "hey!" once things get slightly complicated.) Thus I'm very interested in other approaches. – Emilio Pisanty Feb 5 '16 at 19:41

Here is one way (but I remain interested in other ways to solve the problem).

The indefinite integral can still be reinterpreted as a differential equation, of the form $$\frac{\partial F}{\partial t}(t,t')=f(t,t') \quad \text{under }F(0,t')=0$$ and, though this isn't obvious, it can still be handled as a differential equation by NDSolve - except this time it needs to be treated as a PDE.

To set a more concrete example, consider the integrand $f(t,t')=\cos(t)-\cos(t')$. In this case, simply giving the differential equation to NDSolve and using the PDE syntax (with two independent variables)

NDSolve[
{D[F[t, tt], t] == Cos[t] - Cos[tt], F[0, tt] == 0}
, F
, {t, 0, 50}, {tt, 0, 50}
]


evaluates perfectly fine, returning a two-dimensional InterpolatingFunction which reproduces the correct solution: In terms of what's happening, Mathematica seems to interpret the equation as a convection-diffusion equation with a vanished diffusion term, which can occasionally cause some instabilities and some jumpy results; these can sometimes be eliminated by adding some artificial diffusion (a small multiple of $\nabla^2F$ on the right-hand side) which can help stabilize the solver without appreciably changing the results.

However, this stabilization does not always work, particularly with complicated domains, so there's definitely a good bit of room for improvement.