How to calculate the numerical integral more efficiently?

Question

Rencently, I ecountered the following numerical intergral: $$ \begin{cases} \mathbf I_1=\displaystyle \int_{t}\mathbf N' {\mathbf N'}^{\text T}{\rm d}t\\ \mathbf I_2=\displaystyle \int_{t}\mathbf N'' {\mathbf N''}^{\text T}{\rm d}t \end{cases} $$

where,

$$ \begin{cases} \mathbf N'=[N'_{0,p}(t),\cdots,N'_{n,p}(t)]^{\text T}\\ \mathbf N''=[N''_{0,p}(t),\cdots,N''_{n,p}(t)]^{\text T} \end{cases} $$

$N_{i,p}(t)$ is the $i$-th B-spline basis function of degree $p$ and defined on the knot vector $\mathbf T=\{t_0,\cdots,t_{n+p+1}\}$, which could be compute by the built-in BSplineBasis[{p, T}, i, t]. Furthermore, $N'_{i,p}(t), N''_{i,p}(t)$ are the first derivatives and second derivatives of $N_{i,p}(t)$, respectively.

My trial

n = 6;
knots = {0, 0, 0, 0, 0.350812, 0.509446, 0.648472, 1, 1, 1, 1};
Nder1 = D[BSplineBasis[{3, knots}, #, x] & /@ Range[0, n], x];
I1expr = Outer[Times, Nder1, Nder1];
I1 = NIntegrate[I1expr, {x, 0, 1}]

Nder2 = D[BSplineBasis[{3, knots}, #, x] & /@ Range[0, n], {x, 2}];
I2expr = Outer[Times, Nder2, Nder2];
I2 = NIntegrate[I2expr, {x, 0, 1}]

For this simple case $n=6$, I discovered that NIntegrate[] takes about 0.85s and 0.99sfor $\mathbf I_1$ and $\mathbf I_2$, respectively. In additon, it is very intersting to see that the $\mathbf I_1$ and $\mathbf I_2$ are symmetric matrix.

However, when $n=99$, $\mathbf I_1$ will take about 5 min. Obviously, it's very time-consuming! While for the actual case, $n$ ranges from $100$ to $1000$.

So I would like to know:

• Is there good strategy to speed up it?

---

Update

Some strategies I could come up with:

(1) For the $\mathbf I_1$, which is a matrix with dimensions $\{n+1,n+1\}$

$$I_{i,j}=\int_tN'_{i,p}(t)N'_{j,p}(t){\rm d}t\Rightarrow I_{i,j}=I_{j,i}$$

Namely, $\mathbf I_1$ is a symmetric matrix, we just need to calculate the lower triangular. $$ \begin{pmatrix} I_{0,0}\\ I_{1,0} & I_{1,1}\\ \vdots & \vdots & \ddots\\ I_{n,0} & I_{n,1} & \cdots & I_{n,n} \end{pmatrix} $$

(2)The first derivatives $N'_{i,p}(t)$ has the following identity:

$$ N'_{i,p}(t)=p\left( \frac{N_{i,p-1}(t)}{t_{i+p}-t_i}-\frac{N_{i+1,p-1}(t)}{t_{i+p+1}-t_{i+1}} \right) $$

(3) Local supporting property

If $t\notin [t_i,t_{i+p+1})$, then $N_{i,p}(t)=0$. So I guess some items must be equal to $0$. i.e., there is no need to calculate them painfully.

According to (1) ~ (3), the matrix $\mathbf I_1$ should be a an $(n+1) \times (n+1)$ symmetric banded matrix as shown below.

Answer 1

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

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 *)