# DEigensystem gives x-dependent eigenvalues

Bug introduced in Version 11 or earlier and persisting through 12.1. Reported to Wolfram Technical Support as CASE:4532301.

I am considering the eigenvalue problem associated with the double-well harmonic oscillator. Using DEigensystem,

DEigensystem[-1/2 y''[x] + (-x^2/2 + x^4/4) y[x], y[x], {x, -∞, ∞}, 2]


gives eigenvalues that depend on $$x$$:

{Sqrt[-2 + x^2]/(2*Sqrt[2]), (3*Sqrt[-2 + x^2])/(2*Sqrt[2])

Maybe Mathematica 11 cannot solve for this problem, but why does it provide x-dependent eigenvalues?

EDIT: After bbgodfrey's answer, I'd like to go deep into this problem. The solutions are HeunT functions and the eigenvalues exist. In Dong 2019 you have an in-depth reference to this problem. Mathematica should be able to solve in version 12.1 since it has Heun triconfluent functions.

• Thanks for the reference, the abstract of which is consistent with my answer below. Unfortunately, I do not have access to the article itself. Does it give symbolic solutions for the eigenvalues too, or only for the eigenfunctions? By the way, in my experience DSolve, on which DEigensystem depends, fails to solve many ODEs for which solutions are known.. What troubles me more, though, is that DEigensystem gives spurious answers for {x, -∞, ∞}, instead of returning unevaluated. Perhaps, it is time to report this as a bug. – bbgodfrey Apr 15 '20 at 12:45
• @bbgodfrey They give symbolic solutions for the eigenfunctions, but the eigenvalues are computed numerically. I have the paper, I would not attach here but I can send it to you, if you're interested. I guess reporting a bug would be useful for other users. – Fabio Apr 15 '20 at 14:24
• I am happy to report it as a bug and add the usual bug header, if you like. – bbgodfrey Apr 15 '20 at 15:00
• ok, let's do it! :) – Fabio Apr 15 '20 at 15:18

Version 12.1 also give this strange result, which may be a bug.

DEigensystem[{-1/2 y''[x] + (-x^2/2 + x^4/4) y[x]}, y[x], {x, -∞, ∞}, 2]
(* {{Sqrt[-2 + x^2]/(2 Sqrt[2]), (3 Sqrt[-2 + x^2])/(2 Sqrt[2])},
{E^(-((x^2 Sqrt[-2 + x^2])/(2 Sqrt[2]))),
2^(3/4) E^(-((x^2 Sqrt[-2 + x^2])/(2 Sqrt[2]))) x (-2 + x^2)^(1/4)}} *)


On the other hand, reducing the limits to

DEigensystem[{-1/2 y''[x] + (-x^2/2 + x^4/4) y[x]}, y[x], {x, -4, 4}, 2]


returns unevaluated. It is, however, easy to solve this problem numerically.

snn = NDEigensystem[{-1/2 y''[x] + (-x^2/2 + x^4/4) y[x]}, y[x], {x, -4, 4}, 4];
snn // First
Plot[Evaluate[snn // Last], {x, -4, 4}, PlotRange -> All]
(* {0.147275, 0.872551, 2.12949, 3.59638} *)


Returning now to the symbolic solution problem posed in the question, consider

s = DSolveValue[-1/2 y''[x] + (-x^2/2 + x^4/4) y[x] == lamda y[x], y[x], x]
(* E^((x (3 - x^2))/(3 Sqrt[2])) C[1]^2
HeunT[-(1/2) - 2 lamda, -Sqrt[2], Sqrt[2], 0, -Sqrt[2], x] +
E^(-((x (3 - x^2))/(3 Sqrt[2]))) C[2]
HeunT[-(1/2) - 2 lamda, Sqrt[2], -Sqrt[2], 0, Sqrt[2], x] *)


The HeunT functions, newly defined in Version 12.0, are not bounded for large Abs[x], which may confuse Mathematica. However, if DEigensystem, cannot solve this problem, it should return unevaluated. As noted above, this may be a bug. The first eigenvalue can, however, be obtained by

sr1 = FindRoot[0 == D[s, x] /. {C[1]^2 -> 1, C[2] -> 1, x -> 4},
{lamda, 0.15}, WorkingPrecision -> 45]
(* {lamda -> 0.147235140084093444055886856920046970641521561} *)


which agrees with the first eigenvalue determined by NDEigensystem above to four significant figures. (Note that FindRoot does not converge for smaller WorkingPrecision.) Likewise, the second eigenvalue can be obtained by

sr2 = FindRoot[0 == D[s, x] /. {C[1]^2 -> 1, C[2] -> -1, x -> 4},
{lamda, 1}, WorkingPrecision -> 45]
(* {lamda -> 0.872261197867424491236128918505849507397539635} *)


Larger even and odd eigenvalues also are given by sr1 and sr2, respectively, with larger initial guesses for lamda. Substituting the values for {C[1]^2, C[2], lamda} into s then yields the same curves as above, up to a normalization factor. (Plot requires WorkingPrecison -> 30 for smooth curves.) So, with assistance Mathematica can obtain symbolic solutions for the eigenfunctions, although not for the eigenvalues.

I would welcome the readers' views on whether the behavior of DEigensystem is a bug. (Trace produces very lengthy but not particularly informative output.)

Addendum - Computation Using New Feature of DSolve 12.1

Documentation for DSolve 12.1 (under Scope) describes how to solve Sturm-Liouville problems. It can be applied to the present problem as follows.

newds = DSolveValue[{-1/2 y''[x] + (-x^2/2 + x^4/4) y[x] == lamda y[x],
y'[-4] == 0, y'[4] == 0}, y[x], x, Assumptions -> 0 < lamda < 4];


producing a lengthy Piecewise function that contains the eigenfunction with one constant of integration eliminated and a transcendental equation for lamda (not the actual eigenvalues in this case). This equation is extracted by newds[[1, 1, 2, 1]], which is solved without difficulty for the eigenvalues.

Table[FindRoot[newds[[1, 1, 2, 1]] /. C[1] -> 1, {lamda, n},
WorkingPrecision -> 45], {n, .2, 3.2, 1}] // Values // Flatten // N[#, 6] &
(* {0.147235, 0.872261, 2.12798, 3.59109} *)


as expected.

• Hi, I've added a reference to the analytical solutions of this problem. your analysis is interesting. Mathematica should solve this problem analytically and provide the correct eigenvalues, since the solutions are Heun triconfluent functions – Fabio Apr 15 '20 at 7:15
• What is the reason to put c[1],c[2]$=\pm1$ in the code for sr1,sr2? I feel they can have different magnitude in general? – xiaohuamao Oct 7 '20 at 23:08
• @xiaohuamao I chose C[1] = C[2] and C[1] = - C[2], because I knew that doing so would give symmetric and antisymmetric solutions, consistent with the numerical plot in my answer. To verify this, run FindRoot[{0 == D[s, x] /. x -> 4, 0 == D[s, x] /. x -> -4} /. C[1] -> 1, {{lamda, 0.15}, {C[2], 2}}, WorkingPrecision -> 45], which returns lambda as before and C[2} -> 1.. Similarly, FindRoot[{0 == D[s, x] /. x -> 4, 0 == D[s, x] /. x -> -4} /. C[1] -> 1, {{lamda, 1}, {C[2], 2}}, WorkingPrecision -> 45] returns the second eigenvalue as before and C[2} -> -1.. – bbgodfrey Oct 8 '20 at 1:13
• Thanks. Now I see this is partially because the differential operator is parity even and the Sturm-Liouville theorem. – xiaohuamao Oct 8 '20 at 5:09

The question is still ill-posed. This is most probably not a sensation that the eigensystem is dependent on x. It is a sensation the n in the built-in is not chosen appropriately and as suggested by the Mathematica documentation. Prefer n=4 instead of n=2. The reason this is like the Laplacian a second-order differential equation with on x relying coefficient functions.

So the input is:

DEigensystem[-1/2 y''[x] + (-x^2/2 + x^4/4) y[x],
y[x], {x, -∞, ∞}, 4]


With my version 12.0.0 even n=2 will pose the results shown by others.

The main problem is that the interval between {-2,2} has no real values they are all imaginary. Therefore the system can not be properly solved on the complete interval entered in the DEigensystem.

So the result in the imaginary part {-2,2} is displayed with

Plot[Evaluate[Im@%28[[2]]], {x, -10, 10}, PlotRange -> Full]


The real part:

Plot[Evaluate[%28[[2]]], {x, -10, 10}, PlotRange -> Full]


This is in Mathematica 12.0.0. Increase the n option in yours might be appropriate.

• To what does %28 refer? – bbgodfrey Apr 15 '20 at 21:21
• dear @user2432923, honestly I don't understand your comment at all and how can be useful your suggestion in providing n=4 to find the correct solution. I also claim that the problem is well-posed. The eigenvalues cannot depend on the spatial variable, furthermore from a software perspective we have provided that symbol and the software should interpret this information accordingly. I think also that DEigensystem should give the correct result at least for the eigenfunctions since DSolve does. Likely, the eigenvalues cannot be computed in an analytical way... DEigensystem seems to be buggy. – Fabio Apr 16 '20 at 7:37