# 2D interpolation of function with dependent interpolation border

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

Flatten[Parallelize[
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]
Flatten[Parallelize[
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}]
NIntegrate[
f[x, y]*g[x], {x, min, max}, {y, min, Sqrt[x^2 - a^2]}]

• I might be missing it, but I don't see even a toy version of g in your code. – barrycarter Aug 29 '18 at 16:32
• Exp[-1/10 x] in lower line – satoru Aug 29 '18 at 17:02
• 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. – 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]},
SparseArray[
Transpose[{
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.


0.980202

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


0.0000189374

# Interpolating function on unstructured grid

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

Needs["NDSolveFEM"]
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.

• 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)$... – satoru Aug 31 '18 at 10:25
• 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. – Henrik Schumacher Dec 16 '18 at 10:46
• Thanks. Can't the second solution be turned into the first one using zeroth order interpolation? – chris Dec 16 '18 at 17:23