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?

  • $\begingroup$ Integrate[f[x, y], x] /. x -> 1? $\endgroup$ – Michael E2 Dec 14 '17 at 15:00

You can use Integrate[ifn[x], x] on an InterpolationFunction 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]

Mathematica graphics

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

Mathematica graphics

|improve this answer|||||
  • $\begingroup$ This is spot on, I hadn't realised that Interpolating functions could be integrated in this way, thanks! $\endgroup$ – Jiles Dec 15 '17 at 10:02

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]},
    {t, 0, 2},
    {tau, 0, 5}

A comparison with your f function:

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

enter image description here

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:

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

enter image description here


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}]

enter image description here

|improve this answer|||||
  • $\begingroup$ Oops, I was about to post my answer a couple of hours ago when I got called away. I didn't see your answer until I posted. Oh well... $\endgroup$ – Michael E2 Dec 14 '17 at 18:02
  • $\begingroup$ Thanks for the answer, I hadn't appreciated you could use NDSolve in such a way. As it turns out, your answer inspired me to alter some other code, speeding it up 1000 fold! $\endgroup$ – Jiles Dec 15 '17 at 10:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.