# Eigenfunctions of the Laplacian on an arbitrary mesh

So, I've constructed a mesh over which I'd like to find eigenfunctions of Laplace's equation with a free boundary (a zero Neumann boundary condition along the edge).

Mostly because I figured an electric guitar–shaped Chladni plate would look really cool (but also because I'd like to try to find the resonant frequencies of a lamina this shape).

Normally, to do this sort of eigenvalue problem, I'd do something like this:

ParametricNDSolve[
{Laplacian[phi[x, y], {x, y}] == -w^2*phi[x, y]},
{x, y} ∈ mesh,
w
]


with some boundary conditions in there.

However, while I feel like I understand (well enough) the theory behind what I want to accomplish (at the very least, I can derive the eigenfunctions in coordinate systems where Laplace's equation is solvable by separation of variables), I don't quite know how to go about implementing it in this case.

How do I impose a free boundary for an arbitrary mesh?

This is relatively easy on a rectangle or circle, but on this mesh, I think would have to somehow determine the vector normal to the mesh boundary at each point (so that I can specify $\mathbf{n}\cdot\nabla \phi = 0$). Furthermore, once I've specified the Neumann boundary condition, are there any other boundary conditions I have to keep in mind (I imagine that if I only specify a zero derivative on the edge, NDSolve will give me a constant solution)? Thank you so much in advance for your help, and I promise to post cool plots of irregularly-shaped guitar Chladni plates when I'm done.

Edit: Alright, I guess the question I was asking wasn't the question I should have been asking. To boil it down, we don't really have to worry about the Neumann boundary condition in this case, as it's zero by default. However, Mathematica doesn't actually have built-in eigenvalue analysis for these sorts of problems (thus the necessity of using some low-level code, as utilized to great effect by user21 here).

• Have you seen this: Numerically solving Helmholtz equation in 2D for arbitrary shapes? Not saying it's a duplicate, just pointing out the close relation. – Jens May 28 '15 at 4:57
• Cool! I would help if I could. I am looking forward to the Chladni plates though :-) – MarcoB May 28 '15 at 4:58
• I guess what you need is NeumannValue – Jens May 28 '15 at 5:03
• The Helmholtz example looks like it'll at least be very valuable in the endgame. I'm going to try modifying the code a little to see whether it works with NeumannValue instead of DirichletCondition. – Michael L. May 28 '15 at 5:22
• What you need is a Neumann zero boundary value. Since that is the default, you do not need to give any boundary condition at all. – user21 May 28 '15 at 6:42

As promised:

And here are the nodes:

Now to outline the process. [This first part has no Mathematica.] First, I found an image of a guitar using Google's Image Search. I then went into GIMP (although you can use any image editing software, or even draw the image yourself) and used it as a template to create a silhouette (it's okay if the edges are a little rough). Then, I created a parametric plot using this technique.

img = Binarize[Import["/home/michael/Downloads/test.jpg"]~ColorConvert~"Grayscale"~ImageResize~500~Blur~3]~Blur~3;
lines = Cases[Normal@ListContourPlot[Reverse@ImageData[img], Contours -> {0.5}], _Line, -1];

param[x_, m_, t_] := Module[{f, n = Length[x], nf},
f = Chop[Fourier[x]][[;; Ceiling[Length[x]/2]]]; nf = Length[f];
Total[Rationalize[2 Abs[f]/Sqrt[n] Sin[Pi/2 - Arg[f] + 2. Pi Range[0, nf - 1] t],.01][[;;Min[m, nf]]]]]

