Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
Look into What's New in M10 for a start. –  m0nhawk Jul 31 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 at 13:44

2 Answers 2

up vote 31 down vote accepted

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]]

enter image description here

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

enter image description here

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]]

enter image description here

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]

enter image description here

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}
]
share|improve this answer
    
It is also possible to solve unconstrained eigenvalue problems: In that case no DirichletCondition are applied. –  user21 Aug 1 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 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.

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]]

enter image description here

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"]

enter image description here

mesh = ToElementMesh[boundaryMesh, "MeshOrder" -> 1, "MaxCellMeasure" -> 0.005];
mesh["Wireframe"]

enter image description here

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]

enter image description here

Voilà.

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

share|improve this answer
7  
I'm ridiculously excited! –  Mark McClure Jul 31 at 16:02
5  
Excellent. That is all. –  RunnyKine Jul 31 at 16:06
3  
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 at 16:22
    
Why k=1/10? In fact, it seems that you can suppress the k entirely as t is playing the role of the parameter. –  Mark McClure Jul 31 at 16:44
1  
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 at 16:54

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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