Integrate only one variable of a 2D interpolating function

Question

Lets say I have generated a 2D interpolating function using the code:

da = Flatten[Table[{t, tau, N@Sin[2 (t + 3 tau)] Exp[-2 t - tau]}, {t, 0, 2, 2/100}, {tau, 0, 5, 5/100}], 1]
f = Interpolation@da

and I want to integrate over one of the variables (I realise that in this case the integral can be done analytically, but I want to integrate an interpolating function from NDSolve which is too complicated to post here...). I can do this using something like:

tend = 2;
tauend = 5;
dtau = tauend/100; 
taurange = Range[0, tauend,dtau];
int = Table[NIntegrate[f[t,tau],{t,0,tend}],{tau,taurange}];
intfunc = Interpolation@int;

but this seems somewhat inelegant, and also is pretty slow. Is there a better way to approach such a problem?

Answer 1

You can use Integrate[ifn[x], x] on an InterpolatingFunction ifn with domain {x1, x2} to compute the equivalent of the indefinite integral Integrate[ifn[t], {t, x1, x}]:

da = Flatten[
   Table[{t, tau, N@Sin[2 (t + 3 tau)] Exp[-2 t - tau]}, {t, 0, 2, 
     2/100}, {tau, 0, 5, 5/100}], 1];
f = Interpolation@da;
{{x1, x2}, {y1, y2}} = f["Domain"];

intx = Integrate[f[x, y], x] /. x -> x2;

nintx[y_?NumericQ] := Module[{x}, NIntegrate[f[x, y], {x, x1, x2}]];

Plot[nintx[y], {y, y1, y2}, PlotRange -> All]

Plot[intx, {y, y1, y2}, PlotRange -> All]

Answer 2

I think the best way to do this is to have NDSolve output the desired integral directly. Let's take your example and convert it into a differential equation:

g[t_, tau_] := Sin[2 (t+3 tau)] Exp[-2 t-tau]

base = NDSolveValue[
    {D[u[t, tau], t] == D[g[t, tau], t], u[0, tau] == g[0, tau]},
    u,
    {t, 0, 2},
    {tau, 0, 5}
];

A comparison with your f function:

Row[{
    Plot3D[base[x, y], {x, 0, 2}, {y, 0, 5}, ImageSize->300],
    Plot3D[f[x, y], {x, 0, 2}, {y, 0, 5}, ImageSize->300]
}]

Augmenting the NDSolve with an integral:

{base, integral} = NDSolveValue[
    {
    D[v[t, tau], t] == u[t, tau], v[0, tau] == 0,
    D[u[t, tau], t] == D[g[t, tau], t], u[0,tau] == g[0, tau]
    },
    {u, v[2, tau]},
    {t, 0, 2},
    {tau, 0, 5}
]; //AbsoluteTiming

{0.002504, Null}

Quite a bit faster than your approach. A comparison with your intfunc function:

Row[{
    Plot[integral, {tau, 0, 5}, ImageSize->300],
    Plot[intfunc[tau], {tau, 1, 101}, ImageSize->300]
}]

Addendum

As @MichaelE2 suggest, you could also symbolically integrate the interpolating function. Instead of Integrate, I use Derivative:

h = Derivative[-1, 0][f];
Plot[h[2, tau], {tau, 0, 5}]