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 read here, the discussion about how to solve one dimensional eigenvalue problem. I am wondering, how can one generalize these methods to two dimensions. For example:


or any other potentials.

What about three dimensions? Is there any general procedure?

share|improve this question
I am curious about, what can one obtain by coupling two harmonic oscillators, such a potential: x^2 + y^2 +xy? Is this solvable? I can find the ground state energy analytically, but I want to see the whole eigenstates and the corresponding eigenvalues. –  NoOne Jun 14 '14 at 4:08
Yes, just rotate the coordinate system and the oscillators decouple (normal mode decomposition, works for any potential of degree 2). –  Jens Jun 14 '14 at 4:09
Instead of de-coupling them, is it possible to solve it directly by Mathematica? –  NoOne Jun 14 '14 at 4:13
Yes, I have a version of my answer in the question you linked that applies to two dimensions. But it was written by a student in my class a few years ago, so I can't post it online. Anyway, this means the method works in 2D, only unfortunately quite slow because obviously the matrices are harder to construct and get much bigger. –  Jens Jun 14 '14 at 4:26
Oh Lord! Bad luck! –  NoOne Jun 14 '14 at 4:29

1 Answer 1

up vote 14 down vote accepted

In order to give one possible answer, I'll just take the isotropic harmonic oscillator in 2D and do a finite-difference calculation by discretizing the xy plane with constant spacing a.

Here is the construction of the resulting matrix for the Hamiltonian, h. I assume the origin of our spatial grid (where the potential minimum is) lies at {0,0}, and the number of grid points in all directions from the origin is nX. Each grid step corresponds to a length of a = 0.2. On this grid, I then define the potential energy as a function v[x_, y_], and evaluate it at the grid points to create the 2D matrix vGrid. But this is the easy part.

Now we have to make a matrix out of the Hamiltonian in the space of tuples of {x,y} positions. This means the rows and columns of the matrix that we want to diagonalize are labeled by the entries of a list of tuples, which I call xyList. The Hamiltonian has one term that corresponds to the potential energy, and it is created in the above basis of tuples by taking the entries of the matrix vGrid corresponding to the tuple labeling each diagonal element of the Hamiltonian. This happens in the DiagonalMatrix term.

Finally, we have to construct the kinetic energy, i.e., the Laplacian. This is the only off-diagonal part of the matrix. How far off the diagonal its matrix elements extend depends on the order of the finite-difference approximation by which we replace the second derivatives in the Laplacian. Here I chose the simplest possible approximation.

nX = 20;
a = .2;

v[x_, y_] := 1/2 (x^2 + y^2)

vGrid = Table[v[a i, a j], {i, -nX, nX}, {j, -nX, nX}];

