The basic problem
I tried to simplify the question down to the following. My apologies if I oversimplified it. The basic question is why does NIntegrate[] succeed on this:
NIntegrate[{t^2, t^3}, {t, 0, 1}]
{0.333333, 0.25}
But fail on this:
integrand[t_?NumericQ] := {t^2, t^3};
NIntegrate[integrand[t], {t, 0, 1}]
NIntegrate::inum: Integrand integrand[t] is not numerical at {t} = {0.00795732}.
NIntegrate::inum: Integrand integrand[t] is not numerical at {t} = {0.00795732}.
NIntegrate[integrand[t], {t, 0, 1}]
NIntegrate threads over vectors but does not integrate vector-valued functions
Well, the "why not" question could be answered only by the design team at WRI, but it appears the choice was for NIntegrate to integrate scalar-valued functions only. When the integrand expression is literally a List[], then NIntegrate[] is threaded over the list. It is in effect "listable" over the first argument. NIntegrate will integrate the elements of the list if they are scalar-valued. On the other hand, when the argument has a head other than List, NIntegrate expects the value to be a scalar. If it is not, it complains and gives up.
One may inspect the thread with Trace or TracePrint:
TracePrint[
NIntegrate[{t, {t^2, t^3}}, {t, 0, 1}],
TraceDepth -> 2
]
NIntegrate[{t,{t^2,t^3}},{t,0,1}]
NIntegrate
{NIntegrate[t,{t,0,1}],NIntegrate[{t^2,t^3},{t,0,1}]}
List
NIntegrate[{t^2,t^3},{t,0,1}]
{NIntegrate[t^2,{t,0,1}],NIntegrate[t^3,{t,0,1}]}
{0.333333,0.25}
{0.5,{0.333333,0.25}}
Out[..]=
{0.5, {0.333333, 0.25}}
One might also infer that NIntegrate is evaluated independently on each element of a list from the following:
NIntegrate[{t^2, integrand[t]}, {t, 0, 1}]
NIntegrate::inum: Integrand integrand[t] is not numerical at {t} = {0.00795732}.
NIntegrate::inum: Integrand integrand[t] is not numerical at {t} = {0.00795732}.
{0.333333, NIntegrate[integrand[t], {t, 0, 1}]}
Possible workaround
One may integrate a vector-valued function with NDSolve:
Block[{a = 0, b = 1},
NDSolveValue[{y'[t] == integrand[t], y[a] == 0 integrand[0]},
y[b], {t, a, b}]
]
{0.333333, 0.25}
Well, 0 integrand[0] is a slick way to get a zero vector of the same dimension as integrand[t]. One might prefer to write it explicitly as {0, 0} when it is known that the integrand is of dimension 2, say.
One difference is that NDSolve controls error locally whereas NIntegrate uses a global estimate of the error. Other differences are that NIntegrate has sophisticated handling of singularities and integrates multivariate functions.
