# Finding the eigenfunctions of one and two dimensional Harmonic Oscillator

(Edited) For finding the ground state wave function of:

$H\psi(x) = (-1/2)d^2\psi(x)/dx^2 + (1/2)x^2\psi(x) = E \psi(x)$

I have written:

mOneDSchEq[n_] :=
Table[Switch[i - j, -1, p[x[i]],
0, (10/(n + 1))^2 q[x[i]] - 2 p[x[i]], 1, p[x[i]], _, 0], {i,
n}, {j, n}];

q[x_] := -x^2; p[x_] := 1;
Xarray[n_] := Do[x[i] = -5 + i 10/(n + 1), {i, 0, n + 1}];

EigVec[n_] := Eigenvectors[mOneDSchEq[n]];
lisEigVec = EigVec;
OneEigVec[j_] := Part[Reverse[lisEigVec], j];
y[i_] := Part[OneEigVec, i];
listOfPoints =
Join[{{x, 0}}, Table[{x[i], y[i]}, {i, 1, 35}], {{x, 0}}];
ListPlot[listOfPoints, PlotJoined -> True, PlotRange -> All,
PlotLabel -> "Ground State Wave Function of Harmonic Oscillator",
AxesLabel -> {"x", "y"}]


Which I have obtained the Gaussian, correctly.

The question that came to my mind is that:

Is it possible by knowing the ground stat eigenvalue, i.e. 1/2, solving the Schrödinger equation numerically, and obtain the ground state wave function? in other words, to solve:

$H\psi(x) = (-1/2)d^2\psi(x)/dx^2 + (1/2)x^2\psi(x) = (1/2) \psi(x)$

or

$H\psi(x) = (-1/2)d^2\psi(x)/dx^2 + (1/2)x^2\psi(x) = (3/2) \psi(x)$

So, I wrote:

s = NDSolve[{-(1/2) \[Psi]''[x] + (1/2) x^2(\[Psi][x]) == (1/2) \[Psi][
x], \[Psi][-5] == 0, \[Psi] == 0}, \[Psi], {x, -5, 5}]

Plot[Evaluate[\[Psi][x] /. s], {x, -5, 5}, PlotRange -> All]


BUT, I got nothing. What is the problem?

The other question is that, I was traveling through the website and found an elegant approach to two dimensional Harmonic Oscillator here.

My question is, if again we want to solve the Schrödinger equation numercally and obtain wave functions, now two dimensional, by knowing the eigenvalues, what should we do? For example:

$H\psi(x,y) = (-1/2)(d^2/x^2 + d^2/dy^2)\psi(x,y) + (1/2)(x^2 + y^2)\psi(x,y) = (1) \psi(x,y)$

and

$H\psi(x,y) = (-1/2)(d^2/x^2 + d^2/dy^2)\psi(x,y) + (1/2)(x^2 + y^2+ x y)\psi(x,y) = (0.96) \psi(x,y)$

• In Jens' answer, isn't the 1/(2 a^2) bit there to take into account the factor of 1/2 in front of the laplacian? Also, the Partition is there because he is representing 2d space in a 1d vector (basically, he discretises space, then take the 2d matrix and set the rows one after the other to each other so as to form a 1d vector; the Partition undoes this). I'll write an answer tomorrow if I remember, have time and nobody else does! – acl Jun 23 '14 at 0:16
• Actually I removed "2" in the denominator of 1/(2 a^2), but the answer seemed to me non-sense! and also, why did not plot the two dimensional wave function directly, and he used ListDensityPlot? – SomeBody Jun 23 '14 at 0:21
• Hmm sorry, I didn't understand; can you rephrase? – acl Jun 23 '14 at 0:24
• I removed 2 in the denominator and also 1/2 in the potentials, so I expect to obtain two times of energy, i.e. 2E; but It was not. – SomeBody Jun 23 '14 at 0:26
• to Jens: I tried to comment there, but I was not allowed! – SomeBody Jun 23 '14 at 0:28

