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}]]]
  • $\begingroup$ Do you need every single of these integrals or just a (weighted) sum of them? Some context could help. $\endgroup$ Nov 18, 2020 at 16:36
  • $\begingroup$ 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? $\endgroup$
    – alessandro
    Nov 18, 2020 at 17:02
  • 1
    $\begingroup$ 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. $\endgroup$ Nov 18, 2020 at 17:07
  • 1
    $\begingroup$ However, H-matrices are not built into Mathematica. $\endgroup$ Nov 18, 2020 at 17:08
  • $\begingroup$ 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. $\endgroup$
    – alessandro
    Nov 18, 2020 at 18:18

2 Answers 2


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]);
r1 = MapThread[
   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:

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

  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
  • $\begingroup$ 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. $\endgroup$
    – alessandro
    Nov 18, 2020 at 16:23
  • $\begingroup$ I was almost sure that this would happen! =D $\endgroup$ Nov 18, 2020 at 16:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.