1
$\begingroup$

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}}*)
$\endgroup$
  • 1
    $\begingroup$ Could you explain what the 'Compare plot' is? $\endgroup$ – user21 Jul 14 '16 at 18:00
  • $\begingroup$ How do you account for smaller and larger mesh elements (triangles that have different areas)? $\endgroup$ – Michael E2 Jul 14 '16 at 19:08
  • $\begingroup$ @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). $\endgroup$ – Tom Jul 15 '16 at 2:31
  • $\begingroup$ @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. $\endgroup$ – Tom Jul 15 '16 at 2:32
  • $\begingroup$ @Tom Histogram /@ eigenfunction["ElementMesh"]["MeshElementMeasure"] will show you the distribution of areas. $\endgroup$ – Michael E2 Jul 15 '16 at 2:44
4
$\begingroup$

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["NDSolve`FEM`"];
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}
*)

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

NDSolve`ScaledVectorNorm 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.

|improve this answer|||||
$\endgroup$
  • $\begingroup$ @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. $\endgroup$ – Michael E2 Jul 16 '16 at 19:11
  • $\begingroup$ 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? $\endgroup$ – David Zwicker Jul 7 '17 at 21:31
  • $\begingroup$ Thanks – that did it! $\endgroup$ – David Zwicker Jul 8 '17 at 12:28

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.