# Numerically solving Helmholtz equation in 2D for arbitrary shapes

I would like to solve the Helmholtz equation with Dirichlet boundary conditions in two dimensions for an arbitrary shape (for a qualitative comparison of the eigenstates to periodic orbits in the corresponding billiard systems):

$\Omega =$ some boundary e.g. a circle, a regular polygon etc.

$\nabla^2 u(x,y) + k^2u(x,y) =0 \quad x,y \in \Omega \\ u(x,y) = 0 \quad x,y \in \partial\Omega$

There seems already to be a solution here; unfortunately, I'm not experienced enough with Mathematica to extract a minimal working example from this code and adapt it to my needs.

Is there a way to use NDSolve (or ParametricNDSolve) to calculate the eigenvalues k and corresponding eigenstates for the problem?

As far as I understand Mathematica can handle FEM and finite difference methods, which could be used to solve this kind of equations?

• Look into What's New in M10 for a start. Jul 31, 2014 at 14:53
• If you're looking for insights about the correspondence between the billiard-ball problem and the Helholtz and Schrödinger equations, a great place to go is David Tannor's Introduction to Quantum Mechanics: a Time-Dependent Perspective. Aug 1, 2014 at 13:44
• Version 11 has both symbolic and numeric eigensolvers, here is an overview Aug 11, 2016 at 17:43

I've encapsulated the code of the mysterious user21 into a helmholzSolve command. The code is at the end of this post. It adds very little to user21's code but it does allow us to examine multiple examples quite easily, though it has certainly not been tested extensively and could be improved quite a lot I'm sure. It should be called as follows:

{ev,if,mesh} = helmholzSolve[g_Graphics, n_Integer, opts:OptionsPattern[]];


In this code, g can be a Graphics object, an ImplicitRegion, or a ParametricRegion defining the region in question, n is an integer determining the number of eigenvalues that will be computed, and opts is a list of options to be passed to the discretization functions. It returns ev a list of the computed eigenvalues, if a list of corresponding eigenfunctions represented as InterpolatingFunctions and the mesh for plotting purposes. Using this, we can compute the eigenfunctions of the unit disk is as easy as follows:

{ev, if, mesh} = helmholzSolve[Disk[], 6];
ev
(* Out: {6.80538, 15.7385, 15.7385, 27.477, 27.477, 31.5901} *)


We can visualize the eigenfunctions as follows:

GraphicsGrid[Partition[Table[ContourPlot[if[[k]][x, y], Element[{x, y}, mesh],
PlotRange -> All, PlotPoints -> 50], {k, 1, 6}], 3]] Here's a semi-interesting region:

n = 20;
vertices = Table[(1 + (-1)^k/5) {Cos[2 Pi*k/n], Sin[2 Pi*k/n]}, {k, 1, n}];
g = Graphics[{EdgeForm[Black], Gray, Polygon[vertices]}] And the plot of an eigenfunction:

{ev, if, mesh} = helmholzSolve[g, 6, "MaxCellMeasure" -> 0.005];
Plot3D[-if[][x, y], Element[{x, y}, mesh],
PlotRange -> All, PlotPoints -> 20, Mesh -> All,
MeshStyle -> Opacity[0.3]] Here's an implicitly defined region with a hole:

{ev, if, mesh} = helmholzSolve[
ImplicitRegion[1/4 < x^2 + y^2 && x^4 + y^6 <= 1, {x, y}],
4];
ContourPlot[if[][x, y], Element[{x, y}, mesh],
PlotRange -> All, PlotPoints -> 40] Finally, here's the definition of helmholzSolve.

Needs["NDSolveFEM"];
helmholzSolve[g_, numEigenToCompute_Integer,
opts : OptionsPattern[]] := Module[
{u, x, y, t, pde, dirichletCondition, mesh, boundaryMesh,
nr, state, femdata, initBCs, methodData, initCoeffs, vd, sd,
discretePDE, discreteBCs, load, stiffness, damping, pos, nDiri,
numEigen, res, eigenValues, eigenVectors, evIF},

(* Discretize the region *)
mesh = ToElementMesh[DiscretizeRegion[g], opts],
mesh = ToElementMesh[DiscretizeGraphics[g], opts]
];
boundaryMesh = ToBoundaryMesh[mesh];

(* Set up the PDE and boundary condition *)
pde = D[u[t,x,y], t] - Laplacian[u[t,x,y], {x, y}] +  u[t,x,y] == 0;
dirichletCondition = DirichletCondition[u[t,x,y] == 0, True];

