I have recently been plotting eigenfunctions of the laplacian over the unit square using the NDEigensystem command. However, I have noticed something in the plots which puzzles me.

Below is an image of the first 4 eigenfunction plots. However, the second and third eigenfunction do not resemble the standard known eigenfunctions of the unit square: $u_{21}$ and $u_{12}$, and this has me rather puzzled.

The first 4 Eigenfunctions

With this in mind, I decided to investigate what's going on. I first noticed that the second and third eigenvalues, $\lambda_{2}$ and $\lambda_{3}$, are identical (we know this from the separation of variables solution). This seems to be causing the weird patterns in the eigenfunctions (I checked this for higher order eigenvalues/eigenfunctions too and this appears to be the case); wherever the multiplicity of the eigenvalue is greater than 1, the eigenfunctions being plotted are not the standard ones.

On further inspection I noticed that the second eigenfunction is very roughly $\sin \left (2\pi x \right )\sin \left (\pi y \right ) + \frac{10}{13} \sin \left (\pi x \right ) \sin \left (2 \pi y \right )=u_{21}+ \frac{10}{13} u_{21}$ (see below). This is a linear combination of the two symmetric eigenfunctions $u_{12}$ and $u_{21}$.

Plot of $u_{21}+ \frac{10}{13} u_{21}$

Now this makes complete sense; since $u_{21}$ and $u_{12}$ share the same eigenvalue, then the most general eigenfunction is just a linear combination of the two. However, this has puzzled me somewhat and so I wish to ask a few questions:

Q1) How do I prevent mathematica from combining these eigenfunctions? The reason I ask is that I wish to study the eigenfunction plots over unconvential domains, and I won't be able to tell which plots are combinations of eigenfunctions and which are not for more complicated domains.

Q2) How does mathematica 'know' that the eigenvalues are the same? Of course we know this from separation of varaibles but how does mathematica know that they are the same when it's performing a numerical approximation? Surely the values would be ever so slightly different and so mathematica wouldn't even think to combine the eigenfunctions in the first place.

I know this is a strange question to ask but it has got me rather baffled. Ideally, I'm wanting to adapt the code in mathematica so it produces non-combined eigenfunctions for arbitrary domains.

The (very short) code I have been using to get the plots is:

{vals, funs} = 
  NDEigensystem[{-Laplacian[u[x, y], {x, y}], 
    DirichletCondition[u[x, y] == 0, True]}, 
   u[x, y], {x, y} \[Element] 
    Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}], 4];

  funs[[i]], {x, y} \[Element] 
   Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}], 
  AspectRatio -> Automatic, PlotRange -> All, PlotLabel -> vals[[i]], 
  PlotTheme -> "Minimal"], {i, Length[vals]}]
  • 1
    $\begingroup$ The reason is the same as what I said in my comment to your earlier Q: degeneracies can lead to eigenvectors that are arbitrary orthogonal superpositions of the symmetrized states you expected. If your new Q is about achieving this symmetry given the solutions you plotted, then there are several ways to answer it. Is that what you're asking? $\endgroup$
    – Jens
    Commented Jan 30, 2016 at 20:41
  • 1
    $\begingroup$ In more general domains, you won't even see this problem because degeneracies will be very rare. So maybe we aren't really talking about the square here... what shapes are you really interested in? $\endgroup$
    – Jens
    Commented Jan 30, 2016 at 20:47
  • $\begingroup$ Yes - I wish to plot the standard symmetric states. And you are correct, I'm getting to grips with the square before I begin to look at other domains. I intend to look at a variety of different domains. Initially, I wish to look at polygons (I want to look at the pentagon), and then I wish to look at more complicated 2-D shapes. The question I wish to be answered is whether there is a way of telling mathematica to first spot these degeneracies and secondly, prevent the superposition from happening. $\endgroup$
    – Mr S 100
    Commented Jan 30, 2016 at 20:55
  • $\begingroup$ It's easy to spot the superposition in the case of the square, but I'm not so confident that I would be able to spot the degenerate modes of an 'L' shape for example. $\endgroup$
    – Mr S 100
    Commented Jan 30, 2016 at 20:59
  • $\begingroup$ The L shape is very different from the square because the latter is invariant under the symmetry group $C_{4v}$ which causes the degeneracies. There is no such symmetry in the L. So the goal of your question isn't clear to me. You're asking about a highly symmetric example, where the degeneracies can be removed by reducing the domain. But that's not a Mathematica problem, and it's also unrelated to what you'd see in the L shape. All I can say then is: degeneracies are something you see in the spectrum and not in the wave functions. $\endgroup$
    – Jens
    Commented Jan 30, 2016 at 21:56

1 Answer 1


Although the question singles out the square, it is made clear that the actual applications includes other polygonal shapes as well. This means that it's impossible to give a general answer based on the assumption of separability. The square is separable in Cartesian coordinates, but the pentagon (e.g.) is not.

This is why I'm focusing this answer on the exploitation of point symmetries (reflections, rotations) in more general terms. This means that the results here for the special case of the square will not necessarily be in the product form that one obtains from separation of variables, but instead will conform to the classification of symmetries according to irreducible representations of the relevant point group. That's what you would do in the general case when symmetries are present but separation of variables doesn't work. And that's of course the only case in which numerical solutions using NDEigensystem are really needed in the first place.

For the square, you can exploit the $C_{4v}$ symmetry to convert any wave solution to a symmetrized form. Without going into the mathematical details, here is the prescription. One would have to go more deeply into group theory to explain the methodology.

First define the elements of the symmetry group that leaves the square invariant (I use script letters below, but they don't all show up properly unless you paste this into a notebook):

\[ScriptCapitalG] = {ℰ, Subscript[\[ScriptCapitalC], 
   2], Subscript[\[ScriptCapitalC], 4], 
   Subscript[\[ScriptCapitalC], -4], Subscript[σ, vx], 
   Subscript[σ, vy], Subscript[σ, d1], 
   Subscript[σ, d2]};

\[ScriptCapitalG]ℳ = {{{1, 0}, {0, 1}}, {{-1, 
     0}, {0, -1}}, {{0, -1}, {1, 0}}, {{0, 1}, {-1, 0}}, {{1, 
     0}, {0, -1}}, {{-1, 0}, {0, 1}}, {{0, 1}, {1, 0}}, {{0, -1}, {-1,

\[ScriptCapitalG]ℳℐ = 
  Inverse /@ \[ScriptCapitalG]ℳ;

  TraditionalForm /@ \[ScriptCapitalG]ℳ}]


Here, the names of the elements are listed above the 2D matrices that describe the point operations of rotations and reflections in the $xy$ plane, centered in the middle of the square.

For this group, there is a two-dimensional irreducible representation which causes the degeneracies in the spectrum. The second and third solutions in the question lie in the invariant subspace of that representation. What the question asks for is equivalent to projecting these (asymmetric-looking) functions onto the symmetrized basis functions of that subspace. To do this, one can use the following projection operator, where i = 1,2 is the index of the basis function:

projector[i_] = 
   1/4 Sum[\[ScriptCapitalG]ℳ[[n, i, 1]] ψ /. 
         y} -> ({1, 1}/
           2 + \[ScriptCapitalG]ℳℐ[[n]\
].({x, y} - {1, 1}/2))], {n, 

Now test this for the numerical solutions of the question:

{vals, funs} = 
  NDEigensystem[{-Laplacian[u[x, y], {x, y}], 
    DirichletCondition[u[x, y] == 0, True]}, 
   u[x, y], {x, y} ∈ 
    Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}], 10];

The first projection is

f1 = projector[1][funs[[2]]];

ContourPlot[f1, {x, y} ∈ 
  Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}], AspectRatio -> Automatic,
  PlotRange -> All, PlotLabel -> vals[[2]]]


and the second projection is

f2 = projector[2][funs[[2]]];

ContourPlot[f2, {x, y} \[Element] 
  Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}], AspectRatio -> Automatic,
  PlotRange -> All, PlotLabel -> vals[[2]]]


Note that I got the two symmetrized functions by applying two different projections onto the same eigenstate, funs[[2]]. That's because each of the un-symmetrized states had components along both symmetrized basis states.

Expanded method

The above version of the projector is only one among several others. All of them are needed in order to classify the spectrum fully. So I'll define a more general set of projectors for all the irreducible representations of the group $C_{4v}$ now. This first requires entering the character table for the group, which I call dMatrix. It's stored as an Association where the keys are the conventional names of the irreducible representations (irreps). The values are "matrices" with dimension 1 (for the first four) and 2 (for the last irrep, called Ee).

Clear[applyG, dMatrix, A1, A2, B1, B2, Εe, x, y]

irreps = {A1, A2, B1, B2, Ee};

dMatrix = 
  Association @@ 
   MapThread[(#1 -> #2) &, {{A1, A2, B1, B2, 
      Ee}, {{{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}}, \
{{{1}}, {{1}}, {{1}}, {{1}}, {{-1}}, {{-1}}, {{-1}}, {{-1}}}, {{{1}}, \
{{1}}, {{-1}}, {{-1}}, {{1}}, {{1}}, {{-1}}, {{-1}}}, {{{1}}, {{1}}, \
{{-1}}, {{-1}}, {{-1}}, {{-1}}, {{1}}, {{1}}},

applyG[n_][f_] := (f /. 
       y} -> ({1, 1}/
         2 + \[ScriptCapitalG]ℳℐ[[n]].\
({x, y} - {1, 1}/2))]);

projector[irrep_, row_][f_] := 
   Length[\[ScriptCapitalG]] Sum[
    Conjugate[dMatrix[irrep][[i]][[row, 1]]] applyG[i][f], {i, 

To define the new projectors, I use the same prescription as above, but I factored out the application of each group element in a separate function applyG so that it can be adapted more easily to other simulation domains, if needed. The function projector now also depends on an additional label, irrep.

Now you can apply all these projectors to any state whose symmetries you want to analyze. For example, the solution in func[[5]] initially looks like this:

 funs[[5]], {x, y} ∈ 
  Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}], AspectRatio -> Automatic,
  PlotRange -> All, PlotLabel -> vals[[5]]]


But it is actually a superposition of two different degenerate states, belonging to different one-dimensional irreducible representations: A1 and B1. To do the analysis leading to this result, here is a complete table of all the states obtained from NDEigensystem above:

wavePattern[f_] := 
 Module[{d = 
    Threshold[Table[f, {x, 0, 1, .02}, {y, 0, 1, .02}], .1]},
  If[Max[Abs[d]] < .1, Graphics[{}], 
    ColorFunction -> 
      Blend[{White, Darker@Orange, White}, 2 ArcTan[5 x]/Pi + .5]], 
    ColorFunctionScaling -> False, InterpolationOrder -> 0, 
    Frame -> False, Axes -> False]

