Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

This question already has an answer here:

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}*)
share|improve this question

marked as duplicate by Michael E2, Mr.Wizard Aug 18 '14 at 19:53

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

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. – Mr.Wizard Aug 18 '14 at 18:57
Related: (11772), (14645), (21662), and even: (1731) – Mr.Wizard Aug 18 '14 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. – Mr.Wizard Aug 18 '14 at 19:57

3 Answers 3

up vote 4 down vote accepted

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}]}
share|improve this answer
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. – george2079 Aug 18 '14 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 ) – george2079 Aug 18 '14 at 19:34
@george2079 Good point! I meant this answer as an illustration of my understanding of the cause, not a general purpose solution. – Mr.Wizard Aug 18 '14 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.


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.

share|improve this answer
(+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]. – Alexey Popkov Aug 18 '14 at 19:34
@AlexeyPopkov Thanks, that's an important point - I'll mention it in my answer. – Jens Aug 18 '14 at 19:42
Jens, I believe this question is a duplicate. Will your method apply to this example as well?: (34554) – Mr.Wizard Aug 18 '14 at 19:52
@Mr.Wizard I believe this is equivalent to the "workaround" in my answer to that question. – Michael E2 Aug 18 '14 at 19:54
@Michael I missed that. Looks like you covered everything there already. Now your answer just needs more votes. :-) – Mr.Wizard Aug 18 '14 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.

share|improve this answer
Lower level solutions are always of value. +1 – Mr.Wizard Aug 18 '14 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. – Mr.Wizard Aug 18 '14 at 19:59

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