(* Pre-process the equations to obtain the FiniteElementData in StateData *)
nr = ToNumericalRegion[mesh];
{state} = NDSolveProcessEquations[{pde, dirichletCondition,
u[0, x, y] == 0}, u, {t, 0, 1}, Element[{x, y}, nr]];
femdata = state["FiniteElementData"];
initBCs = femdata["BoundaryConditionData"];
methodData = femdata["FEMMethodData"];
initCoeffs = femdata["PDECoefficientData"];

(* Set up the solution *)
vd = methodData["VariableData"];
sd = NDSolveSolutionData[{"Space" -> nr, "Time" -> 0.}];

(* Discretize the PDE and boundary conditions *)
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

(* Extract the relevant matrices and deploy the boundary conditions *)
stiffness = discretePDE["StiffnessMatrix"];
damping = discretePDE["DampingMatrix"];

(* Set the number of eigenvalues ignoring the Dirichlet positions *)
pos = discreteBCs["DirichletMatrix"]["NonzeroPositions"][[All, 2]];
nDiri = Length[pos];
numEigen = numEigenToCompute + nDiri;

(* Solve the eigensystem *)
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 ;

(* Return the relevant information *)
{eigenValues, evIF, mesh}
]

• It is also possible to solve unconstrained eigenvalue problems: In that case no DirichletCondition are applied. Aug 1, 2014 at 7:10
• Thank you very much. Thats exactly what I needed.Just a minor thing: Adding $Head[g] === RegionIntersection$ and $Head[g] === RegionUnion$ to the If condition allows for an easy creation of the boundary. Just as a comment: For large k values the FEM method is known to have problems due to the so called pollution effect. One has to carefully select the mesh if quantitative correct results are desired. Aug 1, 2014 at 9:37

OK, there is good news and there is bad news. In the current version 10 there is no way to do this directly. That's the bad news. The good news is that finite element framework used within NDSolve is exposed and documented; for maximum "hackability" convenience.

Needs["NDSolveFEM"]

(*Compute an eigenfunction of a Helmholzian on the Koch \
Snowflake*)(*lev controls the level of the approximation.*)
lev = 3;
boundaryLev = lev + 1;
KochStep[{p1_, p2_}] :=
With[{q1 = p1 + (p2 - p1)/3,
q2 = (p1 + (p2 - p1)/3) + RotationMatrix[-Pi/3].(p2 - p1)/3,
q3 = p1 + 2 (p2 - p1)/3}, {p1, q1, q2, q3, p2}];
KochStep[pp : {{_, _} ..}] :=
Join[Partition[Flatten[Most /@ (KochStep /@ Partition[pp, 2, 1])],
2], {pp[[-1]]}];
KochVertices =
Nest[KochStep,
N@{{3 Sqrt/4, 3/4}, {-3 Sqrt/4,
3/4}, {0, -3/2}, {3 Sqrt/4, 3/4}}, boundaryLev];
Graphics[Point[KochVertices]] We then first generate a boundary mesh and then the full mesh:

boundaryMesh =
ToBoundaryMesh["Coordinates" -> KochVertices,
"BoundaryElements" -> {LineElement[
Partition[Range[Length[KochVertices]], 2, 1, {1, 1}]]}];

boundaryMesh["Wireframe"] mesh = ToElementMesh[boundaryMesh, "MeshOrder" -> 1, "MaxCellMeasure" -> 0.005];
mesh["Wireframe"] Now for the real stuff. An eigen value problem can be considered as a transient problem.

k = 1/10;
pde = D[u[t, x, y], t] - Laplacian[u[t, x, y], {x, y}] + k^2 u[t, x, y] == 0;
Γ = DirichletCondition[u[t, x, y] == 0, True];


We use NDSolve as a pre-processor:

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


Extract the finite element data:

femdata = state["FiniteElementData"]

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


Set up the solution and variable data:

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


Discretize the PDE and the boundary conditions:

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


Extract the system matrices:

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


Deploy the boundary conditions:

DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs]


Set the number of X smallest eigen values we would like to compute but ignore the Dirichlet positions.

nDiri = First[Dimensions[discreteBCs["DirichletMatrix"]]];

numEigenToCompute = 5;
numEigen = numEigenToCompute + nDiri;


Solve the eigen system:

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


Creating the interpolating functions:

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


Plotting the 5th eigen mode

Plot3D[Evaluate[evIF[][x, y]], {x, y} ∈ mesh,
PlotRange -> All, Axes -> None, ViewPoint -> {0, -2, 1.5},
Boxed -> False, BoxRatios -> {1, 1, 1/4}, ImageSize -> 612,
Mesh -> All] Voilà.

