Periodic potentials and Bloch waves
This is a completely different approach, using functionality that was already present in Mathematica version 8. I'm posting this because it's quite robust and also allows me to demonstrate how Bloch's theorem can be used here to get a more general class of non-periodic solutions in a periodic potential. This is equivalent to imposing boundary conditions that continuously interpolate between between periodic and anti-periodic boundary conditions (as measured by the phase difference between the function values at the left and right boundary). The original question is a special case of the solutions I obtain here.
The starting point would be a periodic potential $U(x)$ with periodicity interval $\ell$, and look for solutions of the form
$$\psi_k(x)\equiv \exp(i k x) u_k(x)$$
where $u_k(x) = u_k(x+\ell)$ is the periodic part. The additional parameter $k$ is the Bloch wave number, and it gives us periodic solutions $\psi_k$ if $k=0$ (or an integer multiple of $2\pi/\ell$). On the other hand, if $k=\pi/\ell$, the solution is "antiperiodic." Any value $k\in[-\pi/\ell,\pi/\ell]$ leads to a unique solution (that's the "first Brillouin zone").
Inserting this Bloch form into the wave equation, we get an equation for $u_k(x)$ which is always solved with purely periodic boundary conditions. The answer to the original question would simply be to set $k=0$ in the following. But the generalization to $k\ne 0$ is of great practical importance for wave propagation in periodic media, where you think of $U(x)$ as unwrapped along the entire $x$ axis.
Implementation with finite-difference method
To exploit the periodicity of $U(x)$ further, one could now expand it in a Fourier series and find $u_k(x)$ in Fourier space as well. However, I will instead do the calculation directly in real space because this allows me to leverage the built-in FiniteDifferenceDerivative capabilities contained in NDSolve. The essential feature relevant for this question is that FiniteDifferenceDerivative supports the option PeriodicInterpolation -> True. This makes it possible to convert the entire problem into a matrix equation using only version-8 functionality to construct all the matrices.
The equation to be solved is written here in a more generic form than in the question. The derivatives of $V(\theta)$ are really just additional terms in the periodic potential that I call $U$. So I'll only use the latter:
$$-\frac{1}{2}\frac{d^2}{dx^2}\psi_k + U(x) \psi_k = E \psi_k$$
Inserting the Bloch form, the resulting equation is (after dividing out the $\exp(i k x)$ factor):
$$-\frac{1}{2}\left(\frac{d^2}{dx^2}u_k+2i k \frac{d}{dx}u_k\right)+ \frac{k^2}{2} u_k + U(x) u_k = E u_k$$
To solve this eigenvalue equation, use the function spectrum. It's actually very short, but I tried to add a lot of comment, so the code looks a bit long:
spectrum[n_, dim_: 10000][potential_, {var_, varMin_, varMax_}, kBloch_: 0] :=
Module[
{e, v, vRange, dx, grid, potentialGrid, eKin, ePot, min, interpolate},
vRange = varMax - varMin;
interpolate =
ListInterpolation[Append[#, First[#]], {{varMin, varMax}},
PeriodicInterpolation -> True] &;
dx = N[vRange/dim];
grid = Range[varMin, varMax, dx];
eKin = -(1/2)
NDSolve`FiniteDifferenceDerivative[2, grid,
PeriodicInterpolation -> True]["DifferentiationMatrix"] -
I kBloch NDSolve`FiniteDifferenceDerivative[1, grid,
PeriodicInterpolation -> True]["DifferentiationMatrix"];
potentialGrid = Table[potential + kBloch^2/2, {var, Most[grid]}];
(* eKin is periodically interpolated,
so its last element is internally dropped by FiniteDifferenceDerivative, as redundant. Therefore,
I also have to drop the last grid element in potentialGrid. *)
min = Min[potentialGrid];
(* Matrix for the potential is shifted so its minimum entry is zero, guaranteeing that eigenvalues will be sorted in descending order: *)
ePot = DiagonalMatrix[SparseArray[potentialGrid - min]];
{e, v} = Eigensystem[eKin + ePot, -n];
(* Final step: turn vectors on spatial grid back into functions of x by interpolation: *)
Append[
Reverse /@ {e + min, Map[interpolate[#/Max[Abs[#]]] &, v]},
(* In the eigenvalues,
potential offset min was added back to get original energy scale. *)
interpolate[potentialGrid]
]]
In order to test this, let's define a potential identical to the one used in other answers (except that I defined my entire equation with a factor $1/2$ to be closer to usual quantum-mechanical choice of units):
potential[x_] = 1/2 (-D[Cos[x], x, x] + D[Cos[x], x]^2);
{eigenvals, eigenvecs, pot} = spectrum[7][potential[x], {x, -Pi, Pi}];
2 eigenvals
(*
==> {-1.80943*10^-10, 1.64737, 1.64737, 4.54059, 4.54059, 9.51758, 9.51758}
*)
In the function call, the first argument [7] is the number of desired eigenstates (starting from the lowest one). This can optionally be followed by an integer giving the grid size (default is 10000). As a second set of arguments, you provide the potential and periodicity range, in the same way you'd pass it to Plot.
My factor $1/2$ is un-done by multiplying the resulting eigenvalues by $2$ above. The results agree with the other answers for this example, also showing the degeneracies expected due to symmetry.
The above results were obtained by omitting the optional argument kBloch, which gives it the default value 0. Since the periodicity interval is $2\pi$, we'll get anti-periodic boundary conditions for $k=1/2$:
{eigenvals, eigenvecs, pot} = spectrum[7][potential[ x], {x, -Pi, Pi}, 1/2];
2 eigenvals
(*
==> {0.0204952, 1.18926, 2.73105, 2.90459, 6.77308, 6.77789, 12.7628}
*)
Plotting the results
With these results, you can make various types of plots. For each Bloch wave number $k$, you can plot the wave functions. Or you can plot the eigenvalues as a function of $k$ (that's the band structure of the periodic potential).
I'll just show the wave functions for the anti-periodic boundary condition here, but I'll wrap it into a plotting function that can be applied directly to the output of spectrum:
plotSpectrum[{ev_, efn_, potential_}, kBloch_: 0] := Module[
{evReal, evRange, eLow, eHigh, amplitude, eTicks, n = Length[ev],
n1, n2, min, max},
Check[{min, max} = First[efn[[1]]["Domain"]],
Print["plotSpectrum needs the output of spectrum as input"]; Abort[]];
evReal = Re[ev];
evRange = Last[ev] - First[ev];
amplitude = evRange/(n + 2);
eLow = First[ev] - amplitude;
eHigh = Last[ev] + amplitude;
eTicks = {{Automatic, Automatic}, {Automatic, None}};
Show[Plot[Evaluate[
Table[
ev[[i]] + amplitude Re[Exp[I kBloch x] efn[[i]][x]], {i,n}]], {x, min, max},
Epilog -> {Gray, Map[Tooltip[Line[{{min, #}, {max, #}}], #] &, Chop[ev]]},
Frame -> True, Axes -> None, FrameTicks -> eTicks,
PlotStyle -> Map[ColorData[1], Range[n]]],
Plot[potential[x], {x, min, max},
PlotPoints -> 250, ExclusionsStyle -> Automatic,
PlotStyle -> Directive[Thick, Opacity[.5], Orange],
Filling -> -10000], PlotRange -> {eLow, eHigh},
PlotRangePadding -> {{Automatic, Automatic}, {Scaled[.1], Scaled[.05]}},
FrameLabel -> {"\!\(\*FormBox[\(x\),TraditionalForm]\)",
"\!\(\*FormBox[\(\*SubscriptBox[\(ℰ\), \(n\)], \*SubscriptBox[\(ψ\), \(n\)]\),TraditionalForm]\)"}, RotateLabel -> False
]
]
plotSpectrum[{eigenvals, eigenvecs, pot}, 1/2]
Shown in orange is the potential energy. The wave functions are offset by their respective energy eigenvalues (now in my units with a factor $1/2$ in the kinetic energy). The plot only shows the real part of the wave function. The anti-periodicity can be seen in those functions whose real part doesn't happen to vanish at the boundaries. The Bloch wave number 1/2 has to be added as an argument to plotSpectrum. This is useful because it allows you to also plot the purely periodic part $u_k$ of the Bloch wave by setting $k=0$ in the last line.
Bound states of non-periodic potentials
The periodic boundary conditions implemented in spectrum become irrelevant when the potential is much higher than the desired range of eigenvalues. In that case, we're essentially dealing with bound states and the potential can have arbitrary shape. This means you can use the above code to calculate bound states in arbitrary potentials, too. Since the computation is also quite fast, it can be done in real-time inside a manipulate.
For example, consider the quartic potential with asymmetry:
Manipulate[
Show[plotSpectrum@
spectrum[10, 2000][x^4 - α x^2 + β x^3, {x, -5, 4}],
PlotLabel ->
Framed[Row[{"Potential \!\(\*FormBox[\(\[ScriptCapitalV](x)\),
TraditionalForm]\) = ", ("\!\(\*FormBox[\(x\),
TraditionalForm]\)")^4 - α ("\!\(\*FormBox[\(x\),
TraditionalForm]\)")^2 + β ("\!\(\*FormBox[\(x\),
TraditionalForm]\)")^3}], FrameStyle -> None, RoundingRadius -> 7,
Background -> Lighter[Orange]]],
{{α, 2, "Quadratic coefficient α"}, 0,
10}, {{β, 0, "Cubic coefficient β"}, 0, 2}]
So this code can do the same things that NDEigensystem does. But you have to choose the simulation interval such that the desired eigenfunctions decay to zero at the boundaries, which makes the periodic boundary conditions unnecessary. In practice, having periodic boundary conditions is useful even in this type of application, because it will continue to give reasonable results when the interval isn't small enough to make all functions decay completely to zero.
