# Speed up the computation of many similar numerical integrals over 2D regions embded in 3D

I have a set of domains and a set of integrands. I would like to numerically integrate each integrand over each domain. What is the most efficient way to do this? In my case specifically I have 2D domains embedded in a 3D space.

A minimum working example of the sort of problems I want to solve:

params = RandomReal[{1, 2}, {10, 6}];
doms = Triangle /@ RandomReal[{1, 2}, {10, 3, 3}];
expr[a_, b_, c_, x_, y_, z_] = ((a xp + b yp + c zp)/
Sqrt[(x - xp)^2 + (y - yp)^2 + (z - zp)^2]);
MapThread[NIntegrate[Evaluate[expr @@ #1], {xp, yp, zp} \[Element] #2] &,
Transpose[Tuples[{params, doms}]]]

• Do you need every single of these integrals or just a (weighted) sum of them? Some context could help. Nov 18, 2020 at 16:36
• The solution to each one of these integrals goes into a matrix equation which I solve with LinearSolve. If I solved the linear system using a different technique maybe I could get away with using weighted sums? Nov 18, 2020 at 17:02
• Aha, aha! Sounds like something where you want to use the fast multipole method to speed up the matrix-vector product and to use an iterative solver. Have you ever heard of H-matrices? You should definitely read about them! =) However, the precise structure of expr is very crucial for whether the fast multipole/H-matrix approach will work. Nov 18, 2020 at 17:07
• However, H-matrices are not built into Mathematica. Nov 18, 2020 at 17:08
• That does sound like the right thing to do. I was hoping to write something short and simple enough that it couldn't have any bugs (but still fast enough I could check my solutions to non trivial problems) before a went on to do more complicated, faster stuff. Thank you for the suggestions though, I'll check them out. Nov 18, 2020 at 18:18

Here we can't speed up your code but only simplify your code.

params = Permutations[Range[2, 5], {3}];
doms = Triangle /@ Partition[params, 3];
expr[a_, b_, c_] = ((a xp + b yp + c zp)/Sqrt[xp^2 + yp^2 + zp^2]);
NIntegrate[Evaluate[expr @@ #1], {xp, yp, zp} ∈ #2] &,
Transpose[Tuples[{params, doms}]]];
r2 = Table[
NIntegrate[expr @@ param, {xp, yp, zp} ∈ dom], {param,
params}, {dom, doms}] // Flatten;
r3 = Outer[NIntegrate[expr @@ #1, {xp, yp, zp} ∈ #2] &,
params, doms, 1] // Flatten;
r1 == r2 == r3
(* True *)


You may exploit that you integrals depend only linearly on the parameters. Thus it suffices to compute the integral of {xp, yp, zp}/Sqrt[xp^2 + yp^2 + zp^2] only once on each triangle. On my machine (and without enforced parallelization), this leads to speed-up of factor 8.5:

First@AbsoluteTiming[

NIntegrate[Evaluate[expr @@ #1], {xp, yp, zp} \[Element] #2] &,
Transpose[Tuples[{params, doms}]]];

]

First@AbsoluteTiming[

ints = NIntegrate[
{xp, yp, zp}/Sqrt[xp^2 + yp^2 + zp^2],
{xp, yp, zp} \[Element] #
] & /@ doms;
B = Flatten[Outer[Dot, params, ints, 1]];

]
Max[Abs[(B - A)/A]]

> 1.40043
>
> 0.163504
>
> 3.80123*10^-8

• I'm afraid this is one of those unfortunate situations where I oversimplified the given example. In general, these functions are not linear. See my edits. This is, of course, a great idea in cases where linearity can be used though. Nov 18, 2020 at 16:23
• I was almost sure that this would happen! =D Nov 18, 2020 at 16:37