# Numerical solution of Eigenvalue equation

I am interested in finding eigenvalues of Schrödinger-type equations, a prototype example being

$$- w^{\prime \prime } (y) - 6 \operatorname{sech}^2 (y) w(y) + w(y) = \lambda w(y)$$

I posted a question on the same equation in Solving a Sturm-Liouville problem and got great help there, both on symbolic methods as well as a numeric approach which delivers the sought result.

Out of curiosity, I have experimented also NDEigensystem and another alternative, described below, and would like to know what is that I am doing wrong with those, as I am not able to recover the same results.

Using NDEigensystem, setting Dirichlet boundary conditions at $$\pm 20$$ to mimick the real conditions at infinity, $$w(y) \to 0$$ as $$y \to \infty$$ I tried

 L = 20
{eigN2, funcN2} =
N @ DEigensystem[{w[y] - 6 Sech[y]^2 w[y] - w''[y],
DirichletCondition[w[y] == 0, True]}, w[y], {y, -L, L}, 12,
Method -> {"Arnoldi", "Criteria" -> "RealPart"}];
Plot[funcN2, {y, -20, 20}, PlotLegends -> eigN2]


I get a warning message

      ParametricNDSolveValue::bvluc: The equations derived from the boundary
conditions are numerically ill-conditioned. The boundary conditions may
not be sufficient to uniquely define a solution. If a solution is computed,
it may match the boundary conditions poorly.


Is this a working precision issue? I found a post stating that NDEigensystem cannot work with arbitrary precision. Can the issue be circumvented by a better choice of boundary conditions or what else?

In any case, I get only positive eigenvalues (which do exist and are related to oscillatory eigenfunctions, is the boundary conditions definition not satisfactory?). I used

Method -> {"Arnoldi", "Criteria" -> "RealPart"}


as I think it should sort the eigenvalues by magnitude, as I found on another post, but I am not successful (there should be an eigenvalue equal to $$-3$$), as documented on my first post on the equation linked above. The equation as written here is obtained by a change of variables from the original one, as documented in the answer by bbgodfrey.

I have also tried the approach detailed in the following post Find eigen energies of time-independent Schrödinger equation. I know the eigenfunction I am interested is even, so I use the b.c. $$w'(0) = 0$$:

U[x_] := -6*Sech[x]^2
L = 20
pfun2 = ParametricNDSolveValue[{ -ψ''[x] +
U[x] ψ[x] + ψ[x] == Ei ψ[x], ψ'[0.] ==
0., ψ[L] == 0.}, ψ, {x, -L, L}, {Ei}]
Show[Plot[U[x], {x, -5, 5}, PlotStyle -> Directive[Thick, Green]],
Plot[Evaluate@Table[Re[pfun2[Ei][x]], {Ei, -3, 2, 1}], {x, -L, L}],
PlotRange -> {All, {-6, 10}}]

