3

My goal is to integrate the last dimension of the output of an Interpolation function. I am getting behavior that I don't understand (note to readers, I've pasted all the code without comments at the bottom so you can copy and paste just once):

Data to demonstrate, not real problem:

data = Table[{t, RandomReal[{-.1, .1}, 12]}, {t, 0, 1, .01}];
intFunc = Interpolation[data]

This works as I would expect:

Integrate[intFunc[t], {t, 0, 1}]
Last[Integrate[intFunc[t], {t, 0, 1}]]

A priori, this didn't:

NIntegrate[intFunc[t], {t, 0, 1}]

So, force the integrand to be numerical:

ifunN[t_?NumericQ] := Last[intFunc[t]]
NIntegrate[ifunN[t], {t, 0, 1}]

OK, I get an extra warning, but the result is fine (puzzled why Integrate is doing a different job than NIntegrate though).

Here is the behavior that really puzzles me:

Integrate[Last[intFunc[t]], {t, 0, 1}]  (*returns 1/2 ???*)
NIntegrate[Last[intFunc[t]], {t, 0, 1}]  (*returns 0.5 which is consistent, but why?*)

Let's experiment a bit to see if we can figure out where the weird result is coming from:

{val, {reap}} = 
  Reap[NIntegrate[Last[intFunc[t]], {t, 0, 1}, 
    EvaluationMonitor :> Sow[{t, Last[intFunc[t]]}] ]];

ifreap = Interpolation[reap];
Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];

This yields an approximation to the "expected behavior" (except for a factor of 2)

Integrate[
 ifreap[t], {t, InterpolatingFunctionDomain[ifreap][[1, 1]], InterpolatingFunctionDomain[ifreap][[1, 2]]}]

So, the experiment doesn't help me figure out the the behavior.

All Code if you want to copy-and-paste once:

data = Table[{t, RandomReal[{-.1, .1}, 12]}, {t, 0, 1, .01}];
intFunc = Interpolation[data]

Integrate[intFunc[t], {t, 0, 1}]
Last[Integrate[intFunc[t], {t, 0, 1}]]


NIntegrate[intFunc[t], {t, 0, 1}]


ifunN[t_?NumericQ] := Last[intFunc[t]]
NIntegrate[ifunN[t], {t, 0, 1}]


Integrate[
 Last[intFunc[t]], {t, 0, 1}]  (returns 1/2 ???)
NIntegrate[
 Last[intFunc[t]], {t, 0, 
  1}]  (returns 0.5 which is consistent,but why?)


{val, {reap}} = 
  Reap[NIntegrate[Last[intFunc[t]], {t, 0, 1}, 
    EvaluationMonitor :> Sow[{t, Last[intFunc[t]]}]]];


ifreap = Interpolation[reap];
Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];


Integrate[
 ifreap[t], {t, InterpolatingFunctionDomain[ifreap][[1, 1]], 
  InterpolatingFunctionDomain[ifreap][[1, 2]]}]
Craig Carter
  • 4,416
  • 16
  • 28
  • If I had to guess, Integrate is using the exact Hermite polynomial and manipulating it symbolically, while NIntegrate is evaluating it numerically. I don't know why you think 0.5 vs 1/2 is weird - NIntegrate does numerical integration, so of course it deals with machine numbers. – flinty Jul 15 '20 at 16:00
  • 1
    Note that Last[intFunc[t]] equals t, so those integrals are justing integrating t. Maybe Indexed[intFunc[t], -1] is what you're after. – Michael E2 Jul 15 '20 at 16:05
  • Note also that Integrate[intFunc[t], t] returns an antiderivative in the form of an interpolating function equal to $\int_a^t f(\tau),d\tau$, where $f$ is intFunc and $a$ is the beginning of the domain, namely 0 in this case. – Michael E2 Jul 15 '20 at 16:23
  • @MichaelE2, Yes. You are right. Thanks. NIntegrate[Indexed[intFunc[t], -1], {t, 0, 1}] gives the correct result. Would never have thought of using Indexed. – Craig Carter Jul 15 '20 at 17:56
  • I consider this question answered. – Craig Carter Jul 15 '20 at 17:57
  • Related/duplicate: https://mathematica.stackexchange.com/questions/126342/nintegrate-of-a-vector-valued-interpolatingfunction-gives-not-numerical – Michael E2 Jul 15 '20 at 18:13

1 Answers1

3

Note that Last[intFunc[t]] equals t, so the two confusing integrals are just integrating t from 0 to 1, which is why one gets 1/2. What is needed is Indexed[intFunc[t], -1], which extracts a part only when its first argument is a vector.

NIntegrate[Indexed[intFunc[t], -1], {t, 0, 1}]

(*  0.00378552  *)
Michael E2
  • 235,386
  • 17
  • 334
  • 747