NIntegrate receiving a function with NumericQ which returns a list

Question

Can some one please explain why

 f[x_?NumericQ] := {{x, 2}, {-1, -x}}
 NIntegrate[ f[S], {S, 0, 1}]
 (*->NIntegrate::inum: Integrand f[S] is not numerical at {S} = {0.00795732}. >> *)

A possible work around is

f1[x_?NumericQ][part__] := {{x, 2}, {-1, -x}}[[part]]

then use something like

ArrayReshape[ NIntegrate[ f1[S][#[[1]], #[[2]]] & /@  {{1, 1}, {1, 2}, {2, 1}, {2, 2}}, {S, 0, 1}], {2, 2}]

but this is ugly, and does it evaluate f1 four times?

This type of problem appears when f involves matrix calculations and a numerical component (can not be done symbolically). This way, one is always forced to calculate all elements, even if you need only one. So it would advantages if NIntegrate did not evaluate f more times than necessary.

I thank any assistance in advance.

Possible related questions 1, 2 and 3.

Answer 1

Perhaps someone who knows can address the actual limitations of NIntegrate. It seems to me that, despite a couple of examples in the documentation, NIntegrate does not integrate vectors, matrices, and other arrays as arrays per se. Instead it integrates their components individually, effectively mapping NIntegrate onto the components.

The main evidence is the following. If I make f[S] return four component expressions, NIntegrate gets mapped onto each one. Note the f[1] is evaluated first, which would tell NIntegrate that the value of f is a 2x2 array. In the OP's code, f[S] is a single expression, which NIntegrate does not like.

ff[x_?NumericQ] := {{x, 2}, {-1, -x}};
ff[x_] = Array[ff[x, ##] &, {2, 2}];
f[x_] := (Print[ff@x]; ff[x]);
NIntegrate[f[S], {S, 0, 1}]

{{1, 2}, {-1, -1}}

{{ff[S,1,1], ff[S,1,2]}, {ff[S,2,1], ff[S,2,2]}}

NIntegrate::inumr: The integrand ff[S,1,1] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1}}. >>
...
General::stop: Further output of NIntegrate::inumr will be suppressed during this calculation. >>

(* {{NIntegrate[ff[S, 1, 1], {S, 0, 1}], 
     NIntegrate[ff[S, 1, 2], {S, 0, 1}]},
    {NIntegrate[ff[S, 2, 1], {S, 0, 1}],
     NIntegrate[ff[S, 2, 2], {S, 0, 1}]}} *)

I don't claim this is conclusive evidence, but I have not found a way to use NIntegrate to integrate the OP's f.

Workaround

NIntegrate and NDSolve are not equivalent, not in functionality, speed, or accuracy. In this case, NDSolve can address the OP's problem on it's own terms -- that is, it can integrate a matrix differential equation in terms of matrices.

ClearAll[a, f];
f[x_?NumericQ] := {{x, 2}, {-1, -x}};

asol = NDSolveValue[{a'[t] == f[t], a[0] == {{0, 0}, {0, 0}}}, a, {t, 0, 1}];
asol[1]

(* {{0.5, 2.}, {-1., -0.5}} *)

It seems (to me) that if the component functions of the matrix function f could be integrated separately, then NIntegrate is likely to do a superior job finding the integrals. But when the components cannot be integrated independently, one can use NDSolve.