Evaluating Hough functions by using NDEigensystem on the Laplace tidal equation

Question

Currently I am looking into the use of Mathematica to solve the classical tidal equation of M. Laplace:

$$\mathcal{F}\Theta+\gamma\Theta=0$$

whose eigenfunctions $\Theta$ are the Hough functions. Here, following the convention in the paper of Wang et al., I let the Laplace tidal operator $\mathcal{F}$ be

$$\mathcal{F}\Theta=\frac{\mathrm d}{\mathrm d\mu}\left(\frac{1-\mu^2}{\sigma^2-\mu^2}\frac{\mathrm d\Theta}{\mathrm d\mu}\right)-\frac1{\sigma^2-\mu^2}\left(\frac{s}{\sigma}\frac{\sigma^2+\mu^2}{\sigma^2-\mu^2}+\frac{s^2}{1-\mu^2}\right)\Theta$$

where $s$ and $\sigma$ are numerical parameters associated with the Hough function.

The paper by Wang et al. uses a pseudospectral approach with either normalized associated Legendre functions or Chebyshev polynomials as the basis, but I wanted to see if I can use NDEigensystem instead as it needs less programming on my part. Unfortunately my efforts have so far been in vain, and I now wonder what I am doing wrong.

As a concrete example let us take the so-called "diurnal, westward propagating, zonal wave number 1" (DW1) mode, corresponding to $s=1,\sigma=\frac12$, and take three Hough functions:

s = 1; σ = 1/2; m = 3;
tidal = -D[(1 - μ^2)/(σ^2 - μ^2) Θ'[μ], μ] +
        1/(σ^2 - μ^2) (s/σ (σ^2 + μ^2)/(σ^2 - μ^2) + s^2/(1 - μ^2)) Θ[μ];
bounds = DirichletCondition[Θ[μ] == 0, True];
{γ, Θk} = With[{h = Sqrt[$MachineEpsilon]},
               NDEigensystem[{tidal, bounds}, Θ, {μ, -1 + h, 1 - h}, m]];

(I had shifted the ends of the solution interval to avoid singularities, but I am very uncomfortable with doing this, so alternatives would be nice to see.)

Unfortunately, the eigenvalues I get from this are markedly different from the eigenvalues I got with re-implementing the pseudospectral approach (whose code I am embarrassed to share because of all the For loops it has); in particular, I am not able to reproduce figure 1 in the Wang et al. paper.

So I am interested in knowing whether I have wrongly used the functionality of Mathematica, or Mathematica is not equipped to handle the Laplace tidal equation. Ideas and other ways are of course welcome.

Answer 1

Here we need to choose eigenfunctions that satisfy the condition DirichletCondition[Θ[μ] == 0, True], and divide into even and odd.

s = 1; \[Sigma] = 1/2; m = 18; h = $MachineEpsilon;
tidal = -D[(1 - \[Mu]^2)/(\[Sigma]^2 - \[Mu]^2) \
\[CapitalTheta]'[\[Mu]], \[Mu]] + 
   1/(\[Sigma]^2 - \[Mu]^2) (s/\[Sigma] (\[Sigma]^2 + \[Mu]^2)/(\
\[Sigma]^2 - \[Mu]^2) + s^2/(1 - \[Mu]^2)) \[CapitalTheta][\[Mu]];
bounds = DirichletCondition[\[CapitalTheta][\[Mu]] == 0, True];
{\[Gamma], \[CapitalTheta]k} = 
  NDEigensystem[tidal, \[CapitalTheta], {\[Mu], -1 + h, 1 - h}, m];

lst = 3; Do[
 If[Abs[\[CapitalTheta]k[[i]][-1 + h]] <= 10^-2, 
  lst = Flatten[{lst, i}]], {i, 4, m}];

{Plot[Evaluate[
   Table[\[CapitalTheta]k[[i]][\[Mu]], {i, 
     Select[lst, EvenQ]}]], {\[Mu], -1, 1}], 
 Plot[Evaluate[
   Table[\[CapitalTheta]k[[i]][\[Mu]], {i, 
     Select[lst, OddQ]}]], {\[Mu], -1, 1}]}