My problem is to compute $$\int\limits_{x_{min}}^{x_{max}}\int\limits_{y_{min}}^{\sqrt{x^{2}-a^{2}}}W(x,y)g(x)\ dx\ dy$$

Because $W$ have highly oscillatory 2d complex numeric integral inside, dumb Nintgrate will last forever.(It takes 30-50 sec to calculate $W$ in single point). So i decided to compute 2d interpolation of $W$ on rough mesh using the power of parallelization. Naively i tried

Table[{{x, y}, W[x, y]}, {x, min, max, max/32}, {y, min, Sqrt[
        x^2 - a^2], Sqrt[x^2 - a^2]/32}]], 1] 
f = Interpolation[%]

Of course this not work because border of $y$ is not fixed.($W$ blows up if $y>\sqrt{x^2 - a^2}$, i cant put $y_{max}=\sqrt{x_{max}^2 - a^2}$) Additionally i find Parallelize wont work with 2d tables, how to specify parallel computation goes for few y for each fixed x? Is it possible to visualize $W(x,y)$ and how to specify variable y in CountourPlot/Plot3D? Full code with toy $W$ and $g$:

min = 1/10
max = Pi
a = 1/20
g[x_] = Exp[-1/10*x];
W[x_, y_] = Exp[-x*y]
  Table[{{x, y}, W[x, y]}, {x, min, max, max/32}, {y, min, Sqrt[
    x^2 - a^2], Sqrt[x^2 - a^2]/32}]], 1]
f = Interpolation[%]
Plot[NIntegrate[f[x, y], {y, min, Sqrt[x^2 - a^2]}], {x, min, max}]
 f[x, y]*g[x], {x, min, max}, {y, min, Sqrt[x^2 - a^2]}]
  • $\begingroup$ I might be missing it, but I don't see even a toy version of g in your code. $\endgroup$ – barrycarter Aug 29 '18 at 16:32
  • $\begingroup$ Exp[-1/10 x] in lower line $\endgroup$ – satoru Aug 29 '18 at 17:02
  • 6
    $\begingroup$ OK, I see it now. I've edited your code to make it clearer what g is. Please doublecheck to make sure I didn't introduce any errors. $\endgroup$ – barrycarter Aug 29 '18 at 17:34

Quick integration without interpolating function

This discretizes the integration domain into a triangle mesh with maximal triangle circumradius h (the smaller h is the better the approximation of the integral).

h = 2 Pi/40;

min = 1/10
max = Pi
a = 1/20
R0 = ImplicitRegion[y <= Sqrt[x^2 - a^2], {{x, min, max}, {y, min, max}}];
R = DiscretizeRegion[R0, MaxCellMeasure -> {1 -> h}];

Now we compute the values of the integrand on the vertices of the mesh.

vals = ParallelMap[X \[Function] W[X[[1]], X[[2]]] g[X[[1]]], MeshCoordinates[R]];

The vector vals has one entry for each vertex of the mesh. This defines a unique piecewise linear function on the mesh. For computing its integral, we need two things: (i) the triangle areas and (ii) the vertex-triangle adjacency matrix; this is the matrix A which has A[[i,j]] == 1 whenever vertex i is a member of triangle j (and A[[i,j]] == 0 otherwise). This is how to obtain the two and how to compute the integral from them:

A = 
  With[{m = MeshCellCount[R, 2], n = MeshCellCount[R, 0]},
       Join @@ Join @@ MeshCells[R, 2, "Multicells" -> True][[All, 1]],
       Join @@ Transpose[ConstantArray[Range[MeshCellCount[R, 2]], 3]]
       }] -> 1,
    {n, m}
areas = PropertyValue[{R, 2}, MeshCellMeasure];

a = vals.(A.areas)/3.


In our toy model, we can compare with what NIntegrate computes as integral over the discretized domain R:

b = NIntegrate[W[x, y] g[x], {x, y} ∈ R];
a - b


Interpolating function on unstructured grid

Alternatively, you can obtain a second order interpolation function on unstructured grids by employing tools from finite elements:

Rdisc = ToElementMesh[R0, "MeshOrder" -> 2,MaxCellMeasure -> {1 -> h}];
vals = ParallelMap[X \[Function] W[X[[1]], X[[2]]],Rdisc["Coordinates"]];
Winterpol = ElementMeshInterpolation[{Rdisc}, vals];

For example, you can plot the interpolating function with

Plot3D[Winterpol[x, y], {x, y} ∈ Rdisc, PlotPoints -> 50]

These interpolating functions are however infamous for their slow evaluation speed. But if W is as expensive as you say, this may still help a lot.

  • $\begingroup$ Thank you for your work, but this solution focusing on calculation of particular integral(i used it as example) rather than interpolation of $W(x,y)$ - this was the main interest. I need interpolation because i have to do a lot of things with $W(x,y)$... $\endgroup$ – satoru Aug 31 '18 at 10:25
  • $\begingroup$ Hi @chris. I fixed some typos and added some detail. Indeed, I wrote a = (A.areas)/3. instead of a = vals.(A.areas)/3.. No clue how this could happen. Moreover, the syntax highlighting does not know that MeshCells can take options, hence "Multicells" -> True is colored in red. Nontheless, it should work. It was user21 who was so kind to give me the hint that this is possible. $\endgroup$ – Henrik Schumacher Dec 16 '18 at 10:46
  • $\begingroup$ Thanks. Can't the second solution be turned into the first one using zeroth order interpolation? $\endgroup$ – chris Dec 16 '18 at 17:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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