To give another answer for the one-dimensional harmonic oscillator, let's use a different approach based on the NDSolve functionality I alluded to in the linked answer. Edit: I also update the linked answer to include the analogue of this approach in two dimensions.

n = 2000;
a = .02;
grid = N[a Range[-n, n]];
derivative2 =
NDSolveFiniteDifferenceDerivative[2, grid]["DifferentiationMatrix"]


SparseArray[<20009>,{4001,4001}]

potential = Map[(1/2 #^2) &, grid];

hamiltonian = -derivative2/2 +
DiagonalMatrix[SparseArray[potential]];

eigenvalues = Chop[Eigenvalues[hamiltonian, -10]]


{9.5, 8.5, 7.5, 6.5, 5.5, 4.5, 3.5, 2.5, 1.5, 0.5}

v = Chop[Eigenvectors[hamiltonian, -10]];

ListLinePlot[{Abs[v[[-1]]]^2, Abs[v[[-2]]]^2,
Abs[v[[-3]]]^2}, DataRange -> grid[[{1, -1}]],
PlotRange -> {{-4, 4}, All}] Here I used a grid spacing of a = 0.02 and get numerically very exact solutions for the lowest states of the harmonic oscillator.

The matrix representing the second derivatives (derivative2) in the Laplacian is generated using FiniteDifferenceDerivative.

To address some of the other issues in the question:

The initial code in the question didn't produce a result for me. However, since you state you got the desired result, I assume that there is some typo in the question. Definitely, one can improve the first code block by wrapping the generated Hamiltonian matrix in N to make it into a machine-precision matrix that can be diagonalized much faster.

However, the main question seems to have been: why does the differential equation

s =
NDSolve[{-(1/2) ψ''[x] + (1/2) x^2 (ψ[x]) == (1/2) ψ[
x], ψ[-5] == 0, ψ == 0}, ψ[x], {x, -5, 5}];
Clear[x];
ψSol[x_] = ψ[x] /. s[[1, 1]];

Plot[Evaluate[ψSol[x]], {x, -5, 5},
PlotRange -> All]


yield an apparently empty plot? The answer is that the boundary conditions are incorrect if you're looking for a non-trivial solution. The solver actually finds the only possible answer, $\psi(x)\equiv 0$ for all $x$. But this is because you forced the wave function to be zero at two points whereas the ground state by definition has no nodes!

So you should solve the following equation instead:

s =
NDSolve[{-(1/2) ψ''[x] + (1/2) x^2 (ψ[x]) == (1/2) ψ[
x], ψ == 1, ψ' == 0}, ψ[x], {x, -5, 5}];
Clear[x];
ψSol[x_] = ψ[x] /. s[[1, 1]];

Plot[Evaluate[ψSol[x]], {x, -5, 5},
PlotRange -> All] This yields the expected Gaussian. I chose boundary conditions for the function to be 1 and its derivative to be 0 at the origin.

• It is much more complicated than I expected. But, as I asked in the question, if, for example the eigenvalue is given 1/2, why does not my simple approach in the question work? – SomeBody Jun 23 '14 at 1:21
• @JensI assumed -5 and 5 physical infinity, where the wave function vanishes. – SomeBody Jun 23 '14 at 1:42
• But then you can't blame Mathematica for giving you the exact solution satisfying your approximate assumption. – Jens Jun 23 '14 at 1:44
• And in your two dimensional isotropic case, if E = 1 is given for example, is your plot there, reproducible by solving the differential equation for E = 1? – SomeBody Jun 23 '14 at 1:47
• No, that won't work. The real problem is that this is an elliptical boundary-value problem, for which NDSolve is not always able to give the solutions you want. That's why the approach based on matrix diagonalization is still necessary. In one dimension, the boundary-value problem isn't too hard because there is only one independent variable. But when there are more than 2 variables, you have a partial differential equation for which NDSolve` is usually not able to find a solution. – Jens Jun 23 '14 at 1:50