TransWikia.com

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

Mathematica Asked on August 3, 2021

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

2 Answers

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

Answered by cvgmt on August 3, 2021

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[
  
  A = MapThread[
     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

Answered by Henrik Schumacher on August 3, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP