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.

(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[35];
OneEigVec[j_] := Part[Reverse[lisEigVec], j];
y[i_] := Part[OneEigVec[1], i];
listOfPoints = 
Join[{{x[0], 0}}, Table[{x[i], y[i]}, {i, 1, 35}], {{x[36], 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][5] == 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) $

Thanks for your attention!

share|improve this question
    
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 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 at 0:21
    
Hmm sorry, I didn't understand; can you rephrase? –  acl Jun 23 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 at 0:26
    
to Jens: I tried to comment there, but I was not allowed! –  SomeBody Jun 23 at 0:28

1 Answer 1

up vote 2 down vote accepted

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.

n = 2000;
a = .02;
grid = N[a Range[-n, n]];
derivative2 = 
 NDSolve`FiniteDifferenceDerivative[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}]

plot of three functions

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, ψ[5] == 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], ψ[0] == 1, ψ'[0] == 0}, ψ[x], {x, -5, 5}];
Clear[x];
ψSol[x_] = ψ[x] /. s[[1, 1]];

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

gaussian

This yields the expected Gaussian. I chose boundary conditions for the function to be 1 and its derivative to be 0 at the origin.

share|improve this answer
    
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 at 1:21
    
@JensI assumed -5 and 5 physical infinity, where the wave function vanishes. –  SomeBody Jun 23 at 1:42
    
But then you can't blame Mathematica for giving you the exact solution satisfying your approximate assumption. –  Jens Jun 23 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 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 at 1:50

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.