# Numerically solving Helmholtz equation in 2D for arbitrary shapes

I would like to solve the Helmholtz equation with dirichlet boundary conditions in 2 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.

$\triangledown^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, unforuntalty 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. –  m0nhawk Jul 31 '14 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. –  episanty Aug 1 '14 at 13:44

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[[6]][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[[4]][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. –  user21 Aug 1 '14 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. –  Julian S. Aug 1 '14 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[3]/4, 3/4}, {-3 Sqrt[3]/4,
3/4}, {0, -3/2}, {3 Sqrt[3]/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[[5]][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! –  Mark McClure Jul 31 '14 at 16:02
Excellent. That is all. –  RunnyKine Jul 31 '14 at 16:06
If you set numEigenToCompute=6 and then plot -efIF[[6]], you essentially get the MATLink Logo, though that was computed using a finite difference technique applied to a regular triangular grid. –  Mark McClure Jul 31 '14 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. –  user21 Jul 31 '14 at 16:54
Extraordinary work. Beautiful. –  duffymo Aug 1 '14 at 10:04