tocurve[Line[data_], m_, t_] := param[#, m, t] & /@ Transpose[data]

parplot = ParametricPlot[Evaluate[tocurve[#, 40, t] & /@ lines], {t, 0, 0.998},
Frame -> True, Axes -> False, PlotPoints -> 3] /.
Line[l_List] :> {{Blue, Polygon[l]}, {White, Line[l]}}


I recommend setting PlotPoints to something very low in order to decrease the need for mesh refinement near the edges of your lamina (although this will decrease the resolution of your boundary curve; I did not do this in my original question, which gave me a mesh that was far too complex near the edge). Also, closing the loop (letting t go to 1 in our plot bounds) was for some reason giving me a weird irregularity (jaggedness of the boundary) that was almost invisible but was causing the mesh to become extremely (and stubbornly) fine near that point.

Then, use DiscretizeGraphics[] and ToElementMesh[] to convert this polygon to an element mesh.

g = DiscretizeGraphics[parplot]


<<NDSolveFEM
mesh = ToElementMesh[g, MeshQualityGoal -> 0.8, MaxCellMeasure -> 30, "MeshOrder" -> 1]
mesh["Wireframe"]


Using the answer provided here by user21 (check out his answer for a better explanation of what each piece of this code does), we can then find the eigenfunctions of the Helmholtz differential equation (which is the eigenvalue equation for the Laplacian) over our lamina.

pde = D[u[t, x, y], t] - Laplacian[u[t, x, y], {x, y}] + u[t, x, y] == 0;
Γ = DirichletCondition[u[t, x, y] == 0, True];

nr = ToNumericalRegion[mesh];
{state} = NDSolveProcessEquations[{pde, Γ, u[0, x, y] == 0}, u, {t, 0, 1}, {x, y} ∈ nr];

femdata = state["FiniteElementData"]
initBCs = femdata["BoundaryConditionData"];
methodData = femdata["FEMMethodData"];
initCoeffs = femdata["PDECoefficientData"];

vd = methodData["VariableData"];
sd = NDSolveSolutionData[{"Space" -> nr, "Time" -> 0.}];

discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

stiffness = discretePDE["StiffnessMatrix"];
damping = discretePDE["DampingMatrix"];

nDiri = First[Dimensions[discreteBCs["DirichletMatrix"]]];
numEigenToCompute = 10;
numEigen = numEigenToCompute + nDiri;

res = Eigensystem[{stiffness, damping}, -numEigen];
res = Reverse /@ res;
eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]];
eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]];

evIF = ElementMeshInterpolation[{mesh}, #] & /@ eigenVectors;


The above plots can then be created with this:

densityplot[i_] := DensityPlot[Evaluate[evIF[[i]][x, y]], {x, y} \[Element] mesh,
PlotPoints -> 256, PlotRange -> All, ColorFunction -> "TemperatureMap", PlotLabel -> i]
nodeplot[i_] := ContourPlot[Evaluate[evIF[[i]][x, y]] == 0, {x, y} \[Element] mesh, PlotPoints -> 256, PlotRange -> All, PlotLabel -> i]

Show[GraphicsGrid[Table[{densityplot[2i-1], densityplot[2i]}, {i, 1, 5}]], ImageSize -> Full]
Show[GraphicsGrid[Table[{nodeplot[2i-1], densityplot[2i]}, {i, 1, 5}]], ImageSize -> Full]


Thanks so much to user21 for his brilliant code.

• I wanted to show the work that I was able to do using a reference provided by another user, as well as make good on my earlier promise to post Chladni plates once I could find them. I assume that somebody else may have a better method or explanation, or be able to answer my question about resonant frequencies, so I have no intention of accepting my own answer. Would this be better as an addendum to my question instead? – Michael L. May 28 '15 at 15:53
• Sure, I'll write something more explanatory soon! – Michael L. May 28 '15 at 15:59
• Is that better? I also decided that I could use a finer mesh and a more refined boundary, and I think it gave me nicer-looking results too. – Michael L. May 28 '15 at 16:49
• doesn't eigenValues hold the resonant frequencies ? – chris May 28 '15 at 17:30
• Yes, you're right, that was a silly question. I'm concerned though that using this method, all of the resonant frequencies are so close together for the guitar model (mostly within a few hundredths or thousandths of 1). – Michael L. May 28 '15 at 18:15