I have been following the method for numerically solving the Helmholtz equation in this example (the answer by User21) and have come across two problems. I have been implementing the method for a 2x1 rectangle domain. The code I have used is as follows:


boundaryMesh = 
 ToBoundaryMesh[Polygon[{{1, 1/2}, {1, -1/2}, {-1, -1/2}, {-1, 1/2}}]]

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

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

nr = ToNumericalRegion[mesh];
{state} = 
  NDSolve`ProcessEquations[{pde, \[CapitalGamma], u[0, x, y] == 0}, 
   u, {t, 0, 1}, {x, y} \[Element] nr];

femdata = state["FiniteElementData"]

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

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

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

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

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

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

numEigenToCompute = 5;
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;

I then plot the 2nd eigenfunction (for example) using:

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

and I have (attempted) to plot the nodal lines of the eigenfunction using:

ContourPlot[Evaluate[evIF[[2]][x, y]]==0, {x, y} ∈ mesh, 
 PlotRange -> All, Axes -> None, ImageSize -> 612, 
 PlotPoints -> 40]

The first problem is this:

It is well known that the nodal lines of the rectangle have the form $u_{mn}$, but the code seems to only give the eigenfunctions for $u_{m1}$ ($n$ is fixed at $1$) [see the below image]. I would be really grateful if somebody could tell me why this is happening. I am just starting to learn the basics of mathematica, so my knowledge is quite limited. I have been reading over the code but cannot seem to understand why it produces only eigenvalues of this form. In the image below, you can see that the code fails to produce eigenfunctions for $u_{12}$, $u_{21}$ and $u_{13}$.


The second issue is not so much of a problem, but more myself seeking some advice from the experts. Basically, I would like to know whether my contour plot for the nodal lines is the best way to approach this task. I have naively asked Mathematica to plot when the eigenfunction vanishes. Is there a better way of doing this?

I am very new to mathematica and I am not the most experienced of individuals when it comes to mathematics. Hence, if there are obvious flaws in my understanding of mathematica, or the mathematical material, then I would be more than happy for you to point them out.

I look forward to reading your responses.

  • $\begingroup$ The eigenfunctions are sorted by their energy eigenvalue, and $u_{1,2}$ is higher than you probably expected. It is found as evIF[[4]]. So there seems to be no problem. Other eigenfunctions may have degeneracies that will lead to superpositions with nodal lines that aren't straight lines. So I think there is no problem with Mathematica's solutions. $\endgroup$
    – Jens
    Commented Jan 27, 2016 at 18:44
  • $\begingroup$ Thanks for the prompt response Jens. You are right - evIF[4] displays the $u_{12}$ eigenfunction. Please forgive my ignorance, but how would I know which number corresponds to $u_{13}$ for example? Is there a way of knowing this for a general $u_{mn}$? I am a little confused with what you mean by their energy eigenvalues. I know this is probably trivial but this is new material for me. Thanks in advance. $\endgroup$
    – Mr S 100
    Commented Jan 27, 2016 at 19:04
  • $\begingroup$ FEM doesn't use assumptions about separability. But this is what the textbook labeling scheme of the eigenfunction is based on. So the only criterion by which FEM is able to sort the solutions is their eigenvalue. Node numbers based on separability provide a different sorting scheme but not in brute-force approaches such as this. I use "energy" above, but you can omit that term. It's an interpretation only. $\endgroup$
    – Jens
    Commented Jan 27, 2016 at 19:35
  • $\begingroup$ I understand - that makes sense. So essentially what mathematica is doing is sorting the eigenvalues in order of magnitude (is the order smallest to largest?) and that's how it distinguishes the order of the eigenfunctions. Is this correct? Thus, for our ordered eigenvalues $\lambda_{1}<\lambda_{2}<\lambda_{3}< \cdots $ we get eigenfunctions evIF[1], evIF[2], evIF[3],.... So when analysing nonseperable domains, it makes more sense to say the first, second, third eigenfunction, rather than $u_{01}$,$u_{11}$ etc. Am I on the right lines? Thanks Again Jens! $\endgroup$
    – Mr S 100
    Commented Jan 27, 2016 at 20:15
  • $\begingroup$ Yes, that's correct. $\endgroup$
    – Jens
    Commented Jan 27, 2016 at 22:04