How to solve this 2nd-order ODE with singularity?

I tried solving the eigenvalue problem of a 2nd-order ODE $$[b^2(k-2)^2y^2-2b(k-2)(1+2ky)+4k^2+b^2(k-2)3y]f(y) \\- 3b(3by-2)f'(y)\\-(3by-2)^2f''(y)=\lambda f(y)$$ with NDEigenvalue. $$b,k$$ are just real parameters. Note that the long bracket in front of $$f(y)$$ is just a quadratic polynomial of $$y$$. Let's say $$b>0$$ for simplicity. The boundary condition is Dirichlet at infinity or at least very far away to get a completely decayed convergent solution and here I just use Dirichlet at empirically large enough yL and yR.
Thanks to @bbgodfrey, who reminded me of a singularity from $$(3by_\mathrm{sgl}-2)=0$$. I updated much of this post and prefer not showing the very naive old version.

1. It appears from the equation that the singularity renders both the 1st- and 2nd-order derivatives to contribute less and less when approaching $${y_\mathrm{sgl}}-$$, beyond which the solution (possibly??) varies little. Numerically, this holds nice when $$b$$ is small while the turning point will be smoothed and shifted further rightward when $$b$$ is larger (but why?).
2. Numerically, another important and intriguing finding is that increasing $$b$$ has the effect of squeezing the solution towards $$y_\mathrm{sgl}$$, which is not yet intuitively clear to me. Therefore, the larger $$b$$ becomes, the more unstable the numerical solution might be, I guess.
3. But what is the solution to be squeezed? Actually, a parent and simpler equation is $$[b^2(k-2)^2y^2-2b(k-2)(1+2ky)+4k^2]f(y) \\-4f''(y)=\lambda f(y),$$ whose $$n$$th solution is basically related to the $$n$$th Hermite polynomial with $$n-1$$ node(s). From the figures below, one would tentatively guess that it's this to be squeezed. Rather surprising to me, the node structure looks intact.

When $$b$$ is small like 0.03, the solution seems probably not that bad. But for larger $$b$$ like 0.5,1.0 or so, the resultant eigenvalues are very unstable against tuning the MeshOptions. Unfortunately, I don't know if there is an analytical solution. So I hope there could be some way around or even other methods to improve or reassure and it certainly also helps if one can somewhat understand the numerical observations.

    k = 0; b = 0.03; Nless = 20;
ysgl = Solve[(-2 + 3 b y) == 0, y][[1, 1]] // Last; yR = 1.3 ysgl; yL = -Sign[b] 7/Sqrt[Abs@b];
lhs = (b^2 (-2 + k) y (3 + (-2 + k) y) - 2 b (-2 + k) (1 + 2 k y) + 4 k^2) f[y]
- (-2 + 3 b y) (3 b D[f[y], y] + (-2 + 3 b y) D[f[y], {y, 2}]);
bc = DirichletCondition[f[y] == 0, True];
{vals, funcs} =
NDEigensystem[{lhs, bc}, f[y], {y, yL, yR}, Nless,
Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions"
-> {"MaxCellMeasure" -> 0.01, "MeshOrder" -> 2}}}}];
vals


The Norm of first ten solutions plotted when $$b=0.03$$. when $$b=0.06$$

Update of the analytic solution and a puzzling point

It's really nice to have an analytic solution like the one by @bbgodfrey. We can actually first shift $$y\rightarrow y+\frac{2}{3b}$$ in the ODE. Then the singularity $$y_\mathrm{sgl}=0$$. The solution is $$f(y)=\mathrm{e}^{-\frac{1}{3}(k-2)y}y^{\frac{\Delta}{9b}}u(y)$$ where $$\Delta=\sqrt{16(k+1)^2-9\lambda}$$ and $$u$$ is polynomial.

However, the exponential factor in $$f(y)$$ will diverge on the right of $$y_\mathrm{sgl}$$. It looks that analytic solution doesn't allow a physical solution vanishing at $$+\infty$$, although it's good if one sets Dirichlet b.c. exactly at $$y_\mathrm{sgl}$$. But I would expect a solution even for $$y>y_\mathrm{sgl}$$. If not, this is really surprising since the original physical meaning of y is merely a coordinate in $$(−\infty,+\infty)$$, not anything like a radius defined in $$(0,+\infty)$$.

On the other hand, the numerical solution always looks vanishing for $$y>y_\mathrm{sgl}$$. From the analytic solution, @bbgodfrey claims an upper bound of eigenvalues and higher ones are spurious. What I observe in numerics is twofold. For those 'good' eigenvalues, $$f(y_\mathrm{sgl})=0$$ seems sharp. For those 'spurious' eigenvalues, $$f(y_\mathrm{sgl})=0$$ or a turning point seems smoothed. I really have no idea if numeric solution for $$y>y_\mathrm{sgl}$$ is just an artifact or anything else.

