# Numerically integrating a list-valued function [duplicate]

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}*)

• Thanks for the Accept, but I suggest that you wait at least 24 hours to let everyone around the world have a chance to answer. You may like a different answer more if you don't discourage someone from posting by making the issue appear concluded. Aug 18, 2014 at 18:57
• Related: (11772), (14645), (21662), and even: (1731) Aug 18, 2014 at 19:37
• Everyone who voted for my answer or Jens' answer should vote for Michael's earlier answer as it already covered both aspects of the issue. Aug 18, 2014 at 19:57

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

• note this is evaluating each of the list expressions twice for each x .. so it illustrates the issue but isn't likely a good approach for performance reasons. Aug 18, 2014 at 19:01
• that issue is somewhat avoided if you do : foo[x_?NumericQ] :=foo[x]= {x^2, x^3} ( though not entirely in case NIntegrate wants to evaluate the expressions at different points ) Aug 18, 2014 at 19:34
• @george2079 Good point! I meant this answer as an illustration of my understanding of the cause, not a general purpose solution. Aug 18, 2014 at 19:34

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.

• (+1) I can add that if the value at only one point is needed one can get a bit better performance and small memory requirements by replacing {x, 0, 1} with {x, 1, 1}: NDSolveValue[{y'[x] == foo[x], y[0] == {0, 0}}, y, {x, 1, 1}][1]. Aug 18, 2014 at 19:34
• @AlexeyPopkov Thanks, that's an important point - I'll mention it in my answer.
– Jens
Aug 18, 2014 at 19:42
• Jens, I believe this question is a duplicate. Will your method apply to this example as well?: (34554) Aug 18, 2014 at 19:52
• @Mr.Wizard I believe this is equivalent to the "workaround" in my answer to that question. Aug 18, 2014 at 19:54
• @Michael I missed that. Looks like you covered everything there already. Now your answer just needs more votes. :-) Aug 18, 2014 at 19:55

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.

• Lower level solutions are always of value. +1 Aug 18, 2014 at 19:48
• Consider moving (reposting, then deleting) your answer here: (34554) I believe it is equally applicable, and it would be better placed in a non-duplicate Q&A. Aug 18, 2014 at 19:59