Find eigen energies of time-independent Schrödinger equation
You asked for alternative approaches to what you did, so here is one:
A completely different approach to the one-dimensional time-independent Schrödinger equation would be to use matrix techniques. The idea is to eliminate the need for NDSolve entirely. For bound-state problems, you can do this by choosing a basis satisfying the condition of vanishing wave function at infinity and expanding all quantities in that basis.
A particularly nice way of doing this has been pointed out in:
H. J. Korsch and M Glück, "Computing quantum eigenvalues made easy" Eur. J. Phys. 23 (2002) 413-426.
It uses the harmonic-oscillator basis, but also introduces the relation between position and momentum operators on one hand, and the harmonic-oscillator operators on the other hand. This makes it very useful as a teaching method.
I implemented that approach here, with some refinements to satisfy given tolerances and make plots:
The potential energy $U(x)$ is assumed to be a polynomial of finite degree in a single variable. Since we're using a harmonic-oscillator basis to find the eigenstates of the Hamiltonian, a finite-degree polynomial translates to a set of finite powers of the position matrix xOp. By increasing the dimension of these matrices while keeping the powers fixed to their finite values, we can achieve convergence of the energy eigenvalues and eigenvectors near the potential minimum.
The larger the matrix dimension, the more of the low-lying eigenvalues will converge to the the true eigenvalues.
First I define position and momentum operators in the basis of harmonic-oscillator eigenfunctions, and a matrix representing the kinetic energy. The dimension of the matrices is dim:
Define basic matrices
xOp[dim_] :=
SparseArray[{Band[{1, 2}] -> #, Band[{2, 1}] -> #}, {#, #} &@dim] &[
Table[Sqrt[i], {i, dim - 1}]/Sqrt[2]]
eKin[dim_] :=
1/2 SparseArray[{Band[{1, 1}] -> Table[i - 1/2, {i, dim}],
Band[{1, 3}] -> #, Band[{3, 1}] -> #}, {#, #} &@dim] &[
Table[-Sqrt[i (i + 1)/2], {i, dim - 2}]/Sqrt[2]]
Spectrum and eigenfunctions
This is the heart of the computation.
The function called spectrum determines the eigenvalues and eigenvectors of a Hamiltonian with the above kinetic energy and a polynomial-type potential energy term. The number of eigenvalues to be returned is given as the first argument, n, and their maximum error is given by tolerance (the latter can be left out if the default value is sufficient):
spectrum[n_, tolerance_: .00001][potential_, var_] /;
PolynomialQ[potential, var] := Module[
{
x, h, e, v,
c = CoefficientList[potential, var],
min = First@Check[
Minimize[potential, var],
Print[
"No absolute potential minimum found. Make sure the leading \
power has a positive coefficient!"];
Abort[],
Minimize::natt
]
},
FixedPoint[
{#[[1]] + 1,
Eigenvalues[
h = eKin[#[[1]]] +
N[Sum[c[[i]] MatrixPower[xOp[#[[1]]], i - 1], {i, 2,
Length[c]}] + (c[[1]] - min) IdentityMatrix[#[[
1]]]], -n]} &, {n + 1, 2 Reverse@Range[n] - 1},
SameTest -> (Max@Abs[#1[[2]] - #2[[2]]] < tolerance &)];
{e, v} = Reverse /@ Eigensystem[h, -n];
{e + min, ψe /@ v}
]
The second set of arguments above (potential, var) are the potential energy U(x) and the name of the position variable - usually it would be x.
I had to make a change for compatibility with Mathematica version 10 because it throws an error when you try to apply MatrixPower to a numerically singular matrix with power 0. Position matrices can become numerically singular for large dimension, so I have to treat the zeroth power of xOp separately. This wasn't necessary in version 8. Since that term only determines a constant energy offset, it is lumped together with another constant term that I introduce in order to shift the potential minimum to zero (so that the Eigenvalues are sorted in ascending order, starting with the ground state - this shift by min is undone at the end of the calculation).
Output results in position basis
To get the wave function, we implement the sum over basis states to get the wave functions ψe (here, e doesn't have anything to do with parity - the method works independently of parity symmetry). The output format is defined below to yield a shortened representation that doesn't spit out the entire list of coefficients, only listing their number.
Format[ψe[l_List]] := ψe[
Row[{"\[LeftSkeleton]", Length[l], "\[RightSkeleton]"}]]
ψe[l_List][x_] := 1/(Pi^(1/4))Exp[-x^2/ 2] l.(1/Sqrt[2^# #!] HermiteH[#, x] &[Range[Length[l]] - 1])
plotSpectrum[data_, potential_, {var_, min_, max_}] :=
Show[
Apply[
Plot[
Evaluate[#1 + (Last[#1] - First[#1])/(2 + Length[#1])/
2 Through[#2[var]]], {var, min, max},
Prolog -> {Gray, Map[Line[{{min, #}, {max, #}}] &, #1]},
Frame -> True
] &,
data
],
Plot[potential, {var, min, max}]
]
Examples
Check that we get the correct numerical eigenvalues for the harmonic oscillator potential $U(x) = \frac{1}{2}x^2$ in dimensionless units:
spectrum[10][1/2 x^2, x]
{{0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5}, {ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>]}}
The first list contains the eigenvalues, the second list has the eigenfunctions in shortened form, ready to be plotted. Note from the short forms that that the number of coefficients used in ψe is 15 here, even though I specified that I want to find only 10 eigenvalues in spectrum. This is because these 10 eigenvalues will not be accurate to the given tolerance unless the Hamiltonian matrix is diagonalized in a larger Hilbert space of 15 basis functions. Using FixedPoint, the procedure iteratively increased the number of basis function until none of the first 10 eigenvalues changed by more than the given tolerance in the next iteration.
Example 2: quartic oscillator
Now we'll look at the eigenstates of a quartic potential, shown here:
Plot[-3 x^2 + x^4, {x, -2, 2}]
This is a double-well potential, and the low-lying eigenstates will be localized except for a small tunneling coupling. The potential energy is still a polynomial in $x$, so the method should also work here.
ψ1 = spectrum[5][-10 x^2 + x^4, x] // Last // Last
ψe[<<72>>]
Plot[Evaluate[ψ1[x]], {x, -5, 5}]
If we plot all eigenstates together, the amplitudes are compressed, but you can discern the distinction between even and odd states:
With[{v = -10 x^2 + x^4, n = 10},
plotSpectrum[
spectrum[n][v, x], v, {x, -4, 4}]
]
As you can see in the quartic-oscillator example, the number of basis states (72) is larger for more complicated potentials. But for the harmonic oscillator test case, it's no problem to get very high accuracy for highly excited levels:
First[spectrum[100][1/2 x^2, x]]
{0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 61.5, 62.5, 63.5, 64.5, 65.5, 66.5, 67.5, 68.5, 69.5, 70.5, 71.5, 72.5, 73.5, 74.5, 75.5, 76.5, 77.5, 78.5, 79.5, 80.5, 81.5, 82.5, 83.5, 84.5, 85.5, 86.5, 87.5, 88.5, 89.5, 90.5, 91.5, 92.5, 93.5, 94.5, 95.5, 96.5, 97.5, 98.5, 99.5}
As Jens mentioned, the spatially discretize the equation is another alternative for bound state problem. Here is my very simple implementation of this approach. The basic idea is express the equation on a grid. The differentials can be expressed as finite differences. For example, the second order derivative can be expressed as
$$ \frac{d^2\psi}{d x^2}\approx \frac{\psi(x+\Delta x)+\psi(x-\Delta x)-2\psi(x)}{{\Delta x}^2} $$
So the $\frac{d^2}{dx^2}$ operator can be expressed as a tridiagonal matrix(zero boundary condition)
$$ \left[ \begin{matrix} -2 & 1 & & 0\\ 1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ 0 & & 1 & -2 \end{matrix} \right] $$
and potential operator is an diagonal matrix
$$ \left[ \begin{matrix} V(x_1) & & & 0\\ & \ddots & & \\ & & \ddots & \\ 0 & & & V(x_N) \end{matrix} \right] $$
so the equation in matrix form is
$$ \left[ \begin{matrix} 2a_0+V(x_1) & -a_0 & & 0\\ -a_0 & \ddots & \ddots & \\ & \ddots & \ddots & -a_0 \\ 0 & & -a_0 & 2a_0+V(x_N) \end{matrix} \right] \cdot \left[\begin{matrix} \Psi_1\\ \vdots\\ \vdots\\ \Psi_N \end{matrix}\right] = E \left[\begin{matrix} \Psi_1\\ \vdots\\ \vdots\\ \Psi_N \end{matrix}\right] $$
where $a_0=\frac{\hbar^2}{2m(\Delta x)^2}$.
Diagonalize the hamiltonian will give the eigen-states.
Code
TISE1D[U_Function, {xmin_, xmax_}, N0Grid_: 101,
BoundaryCondition_String: "zero"] := Module[
{Δx = (xmax - xmin)/(N0Grid - 1), Hmtx, Tmtx, Vmtx},
Tmtx = -(1/(2 (Δx)^2))
SparseArray[{{i_, i_} -> -2, {i_, j_} /; Abs[i - j] == 1 ->
1}, {N0Grid, N0Grid}];
Vmtx = DiagonalMatrix[U /@ Range[xmin, xmax, Δx]];
Hmtx = Tmtx + Vmtx;
If[BoundaryCondition == "periodic",
Hmtx[[1, -1]] = Hmtx[[-1, 1]] = -(1/(2 (Δx)^2));
];
Sort[Transpose@Eigensystem[Hmtx], (#1[[1]] < #2[[1]]) &]
]
Examples
1.Harmonic oscillator
V1[x_] = 1/2. x^2;
Transpose[TISE1D[Function[{x}, V1[x]], {-10, 10}, 400]][[1, 1 ;; 4]]
(*{0.499921,1.49961,2.49898,3.49804}*)
ListPlot[TISE1D[Function[{x}, V1[x]], {-10, 10}, 400][[1 ;; 4, 2]],
Joined -> True, PlotRange -> {{-5, 5}, All}, DataRange -> {-10, 10},
Axes -> False, Frame -> True]
2.Reflectionless potential
V2[x_] = With[{\[HBar] = 1, m = 1, a = 1.4}, -((\[HBar]^2 a^2)/m) Sech[a x]^2];
For this potential the analytical bound state wave function is $\psi(x)=A*sech(a x)$ and the ground state energy is $E_0=-\frac{a^2 \hbar ^2}{2 m}$.
grd = TISE1D[Function[{x}, V2[x]], {-10, 10}][[1, 2]];
With[{a = 1.4}, {TISE1D[Function[{x}, V2[x]], {-10, 10}, 200][[1, 1]], -1.4^2/2}]
(*{-0.980757,-0.98}*)
Show[ListPlot[Abs[grd]^2, PlotRange -> All, Joined -> True,
DataRange -> {-10, 10}],
Plot[Max[Abs[grd]^2]*Sech[1.4 x]^2, {x, -10, 10}, PlotRange -> All,
PlotStyle -> Red]]
3.potential $V(x)=-\frac{b \hbar ^2}{m}\text{sech}^2(a x)$.
V3[x_] = With[{\[HBar] = 1, m = 1, a0 = 10., a = 0.81, b = 1.94}, -((\[HBar]^2 b)/m) Sech[a x]^2];
bdE = Select[Transpose[TISE1D[Function[{x}, V3[x]], {-10, 10}, 200]][[1]]*27.211, # < 0 &]
(*{-35.1018, -8.64135}*)
ListPlot[TISE1D[Function[{x}, V2[x]], {-10, 10}, 200][[1 ;; 2, 2]],
PlotRange -> All, Joined -> True, DataRange -> {-10, 10},
Axes -> False, Frame -> True]
If the potential is $(\tanh (x)+1) (\tanh (x)-1)$ you can obtain the analytic solution using Mathematica as follows:
[I have omitted some of the detail - e.g. the asymptotic expansions - because the details are analogous to the simple harmonic oscillator case in my previous answer (see above).]
Define the potential.
u[x_] = (1 + Tanh[x]) (-1 + Tanh[x])
Define the differential equation.
diffeq = -(1/2) \[Psi]''[x] + u[x] \[Psi][x] == e \[Psi][x]
Solve the differential equation.
soln = DSolve[diffeq, \[Psi], x][[1, 1]] // Simplify
(* \[Psi] -> Function[{x},
C[1] LegendreP[1, I Sqrt[2] Sqrt[e], Tanh[x]] +
C[2] LegendreQ[1, I Sqrt[2] Sqrt[e], Tanh[x]]] *)
Verify this solution.
diffeq /. soln // FullSimplify
(* True *)
[This is where I have omitted the asymptotic expansion details.] The LegendreQ piece diverges as $x \longrightarrow \infty$, so use C[2] -> 0 to suppress it. The LegendreP piece diverges as $x \longrightarrow -\infty$ with one exception when e -> -(1/2).
So the bound state solution is
\[Psi][x] /. soln /. {C[2] -> 0, e -> -(1/2)}
(* 1/2 C[1] Sqrt[1 - Tanh[x]^2] *)
Here is a Manipulate that illustrates how the solution of the differential equation depends on the energy.
Manipulate[
Plot[LegendreP[1, I Sqrt[2] Sqrt[e], Tanh[x]] // Chop, {x, -5, 5},
AxesLabel -> {"x", "\[Psi]"}], {{e, -1/2, "e"}, -1, 0}]