# Evaluating Hough functions by using NDEigensystem on the Laplace tidal equation

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. • What are the eigenvalues you expect? – user21 Apr 24 at 10:44 ## 1 Answer 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}]} 