I'm trying to get the eigenvalues of a one dimensional time-independent Schrödinger equation,
$-\frac{h^2}{2m_0}\frac{d^2\psi}{dx^2}+U(x)~\psi=Ei~\psi$
where U(x) is some potential and Ei is the eigen energy. I'm trying to find the Ei if given a potential U(x).
For example, if U(x) is the harmonic potential
$U(x)=\frac{1}{2}m_0\omega^2x^2$
then the solution for Ei would be
$Ei=(n+\frac{1}{2})\hbar\omega$.
Here is my try to solve a harmonic potential:
There is a new function introduced in version 9 called ParametricNDSolveValue which can be used to solve parametric differential equations. In its advertisement, there is a demonstration of solving a similar eigenproblem of the finite potential well. There the energy is treated as a parameter, and then the Schrödinger is solved as a parameter of the energy. Then we select only the energies that make the wave function vanish at large distance, which is required by physics. Those selected energies are the eigen energy of the system. Here I use the same way.
\[HBar] = e = m0 = 1.; ω = 1;
U[x_] := 1/2 m0 ω^2 x^2
The potential is symmetric, so that the wave function must be odd or even function.
First if the wave function is odd, then it must goes through zero.
pfun1 = ParametricNDSolveValue[{-(\[HBar]^2/(2 m0)) ψ''[x] +
U[x] ψ[x] == Ei ψ[x], ψ[0.] == 0., ψ[1.] ==
1.}, ψ, {x, -4, 4}, {Ei}]
These are the solutions
Show[Plot[U[x], {x, -5, 5}, PlotStyle -> Directive[Thick, Green]]
, Plot[Evaluate@Table[Re[pfun1[Ei][x]], {Ei, 0., 5., 0.5}], {x, -4,
4}], PlotRange -> {All, {-6, 10}}]
The wave function value at -4 as a function of parameter Ei
Plot[pfun1[Ei][-4], {Ei, 0, 5}, ImageSize -> Medium]
We can select the energies which make the wave function goes to zero, those are the eigen energies. (Note that previously this just said val instead of val1, which gave an error down the road.)
val1 = Map[
FindRoot[pfun1[Ei][-4], {Ei, #}] &, {1, 2, 3}]
(*{{Ei -> 1.50001}, {Ei -> 1.50001}, {Ei -> 3.50169}}*)
Similarly, if the wave function is even, then the first derivative should be 0 at 0.
pfun2 = ParametricNDSolveValue[{-(\[HBar]^2/(2 m0)) ψ''[x] +
U[x] ψ[x] == Ei ψ[x], ψ'[0.] == 0., ψ[1.] ==
1.}, ψ, {x, -4, 4}, {Ei}]
Show[Plot[U[x], {x, -5, 5}, PlotStyle -> Directive[Thick, Green]]
, Plot[Evaluate@Table[Re[pfun2[Ei][x]], {Ei, 0., 5., 0.5}], {x, -4,
4}], PlotRange -> {All, {-6, 10}}]
Plot[pfun2[Ei][-4], {Ei, 0, 5}, ImageSize -> Medium]
The eigen energies can be selected in the same way
val2 = Map[FindRoot[pfun2[Ei][-4], {Ei, #}] &, {1, 2, 3}]
(*{{Ei -> 0.5}, {Ei -> 2.5002}, {Ei -> 2.5002}}*)
So finally the eigen energyies are
eigenEv = val1 /. Rule[_, x_] -> x // Flatten;
eigenOd = val2 /. Rule[_, x_] -> x // Flatten;
Sort[Join[eigenEv, eigenOd]]
(*{0.5, 1.50001, 1.50001, 2.5002, 2.5002, 3.50169}*)
eigen functions are
Show[Plot[U[x], {x, -5, 5}, PlotStyle -> Directive[Thick, Green]]
, Plot[Evaluate@{Table[Re[pfun1[Ei][x]], {Ei, eigenEv}],
Table[Re[pfun2[Ei][x]], {Ei, eigenOd}]}, {x, -4, 4}],
PlotRange -> {All, {-3, 3}}]
which basically agree with the analytical results.
Questions:
Are there other ways to find the eigen energies which are more efficient and accurate than the parameter methods? For example, in this case we can see that the fourth eigen energy should be exactly 3.5 but it gave 3.50169 which is quite a large error. Also in my real system, I'm interested in calculating up to the hundredth eigen energies, this method seems too slow to do that.
If the potential U(x) is not symmetric, how do I set the initial or value? For example, in the harmonic potential case, we set
ψ'[0]=0, ψ[1]=1in the even case andψ[0]=0, ψ[1]=1in the odd case, in a more general potential, how do I set the initial or boundary value?
FindRootline in the code and let it return more roots. $\endgroup$