• What are your boundary conditions? Do you want it to decay towards infinity in both directions? You have $b$ and $k$ in your equation, is one the eigenvalue and one a parameter? As you assign values to both. – KraZug Dec 19 '18 at 9:03
• Actually, the ODE is singular at (-2 + 3 b y) == 0, and this is the source of your problem, I believe., – bbgodfrey Dec 19 '18 at 11:17
• @bbgodfrey Thanks a lot for correcting me! Updated with more efforts. – xiaohuamao Dec 19 '18 at 21:20
• @KraZug Updated with more efforts and details. Yes, decay towards infinity in both directions. $b,k$ just parameters. – xiaohuamao Dec 19 '18 at 21:24
• Ugh, I hate that NDEigensystem requires to only give the left hand side of the equation, which seems so counter to everything else in Mathematica. – KraZug Dec 20 '18 at 7:32

The ODE itself can be solved symbolically. For generality, assume that k and b are undefined and that ll is the unknown eigenvalue. Then,

s = DSolveValue[lhs == ll f[y], f[y], y] // Simplify
(* E^(-(((-2 + k) (-2 + 3 b y))/(9 b))) (-6 + 9 b y)^(Sqrt[16 + 32 k + 16 k^2 - 9 ll]/(9 b))
(C[1] HypergeometricU[(-4 + 9 b - 4 k + Sqrt[16 + 32 k + 16 k^2 - 9 ll])/(9 b),
1 + (2 Sqrt[16 + 32 k + 16 k^2 - 9 ll])/(9 b), (2 (-2 + k) (-2 + 3 b y))/(9 b)] +
C[2] LaguerreL[-((-4 + 9 b - 4 k + Sqrt[16 + 32 k + 16 k^2 - 9 ll])/(9 b)),
(2 Sqrt[16 + 32 k + 16 k^2 - 9 ll])/(9 b), (2 (-2 + k) (-2 + 3 b y))/(9 b)]) *)


To proceed, apply the boundary condition that the solution vanishes at y == 2/(3 b). A bit of numerical experimentation indicates that C[1] -> 0 is required, substantially simplifying s. The second boundary condition is that the solution vanishes for large negative s. Again, some numerical experimentation indicates that the first argument of LaguerreL must be a non-negative integer, which determines the eigenvalues.

As a specific example, consider the first set of parameters in the question.

s /. C[1] -> 0 /. {k -> 0, b -> 3/100}
(* E^(200/27 (-2 + (9 y)/100)) (-6 + (27 y)/100)^(100/27 Sqrt[16 - 9 ll])
C[2] LaguerreL[-(100/27) (-(373/100) + Sqrt[16 - 9 ll]),
200/27 Sqrt[16 - 9 ll], -(400/27) (-2 + (9 y)/100)] *)


Setting the first argument of LaguerreL to an integer yields

Solve[-(100/27) (-(373/100) + Sqrt[16 - 9 ll]) == n, ll, Integers] // Flatten
(* {ll -> ConditionalExpression[(2319 + 2238 n - 81 n^2)/10000,
(n | (2319 + 2238 n - 81 n^2)/10000) \[Element] Integers && n <= 13]} *)


So, there are only 14 eigenvalues, namely

Table[{n, (-3 (-773 - 746 n + 27 n^2))/10000}, {n, 0, 13}]
(* {{0, 2319/10000}, {1, 1119/2500}, {2, 6471/10000}, {3, 519/625}, {4, 399/400},
{5, 2871/2500}, {6, 12831/10000}, {7, 876/625}, {8, 15039/10000}, {9, 159/100},
{10, 16599/10000}, {11, 1071/625}, {12, 17511/10000}, {13, 4431/2500}} *)


These values agree well with the first fourteen eigenvalues given in the question. (The rest of the numerical eigenvalues in the question are spurious.) The plot of the eigenfunction corresponding to the largest eigenvalue is, for example,

coef = s /. {C[1] -> 0, C[2] -> 1} /. {k -> 0, b -> 3/100}
/. ll -> 4431/2500 /. y -> 0
Plot[Evaluate[s/coef /. {C[1] -> 0, C[2] -> 1} /. {k -> 0, b -> 3/100}
/. ll -> 4431/2500], {y, -40, 200/9}, PlotPoints -> 1000,
ImageSize -> Large, AxesLabel -> {y, f}, LabelStyle -> {Bold, Black, 15}]


(The magnitude of coef is arbitrary, but its phase causes the function to be real, allowing it to be plotted.)