xyList = Tuples[Range[1, 2 nX + 1], 2];

 h = DiagonalMatrix[
     SparseArray[(vGrid[[##]] & @@@ xyList) + 2/a^2]] - 
    1/(2 a^2) AdjacencyMatrix[GridGraph[{2 nX + 1, 2 nX + 1}]];]

(* ==> {0.009491, Null} *)

v = Eigenvectors[h, -10];

ListDensityPlot[Partition[v[[10]], 2 nX + 1], 
 PlotRange -> All]

ground state

This is the ground state. The eigenvalues have near-degeneracies:

e = Eigenvalues[h, -10]

==> {3.97742, 3.97742, 3.96769, 3.96769, 2.98745, 2.98244, \
2.98244, 1.99247, 1.99247, 0.997494}

The exact values of the energy should be 4, 3, 2, 1 in descending order.

The degeneracies are the reason why for most states you get superpositions that don't look like the rotationally symmetrized Gauss-Laguerre functions of Gauss-Hermite functions - but they are valid solutions nonetheless (to within the accuracy of the discretization):


state 4

In the construction of h, I used a trick to compute the matrix corresponding to the Laplacian: it invokes the AdjacencyMatrix of a GridGraph. In principle, it may be better to instead use NDSolve'FiniteDifferenceDerivative for this. But I think it's a cute fact that this particular AdjacencyMatrix also works. I may get back to adding an alternative approach using the NDSolve functionality later.

This is a quick-and dirty implementation without bells and whistles, but it shows the main steps. It also shows that the boundary effects at the edge of the grid are not important as long as you only look at states that are low enough in energy so as to not reach that boundary. This is also why time-dependent wave packet solutions in the same potential can be modeled on a finite grid.


The next step is to apply the method to an anisotropic harmonic oscillator, with potential energy $V(x,y) = 1/2 (x^2 + y^2 + x y)$, as was asked for in the comments:

nX = 20;
a = .2;

v[x_, y_] := 1/2 (x^2 + y^2 + x y)

vGrid = Table[v[a i, a j], {i, -nX, nX}, {j, -nX, nX}];

xyList = Tuples[Range[1, 2 nX + 1], 2];

Timing[h = 
   DiagonalMatrix[SparseArray[(vGrid[[##]] & @@@ xyList) + 2/a^2]] - 
    1/(2 a^2) AdjacencyMatrix[GridGraph[{2 nX + 1, 2 nX + 1}]];]

(* ==>{0.009491,Null}*)

v = Eigenvectors[h, -10];

ListDensityPlot[Partition[v[[10]], 2 nX + 1], PlotRange -> All]

(* ==> {0.008648, Null} *)

ground state2

As expected, the ground state now is deformed (elliptical) and oriented along a diagonal.

Edit 2: an even more basic approach

To avoid any advanced functions that may not be present in older versions of Mathematica, here is a practically identical formulation where the finite-difference approximation is constructed in a more straight-forward (i.e. low-level) way. I'm again using the anisotropic oscillator of the last example - everything up to xyList is copied verbatim from above.

nX = 20;
a = .2;

v[x_, y_] := 1/2 (x^2 + y^2 + x y)

vGrid = Table[v[a i, a j], {i, -nX, nX}, {j, -nX, nX}];

xyList = Tuples[Range[1, 2 nX + 1], 2];

Evaluate[hash /@ xyList] = Range[Length[xyList]];
hash[{0, n_}] := hash[{2 nX + 1, n}];
hash[{n_, 0}] := hash[{n, 2 nX + 1}];
hash[{2 nX + 2, n_}] := hash[{1, n}];
hash[{n_, 2 nX + 2}] := hash[{n, 1}];

s1 = 
   Flatten[Table[{{hash[{i, j}], hash[{i, j + 1}]} -> 
       1, {hash[{i, j}], hash[{i + 1, j}]} -> 
       1, {hash[{i, j}], hash[{i, j - 1}]} -> 
       1, {hash[{i, j}], hash[{i - 1, j}]} -> 1}, {i, 1, 
      2 nX + 1}, {j, 1, 2 nX + 1}]]];

h1 = 
  DiagonalMatrix[SparseArray[(vGrid[[##]] & @@@ xyList) + 2/a^2]] - 
   1/(2 a^2) s1;

This is the equivalent to h in the previous examples. It is the Hamiltonian matrix that will now be diagonalized to get the eigenvalues and eigenvectors:

e1 = Eigenvalues[h1, -10]

==> {4.09354, 3.77227, 3.58063, 3.39729, 3.07122, 2.88156, \
2.36968, 2.1824, 1.66728, 0.963587}

v1 = Eigenvectors[h1, -10];

ListDensityPlot[Partition[v1[[4]], 2 nX + 1], 
 PlotRange -> All]


Here I chose to plot an excited state instead of the ground state because it shows nicely that after breaking the rotational symmetry we don't have degeneracies anymore and therefore get eigenfunctions that look like the expected Gauss-Hermit functions in both orthogonal directions (along the diagonals), as one would also obtain by first finding the normal modes exactly. This finite difference algorithm of course doesn't know or care about the existence of rotational symmetry or separability, so it will also work with binding potentials that are not harmonic oscillators.

Regarding three dimensions

The generalization to higher dimensions is in principle quite straightforward, but you'll have to consult some text on finite-difference methods to construct the second-derivate stencil properly. In those cases, the availability of built-in methods becomes really useful. There is a lot more to read about this in the documentation (search for FiniteDifferenceDerivative, it will in fact tell you how to construct a 2D Laplacian).

Edit to generalize my first approach to 3D

To generalize my earlier AdjacencyMatrix trick to three dimensions, you would need to calculate the Hamiltonian as follows:

h = DiagonalMatrix[SparseArray[(vGrid[[##]] & @@@ xyList) + 3/a^2]] - 
    1/(2 a^2) AdjacencyMatrix[GridGraph[{2 nX + 1, 2 nX + 1, 2 nX + 1}]]

Note the factor 3 instead of 2 in the additive constant on the diagonal. The constant on the diagonal is 1/a^2 times the number of dimensions of the grid.

The difference between using AdjacencyMatrix and the last approach is in the boundary conditions. The AdjacencyMatrix method does not implement periodic boundary conditions, whereas the last approach does. However, these boundary conditions are irrelevant in the present calculation because the bound states are assumed to not feel the boundary anyway.

I like the AdjacencyMatrix approach mainly because it is a more descriptive formulation, saying in code what the topology of the laplacian matrix really looks like.

Edit regarding grid spacing

To change the grid spacing without changing the physics, you have to keep the product a*nX large enough to ensure that the eigenstates you're interested in are not affected by the grid boundary. The value a*nX is the largest x or y coordinate that will be sampled, and your potential at those coordinates should be high enough to exclude the wave function of any relevant eigenstate from the boundary.

One could also combine this algorithm with my answer here to insure that the eigenvalues and eigenvectors are sorted in the correct way.

share|improve this answer
Fantastic. But it is very concise. I am going through it, but it would be great if you could add some extra explanation to it. –  NoOne Jun 14 '14 at 5:21
It seems to me NDSolve will be more straightforward. I am looking enthusiastically for how it works. God bless you! –  NoOne Jun 14 '14 at 5:26
Actually, this is pretty much the simplest approach I can think of. By using NDSolve functionality, I only meant that the approximation for the second derivatives in the kinetic energy can be improved using some built-in functions. But the remainder, i.e., the construction and diagonalization of h, will stay the same. I'll hopefully get to that modification tomorrow. –  Jens Jun 14 '14 at 5:38
I would like to go through a Fourier analysis of its wave packet. Is your approach here for finding the evolution of the wave packet, applicable for 2D H. O.? Because there, the potential you have considered is periodic, but here is not. –  NoOne Jun 14 '14 at 5:52
@NoOne No, you do have to modify the potential to make it consistent with the periodic boundary conditions. This could be done artificially while still maintaining a parabolic minimum in the center. But I'll have to let you try to figure that out yourself... in the meantime, could you modify your original question to include the ho potential $x^2 + y^2 + x y$? That's the more doable test case, after all. –  Jens Jun 14 '14 at 5:56

Your Answer


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.