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}}*)
Histogram /@ eigenfunction["ElementMesh"]["MeshElementMeasure"]
will show you the distribution of areas. $\endgroup$ – Michael E2 Jul 15 '16 at 2:44