frames = Table[GraphicsGrid[Partition[Flatten@Table[Table[
        Show[wavePattern[Evaluate[projector[irr, row][funs[[m]]]]], 
         PlotLabel -> Row[{irr, " row", row}]], {row, 
       {irr, irreps}], 3], 
    PlotLabel -> 
     Framed[Style[Row[{"Eigenvalue \[LongEqual] ", vals[[m]]}], 
       Darker[Red]]], Frame -> All], {m, Length[vals]}];



Each frame refers to one of the wave solutions. In each frame, the table shows all the rows of all the irreducible representations. In the corresponding cells, I plot the wave function obtained by the corresponding projection operator, applied to the eigenfunction whose eigenvalue is given at the top. When a cell is empty, the wave solution has no component under that projection. The decomposition for funs[[5]] is the frame with energy 78.111, and you can see that the cells for representations A1, B1 are filled. By looking at the character table in dMatrix, you can reconstruct how the different operations of the group will change the sign of the symmetrized components.

When more than one cell of the table is filled, it means that the given state is a superposition of the functions shown in those cells. I.e., if you were to simply add the wave functions of the two cells for a state with energy 78.111, you would get back the original state (e.g., funs[[5]]).

What this analysis shows is that whereas the degeneracies of eigenvalues 49.35, 128.338 and 167.855 are caused by the symmetry group of the square (its two-dimensional representation Ee), the other degeneracies are "accidental" in that they would not survive any perturbation that deforms the square to a non-separable shape while preserving the invariance under all symmetry operations. This could e.g. be done by truncating all the corners at 45 degree angles, or rounding them.

Accidental degeneracies are called "accidental" because they correspond to degenerate functions that belong to two different irreducible representations. Functions belonging to the same irrep transform in the same way under symmetry operations (e.g., they may or may not change sign under each of the possible reflections at the four symmetry lines of the square). When two such functions with different transformation behavior under symmetries still have the same eigenvalue, then that degeneracy is not caused by the known symmetry group -- this makes them accidental. The "non-accidental" degeneracies are immune to symmetry-preserving perturbations, the accidental ones aren't. However, it sometimes happens that accidental degeneracies aren't really accidental, but due to hidden symmetries that we just haven't considered. Of course in the square, what I call accidental degeneracies have perfectly simple explanations when looking at the problem with separation of variables. But as I said in the introduction, my remarks concern features that don't assume integrability.

This set of functions (projector etc.) can now in principle also be applied to other groups, by changing the definition in the first two lines, i.e., the group elements and their matrices. You may also have to change applyG if the center of symmetry is not the point $(\frac12,\frac12)$.

Finally, it should be noted that this answer assumed the wave solutions to be given as stated in the question. A cleverer use of symmetry would of course be to reduce the computational domain to a wedge shaped subset on which either Dirichlet or Neumann boundary conditions can be imposed. That will automatically lead to the properly symmetrized eigenfunctions. What I did here is "post-processing," which can be useful because the brute-force calculation for the full domain is easier to program.

  • $\begingroup$ This is excellent - thank you! I've just tested it on the 2nd and 3rd eigenfunction and it works. I'll have to invest the time into finding out exactly why this works (and the group theory behind it). However, I set the eigensystem to find the first 15 functions, and tested the projection on the 5th and 6th modes (which both have eigenvalue ~98.69) and the projection isn't working for me here. Am I correct in adjusting it to f5 = projector[1][funs[[5]]]? Will this also work in the case of eigenvalue multiplicity greater than 2; i.e the 31st, 32nd, 33rd mode (all share eigenvalue 493.48)? $\endgroup$
    – Mr S 100
    Commented Jan 30, 2016 at 23:01
  • $\begingroup$ The projector in the case funs[[5]] returns numerically zero, indicating that these states don't belong to the 2D irreducible representation onto which I project. But there are other (1D) representations, each with their own projection operators. You would have to construct them using the rules of group theory and apply them to each eigenstate to isolate the symmetrized components. You'll need a character table for the group and repeat what I did above - it can be automated, but I'll have to leave that to you or someone else for now. $\endgroup$
    – Jens
    Commented Jan 30, 2016 at 23:56
  • $\begingroup$ Amazing! Thank you for all your help – your knowledge has been invaluable. As I know nothing about group theory, it's probably best that I first read up on the mathematics behind this post in order to understand your code. From what I can see, it works – and works well. I do have one final question: you say that the analysis shows that the degeneracies at eigenvalue 98.71 are "accidental". This is obvious as the frame in the table returns irreducible representations which are not of the standard form. Is this something that just need to be spotted on sight? Thank you so much for your help. $\endgroup$
    – Mr S 100
    Commented Jan 31, 2016 at 12:43
  • $\begingroup$ I added another paragraph on accidental degeneracies. $\endgroup$
    – Jens
    Commented Jan 31, 2016 at 19:05
  • $\begingroup$ Excellent! Thanks for clearing that up for me. I would just like to say thanks for taking the time to answer my question and having patience with me and my ignorance. This is honestly a superb answer. I now look forward to learning about the theory behind all this! $\endgroup$
    – Mr S 100
    Commented Jan 31, 2016 at 20:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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