I have been working with NDEigensystem in order to find the resonant frequencies of different shaped drums. I would also like the coefficients of the different eigenfunctions so that each of these frequencies has a corresponding intensity (the coefficients squared). The only problem is that the NIntegrate command with my two interpolating functions is running really slow. First, here is my function to model a struck drum:

pluck[{{f_, g_}, {umin_, umax_}}, Nbordpts_][{Ax_, Ay_}] := 
  Module[{hmax, d, l, recurt, constructxy, oneline, datapts, recurtpoints, mesh},
  hmax = 1;
  d = .05;
  (* l:[0,1]->R^2 *)
  l[{x1_, y1_}, {x2_, y2_}][t_] := {(1 - t) x1 + t x2, (1 - t) y1 + t y2};
  (* Constructs points for a line from (0,0) to (1,hmax) *)
  recurt[{n_, t_, z_}] := {n + 1,
   1 - (1 - t)/(hmax + d/(n + 1) - z) (2 d)/(n + 1),
   hmax - d/(n + 1)};
  recurtpoints = Join[NestList[recurt, {0, 0, 0}, 10],
   {0, (1 - (2 d)/(hmax + d)) i/4, (hmax + d) (1 - (2 d)/(hmax + d)) i/4},
   {i, 1, 3}]];
  (* constructxy: constructs points in R^3 from points recursively determined previously
  in recurt using a border point calculated from coordinate function heads f,g *)
  constructxy[U_][{n_, t_, z_}] := {l[{f[U], g[U]}, {ax, ay}][t], -z};
  (* oneline: function to create set of points {{x,y},z} along one line from
  {f[U],g[U]} to {ax,ay} - f,g need to be the *heads* of the coordinate functions
  e.g Cos,Sin NOT Cos[u],Sin[u] *)
  oneline[Num_][U_] := (Map[constructxy[U], recurtpoints]);
  (* Essentially create these interpolation points for a bunch of lines from a
  bunch of different border points and put them all into one set, then
  interpolate hmax is Absolute value of maximum depth (default is 1)*)
  datapts = (DeleteDuplicatesBy[Append[Flatten[Map[oneline[10],
   Range[umin, umax, (umax - umin)/Nbordpts]], 1], {{ax, ay}, -hmax}],
   First] /. {ax -> Ax, ay -> Ay});
  (* This requires a random addition of a small amount to make sure that
  Mathematica does not mistake a triangle in the mesh as having zero area *)
  mesh = NDSolve`FEM`ToElementMesh[datapts[[All, 1]] + RandomReal[{-10^-9, 10^-9},
  NDSolve`FEM`ElementMeshInterpolation[{mesh}, datapts[[All, 2]], InterpolationOrder -> 1]];

This function creates a list of points in the form {{x,y},z} and then constructs a first order interpolation of them.

{circlevals,circlefuncs}=NDEigensystem[{-Laplacian[f[x, y], {x, y}],
 DirichletCondition[f[x, y] == 0, True]}, {f}, {x, y} \[Element]Disk[], 40];
circlepluck=pluck[{{Cos,Sin},{0,2 Pi}},100][{0,0}];

Now the list {circlefreqs,circlefuncs} contains the eigenfrequencies and the eigenfunctions for the region of the unit disk, and circlepluck contains an interpolating function that models a circular drum struck at its center. Here is the command that consistently takes over two minutes (everything else is fast on its own):

 NIntegrate[circlepluck[x, y]*circlefuncs[[i]][x, y], {x, y} \[Element]Disk[],
 {i, 1, Length[circlefuncs]}]

I have tried changing the interpolation order in NDEigensystem as well as many of its other options such as making the discretizing of the region more precise and turning off VectorNormalization. Also, I have tried different methods for NIntegrate (like "Trapezoidal", "MonteCarlo", among others) and changing the precision/accuracy goals, but all of them are either too far away from the actual values or too slow.

I do not need the values to be incredibly accurate, but I would like their relative ratios to stay somewhat accurate (To explain slightly more broadly: I will be running this code for many different borders and I would like to find out which border sounds the best).

Any help would be much appreciated, thank you.

  • $\begingroup$ Try to extract the mesh from the eigenfunction and use that for integration. Something like: NIntegrate[ circlepluck[x, y]*circlefuncs[[i]][x, y], {x, y} \[Element] circlefuncs[[1]]["ElementMesh"]] - this would switch off the adaptive mesh refinement. $\endgroup$
    – user21
    Jun 27, 2016 at 17:24
  • $\begingroup$ Also, have a look at the "VectorNormalization" option of NDEigensystem. $\endgroup$
    – user21
    Jun 27, 2016 at 17:25
  • $\begingroup$ Related: "How to integrate functions of linearly interpolated data?" $\endgroup$ Jun 27, 2016 at 18:06
  • $\begingroup$ Extracting the mesh from the eigenfunction did speed it up a ton, thank you. Do you have any suggestions on what to do with the "VectorNormalization" option? From what I am reading, it looks like I can define my own function to be used to normalized, is that correct? $\endgroup$
    – Tom
    Jun 27, 2016 at 18:07
  • $\begingroup$ @Ton, yes that is correct. If you are checking the integral of the eigenfunctions, it might be possible to tweak the "VectorNormalization" such that NDEigensystem returns the eigenfunctions normalized to the integral. I am not sure if this is really possible and if this is really important. Just wanted to point it out. $\endgroup$
    – user21
    Jun 27, 2016 at 18:28

1 Answer 1


Making my comment an answer:

Try to extract the mesh from the eigenfunction and use that for integration. Something like:

NIntegrate[ circlepluck[x, y]*circlefuncs[[i]][x, y], {x, y} \[Element] circlefuncs[[1]]["ElementMesh"]]

This would switch off the adaptive mesh refinement.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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