Addendum: Special case of k = 2

A number of comments on this page are associated with the question of whether eigenvalues exist for ll > 16 (k + 1)^2 /9. After substantial analysis, I am confident that there are none. To illustrate, consider the particularly simple case, k = 2. which has no eigenvalues at all, despite the fact that NDEigensystem asserts their existence.

Clear[b, k];
lhs = (b^2 (-2 + k) y (3 + (-2 + k) y) - 2 b (-2 + k) (1 + 2 k y) + 4 k^2) f[y]
- (-2 + 3 b y) (3 b D[f[y], y] + (-2 + 3 b y) D[f[y], {y, 2}]);
lhs /. k -> 2
DSolve[% == ll f[y], f[y], y] // Flatten

(* 16 f[y] - (-2 + 3 b y) (3 b f'[y] + (-2 + 3 b y) f''[y]) *)
(* {f[y] -> C[1] Cos[(Sqrt[-16 + ll] Log[-2 + 3 b y])/(3 b)] +
C[2] Sin[(Sqrt[-16 + ll] Log[-2 + 3 b y])/(3 b)]} *)


For ll > 16, neither the Sin nor the Cos term vanish at large negative y. And the two terms clearly cannot cancel, so the large negative y boundary condition cannot be satisfied. Moreover, f[y]oscillate every faster as y approaches, 2/(3 b), which also is not acceptable. Nonetheless, NDEigensystem, employed as in the question (but with yR = ysgl to avoid integrating through the singularity, and b = 1/10), predicts eigenvalues

(* {29.0571, 30.6642, 33.1095, 36.353, 40.3784, ...} *)


The eigenfunction corresponding to the first eigenvalue is

Plot[First@funcs, {y, yL, yR}, PlotRange -> All, ImageSize -> Large,
AxesLabel -> {y, f}, LabelStyle -> {Bold, Black, 15}]


Simply stated, numerical methods typically do not handle singularities well. For completeness, consider ll < 16.

lhs /. k -> 2;
Collect[DSolve[% == ll f[y], f[y], y, Assumptions -> (ll < 16)] //
Flatten // TrigToExp, Power[(2 - 3 b y), _], Simplify] /.
{(C[1] - I C[2]) -> 2 C[3], (C[1] + I C[2]) -> 2 C[4]}

(* {f[y] -> (2 - 3 b y)^(-Sqrt[16 - ll]/(3 b)) C[3] +
(2 - 3 b y)^( Sqrt[16 - ll]/(3 b)) C[4]} *)


A bounded solution as y approached 2/(3 b) requires C[3] = 0, and a bounded solution as y approaches -Infinityrequires C[4] = 0. No eigenvalues here either.

Finally, note that the leading terms from

AsymptoticDSolveValue[lhs == ll f[y], f[y], {y, 2/(3 b), 1},
Assumptions -> (y < 2/(3 b))]

(* (-(2/(3 b)) + y)^(Sqrt[16 + 32 k + 16 k^2 - 9 ll]/(9 b)) C[1] +
(-(2/(3 b)) + y)^(-Sqrt[16 + 32 k + 16 k^2 - 9 ll]/(9 b)) C[2] *)


also oscillate ever faster as y approaches 2/(3 b) when ll > 16 (k + 1)^2 /9, which is unacceptable for an eigenfunction.

• Thanks for the answer! I want to ask if there's any mathematical argument for taking the b.c. as vanishing at the singular point. Is it standard practice? Only by looking at the numerical solutions, this is not exactly fulfilled. But this might just be an artifact? – xiaohuamao Dec 20 '18 at 3:46
• @xiaohuamao Taking the boundary condition at the singularity to be zero is not necessary. In fact, you can choose whatever boundary condition you wish. It all depends on what you are trying to accomplish. I perhaps should have said that, of the two independent solutions, I chose the one that was not singular there. In fact, if you choose the smaller of the two solutions at the singular point, even if it too is singular, you obtain more eigenvalues. – bbgodfrey Dec 20 '18 at 4:55
• Let me confirm. You mean HypergeometricU is singular there, right? BTW, your two important numerical-experimentation indications are valid for general $k,b$, right? Not for special ones like $k=0$. – xiaohuamao Dec 20 '18 at 5:01
• @xiaohuamao Numerical experimental cannot be exhaustive, and my analysis certainly was not. I can, however, say with confidence that one solution is, itself, singular at the singular point, and the other vanishes there. – bbgodfrey Dec 20 '18 at 5:09
• @bbgodfrey, I can get good agreement with these eigenvalues when I apply the same bcs via the Evans function method, but they continue to be regularly spaced after that. – KraZug Dec 20 '18 at 9:14