# NDEigensystem to solve differential equation

When trying to solve the differential equation

radialEqdouble[k_] = f''[u] + k*u^2*f[u] - u^4*f[u]


where k is a constant, I am able to produce a plot consistent with the expected double well wavefunction.

radialEqdouble[k_] = f''[u] + k*u^2*f[u] - u^4*f[u]
Simplify[radialEqdouble[k] /. f -> (ψ[ArcTan[#]] &) /.
u -> (Tan[ξ]), Pi/2 > ξ > -Pi/2]
{evdouble4, efdouble4} =
DirichletCondition[ψ[ξ] == 0,
True]}, ψ[ξ], {ξ, -Pi/2, Pi/2}, 1,
Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}},
"Eigensystem" -> {"Arnoldi", MaxIterations -> 80000}}];
fdoub4[x_] = efdouble4[[1]] /. ξ -> x
efdoub4[u_] = efdouble4[[1]] /. ξ -> ArcTan[u]


Plotting this:

Plot[{radξdoub4[fdoub4] -
evdouble4[[1]]*fdoub4[ξ], -evdouble4[[1]]*
fdoub4[ξ]}, {ξ, -Pi/2, Pi/2}, PlotRange -> All]


Properly produces

But when I increase the coefficient k to any value greater than 8, as depicted here:

{evdouble10, efdouble10} =
DirichletCondition[ψ[ξ] == 0,
True]}, ψ[ξ], {ξ, -Pi/2, Pi/2}, 1,
Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}},
"Eigensystem" -> {"Arnoldi", MaxIterations -> 80000}}];
fdoub10[x_] = efdouble10[[1]] /. ξ -> x
efdoub10[u_] = efdouble10[[1]] /. ξ -> ArcTan[u]


and then plot,

Plot[{radξdoub10[fdoub10] -
evdouble10[[1]]*fdoub10[ξ], -evdouble10[[1]]*
fdoub10[ξ]}, {ξ, -Pi/2, Pi/2}, PlotRange -> All]


It produces this plot:

Why is the double well wavefunction behavior lost once the coefficient becomes greater than 8? Why do the error bars blow up? I feel like it might have something to do with my "Method" options in NDEigensystem but have been unable to locate the issue. This is not the behavior I would expect out of NDEigensystem, as the left hand and right hand sides have way too large of a difference between them.

• When you define the following functions using set delay it works: rad[Xi]doub10[q_] := radial[Xi]double[10] /. [Psi] -> q fdoub10[x_] := efdouble10[[1]] /. [Xi] -> x efdoub10[u_] := efdouble10[[1]] /. [Xi] -> ArcTan[u] Somewhere precision is lost, however, I can not explain why the computation inside Plot is more precise than if you precompute it. Has someone an explanation? – Daniel Huber Oct 28 at 10:19
• @DanielHuber Changing to set delayed doesnt seem to work, as the efdouble10 wavefunction is still not a double well wavefunction in set delayed, its actually computed before even changing those functions to set delayed. So I really think its an issue with my use of NDEigensystem – BOUNCE Oct 28 at 16:18
• What should it be? I get a symmetrical function, but not the ground state though. What do you expect? – Daniel Huber Oct 28 at 20:26
• @DanielHuber I am expecting a ground state double well wavefunction, like the first image when k=4 – BOUNCE Oct 28 at 21:05
• According to the documentaton, "Eigenvalues are sorted in order of increasing absolute value." For k greater than about 4.1, the first symmetric eigenfunction has a corresponding eigenvalue that is smaller in absolute value than that of the first antisymmetric eigenfunction. Compute the first two eigenfunctions to see what I mean. I can write this as an answer with a bit more detail, if that would be helpful. – bbgodfrey Oct 28 at 21:23

The antisymmetric eigenfunction in the question is replaced by a symmetric eigenfunction in the question as k is increased due to the way that NDEigensystem sorts eigenvalues.

To see this, compute the first two solutions in each case, instead of just the first one. For k = 4,

{evdouble, efdouble} = NDEigensystem[{radialξdouble[4],
DirichletCondition[ψ[ξ] == 0, rue]}, ψ[ξ], {ξ, -Pi/2, Pi/2}, 2,
Method -> {"SpatialDiscretization" -> {"FiniteElement", {"MeshOptions" ->
{"MaxCellMeasure" -> 0.001}}}, "Eigensystem" -> {"Arnoldi", MaxIterations -> 80000}}];
evdouble
(* {1.24792, 1.71035} *)
Plot[efdouble, {ξ, -Pi/2, Pi/2}]


However, for k = 8, the corresponding solutions are

{0.13202, -1.41966}


So, symmetric and antisymmetric solutions are computed for both values of k. However, using the rule "Eigenvalues are sorted in order of increasing absolute value", NDEigensystem lists the antisymmetric eigenfunction first for k = 4 and second for k = 8. Solving for only one eigenfunction, as in the question, then creates the appearance that the antisymmetric eigenfunction has been replaced by the symmetric one.

Addendum: Finding the "right" antisymmetric eigenfuntion

If, however, we are seeking not just any antisymmetric eigenfunction but the one with no oscillations, more searching must be done. First, use ψ[0] == 0, as a boundary condition to eliminate all symmetric eigenfunctions. Then search for the eigenfunction with the smallest Abs[ψ'[0]] to identify the one with fewest oscillations. This is accomplished for k = 4 by

{evdouble, efdouble} = NDEigensystem[{radialξdouble[4],
DirichletCondition[ψ[ξ] == 0, True]}, ψ[ξ], {ξ, 0, Pi/2}, 4,
Method -> {"SpatialDiscretization" -> {"FiniteElement", {"MeshOptions" ->
{"MaxCellMeasure" -> 0.001}}},
Eigensystem" -> {"Arnoldi", MaxIterations -> 80000}}];
efdouble[[Ordering[Abs[D[efdouble, ξ] /. ξ -> 0], 1]]]/Sqrt[2];
Plot[Evaluate[-Sign[D[%, ξ] /. ξ -> 0]*%], {ξ, 0, Pi/2}]


which picks out the first eigenfunction provided by NDEigensystem, as expected. And, for k = 8, the code yields

which is the fourth eigenfunction provided by NDEigensystem. Larger values of k would require searching over even more eigenfunctions.