## General 

When given an array of integrands `NIntegrate` is run separately over each array element. That is not necessary though, the core `NIntegrate` integration strategies can work with any integrands as long as the error estimates are real numbers.

The motivation for implementing [`ArrayOfFunctionsRule`](https://github.com/antononcube/MathematicaForPrediction/blob/master/Misc/ArrayOfFunctionsRule.m) is to provide a significant speed-up for integrands that are arrays of functions. That is achived by evaluating all functions with the same integration rule abscissas and weights. (Depending on array sizes between 10 and 100 times speed-up is achieved.)

As mentioned in the comments these are related posts/answers: ["NIntegrate over a list of functions"](http://mathematica.stackexchange.com/a/81436/34008), ["How to avoid repetitive calculation when doing numerical integral?"](http://mathematica.stackexchange.com/a/120737/34008). 

For more details how rules like `ArrayOfFunctionsRule` are implemented see ["How to implement custom integration rules for use by NIntegrate?"](http://mathematica.stackexchange.com/q/118324/34008).

## Performance comparison 

The code below uses the definitions in the question.

Load the code for the integration rule `ArrayOfFunctionsRule`:

    Import["https://raw.githubusercontent.com/antononcube/\
    MathematicaForPrediction/master/Misc/ArrayOfFunctionsRule.m"]

   
Call `NIntegrate` with the new rule:

    res1 =
        NIntegrate[1, {x, 0, 1}, 
          Method -> {"GlobalAdaptive", "SingularityHandler" -> None, 
            Method -> {ArrayOfFunctionsRule, "Functions" -> I1expr}}]; // AbsoluteTiming

    (* {0.021779, Null} *)

Compare with the standard `NIntegrate` call:

    I1 = NIntegrate[I1expr, {x, 0, 1}]; // AbsoluteTiming    

    (* {0.503172, Null} *)

We see that `ArrayOfFunctionsRule` provides 20 times speed-up (for the functions defined in the question.)   

Verify agreement of the results:

    Norm[res1 - I1, 2]

    (* 4.39568*10^-7 *)

See the options of `ArrayOfFunctionsRule`. With the option "ErrorsNormFunction" different norms can be used to compute the integration errors.

Note that the rule has to be used with a strategy specification that has the option "SingularityHandler" -> None, and that the rule does not work correctly with ranges that have infinity.

### Larger matrices

This makes a larger matrix (based on the integrands in the question):

    funcsExpr = I1expr;
    funcsExpr = Table[funcsExpr, {100}, {10}];
    funcsExpr = Partition[Flatten[funcsExpr], 100];
    Dimensions[funcsExpr]    
    (* {490, 100} *)

This integration is ~ 100 times faster than the default options one:

    res1 = NIntegrate[1, {x, 0, 1}, 
       Method -> {"GlobalAdaptive", "SingularityHandler" -> None, 
         Method -> {ArrayOfFunctionsRule, 
           "Functions" -> funcsExpr}}]; // AbsoluteTiming
    
    (* {5.4292, Null} *)

Integration with `NIntegrate`'s default method:

    I1 = NIntegrate[funcsExpr, {x, 0, 1}]; // AbsoluteTiming
    
    (* {538.901, Null} *)

Adherence verification:

    Norm[res1 - I1, 2]    
    (* 6.31739*10^-6 *)