Numerically integrating a list-valued function [duplicate]

Question

I want to NIntegrate a List valued function foo[x] which is only defined for numerical arguments x. Here is a simple toy example:

foo[x_?NumericQ] := {x^2, x^3};
NIntegrate[foo[x], {x, 0, 1}]

It gives the following error

NIntegrate::inum: Integrand foo[x] is not numerical at {x} = {0.00795732}. >>
NIntegrate::inum: Integrand foo[x] is not numerical at {x} = {0.00795732}. >>.

Note that the symbolic definition fooS works fine:

fooS[x_] := {x^2, x^3};
NIntegrate[fooS[x], {x, 0, 1}]
(*{0.333333, 0.25}*)

What causes the error?

It seems to be crucial that the function foo is List valued. When I NIntegrate each component separatly, the error does not occur:

foo1[x_?NumericQ] := foo[x][[1]];
foo2[x_?NumericQ] := foo[x][[2]];
{NIntegrate[foo1[x], {x, 0, 1}], NIntegrate[foo2[x], {x, 0, 1}]}
(*{0.333333, 0.25}*)

and even

NIntegrate[{#1, #2} & @@ {foo1[x], foo2[x]}, {x, 0, 1}]
(*{0.333333, 0.25}*)

Answer 1

This seems like a bug, or at least a "glitch" in NIntegrate. I believe that it expects the evaluated structure of the integrand to match when given symbolic and numeric input. I imagine that it looks at the structure of the output of foo[x] and then sets up the rest of the computation based on that; when it then get a List output from e.g. foo[0] it fails to match the expected pattern and errors.

By providing an integrand that is explicitly a two element list in its evaluated (symbolic) form we eliminate the error:

foo[x_?NumericQ] := {x^2, x^3};

NIntegrate[{Indexed[foo[x], 1], Indexed[foo[x], 2]}, {x, 0, 1}]

{0.333333, 0.25}

Note that unlike Part Indexed will remain unevaluated until its first argument is a List.

The evaluated symbolic form:

{Indexed[foo[x], 1], Indexed[foo[x], 2]}

{Indexed[foo[x], {1}], Indexed[foo[x], {2}]}

Answer 2

The explanation was already delivered by Mr.Wizard, but I would like to add that there is a similar capability to the Indexed approach already built in to NDSolve or NDSolveValue. So we can leverage NDSolve instead of NIntegrate as follows:

foo[x_?NumericQ] := {x^2, x^3};
NDSolveValue[{y'[x] == foo[x], y[0] == {0, 0}}, y, {x, 0, 1}][1]
(* ==> {0.333333, 0.25} *)

Here, the initial condition y[0] == {0, 0} tells Mathematica that the integral y[x] is a vector function.

This approach can also be adapted for version 8 by using NDSolve. It doesn't suffer from the performance hit due to double evaluation of the integrand as in Mr Wizard's answer.

Another side-benefit is that the result of NDSolveValue yields a function that can give you the value of the integral for all intermediate values between x=0 and x=1, just by replacing the argument [1] at the end of the command by a different argument.

Edit

What I describe as a potential advantage of NDSolve could also become a disadvantage, if you're really only interested in the end point of the integration interval and don't want to store the intermediate values of the integral y[x]. As Alexey mentions in the comment, the result of NDSolve will remain unchanged if we replace the specification of the integration domain by the apparently meaningless {1, 1}, because the integration will still be done over the interval {0, 1} if we leave the initial condition y[0] == {0, 0} unchanged. However, the storage requirement for the function y[x] over the domain {1, 1} is potentially much reduced.

Answer 3

I'd argue you should reformulate your function to avoid the issue -- however if a straightforward quadrature scheme will work for you (ie. you don't need adaptive schemes etc ) you can do a direct evaluation:

Simpsons rule for example:

 foo[x_?NumericQ] := ({x^2, x^3});
 np = 99;(*assumed odd for simpsons rule*)
 a = 0;
 b = 1;
 wt = (b - a)/(3 (np - 1)) Join[{1}, 
      Flatten@ConstantArray[{4, 2}, (np - 1)/2 - 1], {4, 1}] // N;
 x = a + (b - a) (Range[0, np - 1] /(np - 1)) // N;
 wt.# & /@ Transpose[ foo /@ x ]

{0.333333, 0.25}

Note that if you do need advanced functionality of NIntegrate you likely want it to work independently on each expression.