If you are interested in eigenmode analysis for structural mechanics, have a look here

• I'm ridiculously excited! Jul 31, 2014 at 16:02
• Excellent. That is all. Jul 31, 2014 at 16:06
• If you set numEigenToCompute=6 and then plot -efIF[], you essentially get the MATLink Logo, though that was computed using a finite difference technique applied to a regular triangular grid. Jul 31, 2014 at 16:22
• Hopefully, a future version will include some functionality do compute this directly and make this monkey business irrelevant. But with the low level FEM functions one can do a lot of stuff regardless if the functionality is already implemented or not. Jul 31, 2014 at 16:54
• Extraordinary work. Beautiful. Aug 1, 2014 at 10:04

# As of v10.3, this is now possible

using the built-in commands DEigensystem and NDEigensystem functions. These work pretty much as you hope they would, along the lines of

DEigensystem[ℒ[u[x, y, …]], u, {x, y, …} ∈ Ω, n] gives the n smallest magnitude eigenvalues and eigenfunctions for the linear differential operator ℒ over the region Ω.

Hat-tip to Mark McClure for pointing these out to me.

As a simple illustration,

{eigenvalues[region], eigenfunctions[region]}=NDEigensystem[
{-Laplacian[u[x,y], {x,y}],
DirichletCondition[u[x,y]==0, True]}
, u[x,y]
, {x,y} ∈ region
, 6
];
Grid[Partition[Table[
Show[{
ContourPlot[
eigenfunctions[region][[j]]
, {x,y} ∈ region
, Frame -> None, PlotPoints -> 60, PlotRange -> Full
, PlotLabel -> eigenvalues[region][[j]]
],
RegionPlot[region, PlotStyle -> None, BoundaryStyle -> {Black, Thick}]
}]
,{j, 1, Length[eigenvalues[region]]}], 3]]


produces with, for example, a disk,

region = Disk[]


plots like or with the star from above,

With[{n=20},
region=Polygon[Table[(1+(-1)^k/5) {Cos[2 Pi*k/n],Sin[2 Pi*k/n]},{k,1,n}]]
]; or for the implicit region with a hole,

region = ImplicitRegion[1/4 < x^2 + y^2 && x^4 + y^6 <= 1, {x, y}]; and even for the Koch snowflake of user21, defined as

lev = 3;
boundaryLev = lev + 1;
KochStep[{p1_, p2_}] :=
With[{q1 = p1 + (p2 - p1)/3,
q2 = (p1 + (p2 - p1)/3) + RotationMatrix[-Pi/3].(p2 - p1)/3,
q3 = p1 + 2 (p2 - p1)/3}, {p1, q1, q2, q3, p2}];
KochStep[pp : {{_, _} ..}] :=
Join[Partition[Flatten[Most /@ (KochStep /@ Partition[pp, 2, 1])],
2], {pp[[-1]]}];
KochVertices =
Nest[KochStep,
N@{{3 Sqrt/4, 3/4}, {-3 Sqrt/4,
3/4}, {0, -3/2}, {3 Sqrt/4, 3/4}}, boundaryLev];
region=Polygon[KochVertices];


the first nine eigenfunctions look like or zooming in to the fourth one, using

Plot3D[
Quiet[eigenfunctions[region][]]
, {x, y} ∈ region
, PlotRange -> Full
, ImageSize -> 600
, Boxed -> False, Axes -> False
, BoxRatios -> {1, 1, 0.25}
, MeshFunctions -> {#3 &}
, PlotPoints -> 100
]


you get This is in no way meant to underappreciate the really good work by Mark and user21, of course - just to say that you'd hope for this functionality to be built-in and, since recently, it finally is.

• Do you know what methods are used in the solving of the eigenvalue problem? Are they mesh based methods? Aug 27, 2016 at 19:54
• @Beni I don't really, and indeed the usual WRI position is rather closed to giving too many details on the implementation. That said, NDEigensystem is almost certainly riding on top of the same Finite Element API that user21's code uses, which can be seen from the documentation of its Method option. (DEigensystem, on the other hand, is quite likely a very different beast, and like most of Mathematica's symbolic integration is probably very different to human methods of solution.) Aug 27, 2016 at 20:08
• @BeniBogosel, I only saw this comment now. NDEigensystem basically follows what was described in the post above. So, yes, it's mesh based and used the finite element method for discretization. Feb 2, 2017 at 13:31