3
$\begingroup$

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["DifferentialEquations`InterpolatingFunctionAnatomy`"];

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["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Integrate[
 ifreap[t], {t, InterpolatingFunctionDomain[ifreap][[1, 1]], 
  InterpolatingFunctionDomain[ifreap][[1, 2]]}]
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – flinty
    Commented Jul 15, 2020 at 16:00
  • 1
    $\begingroup$ Note that Last[intFunc[t]] equals t, so those integrals are justing integrating t. Maybe Indexed[intFunc[t], -1] is what you're after. $\endgroup$
    – Michael E2
    Commented Jul 15, 2020 at 16:05
  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    Commented Jul 15, 2020 at 16:23
  • $\begingroup$ @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. $\endgroup$ Commented Jul 15, 2020 at 17:56
  • $\begingroup$ I consider this question answered. $\endgroup$ Commented Jul 15, 2020 at 17:57

1 Answer 1

3
$\begingroup$

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  *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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