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


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.

  • 2
    $\begingroup$ What are the eigenvalues you expect? $\endgroup$
    – user21
    Apr 24, 2019 at 10:44

1 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}];

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



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.