Numerically solving Helmholtz equation in 2D for arbitrary shapes

Question

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?

Answer 1

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["NDSolve`FEM`"];
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 *)
  If[Head[g] === ImplicitRegion || Head[g] === ParametricRegion,
    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} = NDSolve`ProcessEquations[{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 = NDSolve`SolutionData[{"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 *)
  load = discretePDE["LoadVector"];
  stiffness = discretePDE["StiffnessMatrix"];
  damping = discretePDE["DampingMatrix"];
  DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs];

  (* 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}
]

Answer 2

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.

Let's start with a region that @MarkMcClure would consider interesting.

We load our favorite package:

Needs["NDSolve`FEM`"]

(*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} = NDSolve`ProcessEquations[{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 = NDSolve`SolutionData[{"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

Answer 3

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[3]/4, 3/4}, {-3 Sqrt[3]/4, 
      3/4}, {0, -3/2}, {3 Sqrt[3]/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][[4]]]
 , {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.