Problem with NDEigensystem

Question

Consider the Sturm-Liouville problem

$$y''(x) + \lambda x^2 y=0, \ y(0)=0,\ y(1)=0$$

The analytical solution is given by

$$\lambda_n=4\alpha_n^2, \ y_n(x)=\sqrt{x}J_\frac14(\alpha_nx^2)$$

where $\alpha_n$is the $n$-th zero of the Bessel function of the first kind of order 1/4. Thus, $\lambda_1\approx30.93$ and $y_1(x) \approx \sqrt{x} \,J_\frac14(2.78x^2)$

{vals, funs} = 
  NDEigensystem[
    {-(y''[x]/x^2), 
     DirichletCondition[y[x] == 0, x == 0], 
     DirichletCondition[y[x] == 0, x == 1]}, 
    y[x], {x, 0, 1}, 1]

yields vals = 16.1035 which is clearly incorrect.

Answer 1

That is an interesting case I have not seen before. You can work around this by using a time dependent formulation like so:

{vals, funs} = 
 NDEigensystem[{x^2*D[y[t, x], {t, 1}] + Laplacian[y[t, x], {x}] == 0,
    DirichletCondition[y[t, x] == 0, x == 0 || x == 1]}, {y[t, x]}, 
  t, {x, 0, 1}, 1];

vals
{30.933443250231416`}

I suspect that the original formulation has numerical issues close to x=0. Might be worth putting in the possible issues section, what do you think?

Answer 2

This is not an answer, just long comment. This verifies that eigenvalues from NDEigensystem are not correct.

ClearAll[y,x,lam];
sol=y[x]/.First@DSolve[y''[x]+lam*x^2*y[x]==0,y[x],x];
eq1=(sol/.x->0)==0;
eq2=(sol/.x->1)==0;
c=C[1]/.First@Solve[eq1,C[1]]; (*Solve for C[1]*)
eq2=eq2/.C[1]->c (*plug in second equation*)

eq2=eq2/.C[2]->1 (*divide by C[2]*)

Plot[eq2,{lam,0,200}]

There are two zeros, one near 30 and one near 140

FindRoot[eq2,{lam,30}]//Chop

4*(N@BesselJZero[1/4,1])^2

FindRoot[eq2,{lam,140}]//Chop

4*(N@BesselJZero[1/4,2])^2