# Approximating Integral under 2D Interpolation function

First, I have already been working on this problem, and have asked a similar question previously here.

I am trying to find the integrals of 2D interpolation functions that are output from NDEigensystem[]. I will be doing this integral hundreds if not thousands of times, so it is necessary to speed it up as much as possible. I do not care if the value is a scalar multiple away from the most exact value produced by NIntegrate[] and it can be (slightly) off from the a scalar multiple of the most exact value. Essentially, over many different functions the more exact value and the quick approximation should have a roughly linear relationship.

Right now the fastest NIntegrate[] function that I'm using is

fastint[eigenfunction_]:=Abs[NIntegrate[
eigenfunction[x,y],{x,y}\[Element]eigenfunction["ElementMesh"],
Method->{Automatic,"SymbolicProcessing"->0}]];


And the faster approximation I am using is just a sum over the "ValuesOnGrid" for the interpolation function

fastapprox[eigenfunction_]:=Abs[Total[eigenfunction["ValuesOnGrid"]]];


These two functions form a roughly linear relationship when the eigenfunctions are determined by a highly symmetric region (e.g. Disk[]). But when I do this integral with an asymmetric region, I obtain a plot similar to:

which is obviously not linear.

I feel like I am forgetting to incorporate something, but I do not know what. Again, the end goal is to either speed up NIntegrate[] or develop an approximation method that is faster while maintaining the linear relationship with a more exact value.

Any help is much appreciated, thank you.

Edit: I have developed code for a basic Reimannian approximation that is much faster than NIntegrate[] and still is able to get the results I am hoping for:

reimanapprox[region_,numberdiv_:50][functionlist_]:=
Module[{Dx, Dy, area, squaregrid, i, inside, heights, xmin, xmax, ymax, ymin, coord},
coord = region["Coordinates"];
xmax = Max[coord[[All, 1]]];
xmin = Min[coord[[All, 1]]];
ymax = Max[coord[[All, 2]]];
ymin = Min[coord[[All, 2]]];
Dx = (xmax-xmin)/numberdiv;
Dy = (ymax-ymin)/numberdiv;
area = Dx*Dy;
squaregrid = Flatten[Table[{x, y},
{x, xmin, xmax, Dx},
{y, ymin, ymax, Dy}], 1];
inside = {};
For[i = 1, i <= Length[squaregrid], i++,
If[RegionMember[region, squaregrid[[i]]],
AppendTo[inside, squaregrid[[i]]],
False]];
heights = Table[
Map[eigenfunctions[[i]][#[[1]], #[[2]]] &, inside],
{i, Length[eigenfunctions]}];
area*Map[Total, heights]]


Which then can do many of these integrals much faster than NIntegrate[] (efunc is just a list of (30) interpolation functions that I am testing which are all over the region eregion).

RepeatedTiming[fastapprox[eregion][efunc]

(*Returns: {1.00, {-0.86141, -0.129498, -0.194252, 0.357529, 0.247833, 0.0704943, 0.0505611, ..., 0.00346365}}*)

RepeatedTiming[Table[fastintegral[efunc[[i]]], {i, Length[efunc]}]]

(*Returns: {21.1, {-0.861569, -0.129423, -0.194507, 0.357675, 0.248136,0.0706918, 0.050766,..., 0.00274389}}*)

• Could you explain what the 'Compare plot' is? – user21 Jul 14 '16 at 18:00
• How do you account for smaller and larger mesh elements (triangles that have different areas)? – Michael E2 Jul 14 '16 at 19:08
• @user21 The 'Compare plot' is a ListPlot whose points have the precise answer as the x-value and the approximation as the y-value. (If the approximation was good enough, this plot would look roughly linear). – Tom Jul 15 '16 at 2:31
• @MichaelE2 I do not account for different sized mesh elements. I was hoping that they would all be close enough to the same size that accounting for them would simply amount to multiplying each term by the same factor. – Tom Jul 15 '16 at 2:32
• @Tom Histogram /@ eigenfunction["ElementMesh"]["MeshElementMeasure"] will show you the distribution of areas. – Michael E2 Jul 15 '16 at 2:44

Well, since an example NDEigensystem was not provided, I went ahead an took one from the docs.

{vals, funs} = NDEigensystem[
{-Laplacian[u[x, y], {x, y}],
DirichletCondition[u[x, y] == 0, True]},
u,
{x, y} ∈ RegionDifference[Cuboid[{-3, -3/2}, {0, 3/2}], Disk[]],
10,
Method -> {"SpatialDiscretization" ->
{"FiniteElement", {"MeshOptions" -> {"MaxCellMeasure" -> {"Area" -> 0.005}}}}}];


We can compare an "accurate" NIntegrate[] with a quadratic-order triangle rule. Warning: The triangle rule assumes the quadratic points of the element mesh are the midpoints of the sides of the triangles. This may not be true everywhere -- on a curved boundary, especially. In the present example, it can be seen not to introduce a large error.

(accurate = Table[
NIntegrate[
eigenfunction[x, y], {x, y} ∈ eigenfunction["ElementMesh"],
Method -> {Automatic, "SymbolicProcessing" -> 0}],
{eigenfunction, funs}]; "accurate") // RepeatedTiming

Needs["NDSolveFEM"];
emesh = funs[[1]]["ElementMesh"];
indices = emesh["MeshElements"] /.
{TriangleElement[t_]} :> Flatten@t[[All, 4 ;; 6]]; (* midpoints of sides *)
(fast = Table[
Total[
Total[Partition[eigenfunction["ValuesOnGrid"][[indices]], 3], {2}] *
First@eigenfunction["ElementMesh"]["MeshElementMeasure"] / 3
],
{eigenfunction, funs}]; "fast") // AbsoluteTiming
(*
{2.23, "accurate"}
{0.011379, "fast"}
*)


Okay, so it's fast. Let's now examine its accuracy.

fast
accurate
(*
{2.08897, -3.87294*10^-7, -0.310354, 0.749835, -0.0000232558,
0.47398, 6.61989*10^-7, -4.84232*10^-6, -0.679196, -0.290395}

{2.08898, -3.91154*10^-7, -0.310374, 0.749833, -0.0000236342,
0.473997, 2.35524*10^-6, -3.86783*10^-6, -0.679223, -0.290408}
*)

NDSolveScaledVectorNorm[Infinity, {10^-4, 10^-5}][fast - accurate, accurate]
(*  0.483535  *)


NDSolveScaledVectorNorm returns an error scaled according the the {precision_goal, accuracy_goal} to be less than 1 if the error fast - accurate meets the goals, as defined by Mathematica. In short, the integrals above meet the goals corresponding to PrecisionGoal -> 4 and AccuracyGoal -> 5.

• @Tom You're welcome. -- Note: if the domain is a rectangle and the mesh element type QuadElement, then a different integration rule would be needed. Mean[eigenfunction["ValuesOnGrid"]] times the area might give a decent approximation. The weights wouldn't quite be right for the different kinds of mesh points, but it would be very fast, at the least. – Michael E2 Jul 16 '16 at 19:11
• This does not seem to work on Mathematica 11, where I receive the error Part::pkspec1: The expression {NDSolveFEMTriangleElement[{<<2291>>}]} cannot be used as a part specification. Is there any way to adapt this to version 11? – David Zwicker Jul 7 '17 at 21:31
• Thanks – that did it! – David Zwicker Jul 8 '17 at 12:28