val2 = Map[FindRoot[pfun2[Ei][-L], {Ei, #}] &, {1, 2, 3}]
{{Ei -> 1.02897}, {Ei -> 2.00241}, {Ei -> 3.00088}}


and there is no negative eigenvalue, nor the related eigenfunction.

Any hint would be greatly appreciated.

• What version of Mathematica do you use? I use version 12 for macos, and running the code from your first code block works just fine... Dec 23 '20 at 15:09
• I have 11.2. So you run it, without warnings and getting the $-3$ eigenvalue as well? Dec 23 '20 at 15:17
• Yes, it runs fine without any warnings or errors. And the smallest numerical eigenvalue seems to be -2.90075 which is not very far from -3. (Certainly this will come closer to -3 if you refine the discretization.) Dec 23 '20 at 15:20
• Thanks for this. The question is then, is there anything that could be done on 11.2 to get there? This is puzzling though, I wonder if the same underlying reason results in the second method I used failing. Dec 23 '20 at 15:24
• @Smerdjakov I am also got eigenvalue of -2.97433 with v.12.0.0 and without bc  DirichletCondition[w[y] == 0, True] at L=10. So let try L=10 and remove bc and option Method -> {"Arnoldi", "Criteria" -> "RealPart"} in v. 12.2 Dec 23 '20 at 15:39

Ah, this kind of $$\operatorname{sech}^2$$-based equations with decaying boundary conditions are familiar to me, as they cropped up during during my PhD (particularly this article),

And I happen to have written a package to solve eigenvalue equations of this type, using the method that I used during my PhD. This uses the Evans function, and these kind of conditions are included in my implementation.

The Evans function is an analytic function whose roots correspond to the eigenvalues. Some details are available at these two questions, my github page and this PDF, or search for CompoundMatrixMethod to see my previous answers here.

Install the package (also available on my github page):

Needs["PacletManager"]
PacletInstall["CompoundMatrixMethod",
"Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]


Needs["CompoundMatrixMethod"]
Clear[L]


First we need to convert the system into the appropriate form:

U[x_] := -6*Sech[x]^2
sys[L_] = ToMatrixSystem[-ψ''[x] + U[x] ψ[x] + ψ[x] == Ei ψ[x], {}, ψ, {x, -L, L}, Ei]


where I haven't specified a domain length or any boundary conditions; the way I have the package set up is to assume that the solution should decay if no boundary conditions are given (see the notebook of examples in the github).

Now we can evaluate the Evans function for a given value of $$Ei$$ (and domain length $$L$$), say $$Ei=-5$$ and $$L=10$$:

Evans[-5, sys[10]]
(* 0.424581 *)


This is not zero, so $$Ei=-5$$ is not an eigenvalue of the system. So we can use FindRoot:

FindRoot[Evans[Ei, sys[10]], {Ei, -5}]
{Ei -> -3.}


If we plot the Evans function, we can see that it is smooth and has the root at $$Ei=-3$$.

Plot[Evans[Ei, sys[10]], {Ei, -5, 0}]


Note that the values of this function are around 1, as the Compound Matrix Method that my package uses helps to remove the exponential growth that otherwise is present.

If you vary L, you can see that the eigenvalue converges to -3 when you get to about $$L=5$$:

Ei /. Table[FindRoot[Evans[Ei, sys[L]], {Ei, -5}], {L, 1, 10}]

(* {-3.40324, -3.01715, -3.00035, -3.00001, -3., -3., -3., -3., -3., -3.} *)


This doesn't give you the eigenvectors immediately, I have not worked out an easy way to code that into the package. But for this second order equation you can integrate from zero to L with your one boundary condition (that it is even) and an arbitrary value:

eigfun = NDSolveValue[{-ψ''[x] + U[x] ψ[x] + ψ[x] == Ei ψ[x], ψ[0] == 1, ψ'[0] == 0} /.  Ei -> (-3), ψ, {x, 0, 10}]

Plot[eigfun[x], {x, 0, 5}, PlotRange -> All]


Regarding including numerical interpolation functions, my code with the decaying BCs doesn't currently work for non-analytic functions, but if you give boundary conditions it works fine:

uu = FunctionInterpolation[U[x], {x, -10, 10},
InterpolationPoints -> 200, InterpolationOrder -> 5] // Quiet;

sysExact[L_] = ToMatrixSystem[-ψ''[x] + U[x] ψ[x] + ψ[x] == Ei ψ[x], {ψ'[0] == 0, ψ[L] == 0}, ψ, {x, 0, L}, Ei];
sysIntpol[L_] = ToMatrixSystem[-ψ''[x] + uu[x] ψ[x] + ψ[x] == Ei ψ[x], {ψ'[0] == 0, ψ[L] == 0}, ψ, {x, 0, L}, Ei];


These give the same result:

FindRoot[Evans[Ei, sysExact[10]], {Ei, -5}]
FindRoot[Evans[Ei, sysIntpol[10]], {Ei, -5}]


• thanks for pointing out this nice package, it works very well, much more elegant than my caveman approach. I am going through the theory now, the paper you suggested on your GitHub page. Just a question, is the package expected to work if the potential $U$ in the definition is not given in closed-form, but as an interpolating function coming out of a NDSolveValue solution? I get errors but I am not so sure if they are package-related or else. The code in your repository is out of my current Mathematica understanding so I cannot check it. Dec 30 '20 at 8:36
• Yeah, my code is somewhat archaic. In theory it should work with an interpolation function, although it looks like that breaks the decaying boundary conditions at the moment, because I use Limit. That is fixable though. If you try to put in the zero derivative at 0 and zero value at L, then it should work with the interpolation function Dec 30 '20 at 13:36