# NIntegrate fails when integrating over a list from an external (MathLink) function

NIntegrate fails when integrating over a list from an external (MathLink) function. For simplicity, consider an external function f[x] that returns the list {x,2x}. In Mathematica, the function would be defined as f[x_]={x,2x}. In reality I have a much more complicated function that cannot be implemented inside Mathematica. NIntegrate works fine if the external function f[x] returns a single number, not a list.

Install["f"]
f[.5]
(* {0.5, 1.} *)
NIntegrate[f[x], {x, 0, 1}]
(* Integrand f[x] is not numerical at {x} = {0.00795732}. >>*)
(* Integrand f[x] is not numerical at {x} = {0.00795732}. >>*)


This is the file f.tm that defines the external function f[x]:

#include "mathlink.h"

:Begin:
:Function:       f
:Pattern:        f[x_?NumericQ]
:Arguments:      {x}
:ArgumentTypes:  {Real}
:ReturnType:     Manual
:End:

void f(double x)
{
double y;
y = x;
y = 2 * x;

}

int main(int argc, char* argv[])
{
return MLMain(argc, argv);
}

• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey May 9 '15 at 19:12

NIntegrate can handle vector integrands if they are explicitly in vector form:

f[x_] := {x, 2 x}
NIntegrate[f[x], {x, 0, 1}]
(* {0.5, 1.} *)


When the function is a black box which will only evaluate for numeric arguments, NIntegrate no longer sees the list structure and complains about not getting a scalar value:

g[x_?NumericQ] := {x, 2 x}
NIntegrate[g[x], {x, 0, 1}]
(* NIntegrate::inum: Integrand g[x] is not numerical ... *)


We can still integrate g by supplying the components separately:

g1[x_?NumericQ] := g[x][]
g2[x_?NumericQ] := g[x][]
NIntegrate[{g1[x], g2[x]}, {x, 0, 1}]
(* {0.5, 1.} *)


f[x] should return a single value not a pair. Your example data is of the form

data = {#, 2 #} & /@ Range[0, 1, .05];

Clear[f]

f = Interpolation[data];


Note that f[x] is single-valued

f[.5]


1.

Using NIntegrate

NIntegrate[f[x], {x, 0, 1}]


1.

However, you can Integrate an InterpolatingFunction

Integrate[f[x], {x, 0, 1}]


1.

NIntegrate[f[x], {x, 0, 1}] ==
Integrate[f[x], {x, 0, 1}] ==
Integrate[2 x, {x, 0, 1}]


True