Return to Revisions

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 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", "How to avoid repetitive calculation when doing numerical integral?".

For more details how rules like ArrayOfFunctionsRule are implemented see "How to implement custom integration rules for use by NIntegrate?".

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"]

Convert the integrand matrix into a matrix of functions:

I1exprFS = 
  Map[Function[{fx}, Function[Evaluate[fx /. x -> #]]], I1expr, {2}];

Call NIntegrate with the new rule:

res1 =
    NIntegrate[1, {x, 0, 1}, 
      Method -> {"GlobalAdaptive", "SingularityHandler" -> None, 
        Method -> {ArrayOfFunctionsRule, "Functions" -> I1exprFS}}]; // 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} *)

Converting to functions:

funcs = 
  Map[Function[{fx}, Function[Evaluate[fx /. x -> #]]], 
   funcsExpr, {2}];

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

res1 = NIntegrate[1, {x, 0, 1}, 
   Method -> {"GlobalAdaptive", "SingularityHandler" -> None, 
     Method -> {ArrayOfFunctionsRule, 
       "Functions" -> funcs}}]; // 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 *)