On computation of Hough functions

Hough functions are the eigenfunctions of the Laplace’s tidal equation governing fluid motion on a rotating sphere with a resting basic state. Several numerical methods have been used in the past. In this paper, we compare two of those methods: normalized associated Legendre polynomial expansion and Chebyshev collocation. Both methods are not widely used, but both have some advantages over the commonly-used unnormalized associated Legendre polynomial expansion method. Com5 parable results are obtained using both methods. For the first method we note some details on numerical implementation. The Chebyshev collocation method was first used for the Laplace tidal problem by Boyd (1976) and is relatively easy to use. A compact MATLAB code is provided for this method. We also illustrate the importance and effect of including a parity factor in Chebyshev polynomial expansions for modes with odd zonal wavenumbers.


Introduction
Hough functions are the eigenfunctions of the eigenvalue problem of the following form: where F is a linear differential operator, the Laplace tidal operator, defined as with µ = sin φ ∈ [−1, 1], φ being the latitude, s the zonal wave number, and σ the dimensionless frequency normalized by 2 ( the earth's rotation rate), while is Lamb's parameter (Andrews et al., 1987, p. 154), with a being the earth's radius, g the acceleration due to the earth's gravity, and h the so-called equivalent depth.
Several numerical methods have been used to solve the eigenvalue problem for the Laplace tidal equation in the past.Hough (1898) pioneered the solutions of the Laplace tidal equations using spherical harmonic expansion, or equivalently spherical harmonic Galerkin method, so eigenfunctions of the eigenvalue problem (Eq. 1) that describe the latitudinal dependence are often called Hough functions (Flattery, 1967;Longuet-Higgins, 1968;Lindzen and Chapman, 1969).Each function of latitude and longitude is expanded as a Fourier series in longitude using the usual Fourier functions, cos(sλ) and sin(sλ), where s, an integer, is the "zonal wave number", and λ is the longitude.Each longitudinal trigonometric function is multiplied by a latitudinal basis function which depends on the zonal wave number s. Hough and his successors used a latitudinal basis of unnormalized associated Legendre polynomials (ALPs).Both Kato (1966) and Flattery (1967) used the method of continued fractions to solve for eigenvalues one by one with iterations.This is not the most convenient method to work with, and some eigenvalues could be missed.Chen and Lu (2009) also discussed calculation of Hough functions following the same original formulation without showing any details on numerical procedures.

H. Wang et al.: Computation of Hough functions
Computation of Hough functions based on expansion in terms of normalized ALPs was first used by Dikii (1965).It was later elaborated in a note by Groves (1981), along with a method of evaluating related wind functions.Jones (1970) used group-theoretical methods to obtain a matrix representation of Hough functions by expanding in normalized spherical harmonics.
Although it is closely related to the original method of expansion in terms of unnormalized ALPs, expansion in terms of the normalized ALPs leads to two symmetric matrices for symmetric and anti-symmetric modes.This has both computational and conceptual advantages over the original expansion in unnormalized ALPs: (1) the eigenvalue problem of symmetric matrix can be solved very accurately by the Jacobi method (e.g., Demmel and Veselić, 1992), and (2) symmetry guarantees that all of the "eigenvalues are real and that there is an orthonormal basis of eigenvectors" (Golub and Van Loan, 1996, p. 393).
There is also another way of computing Hough functions or global normal modes, such as Longuet-Higgins (1968), Kasahara (1976), andŽagar et al. (2015), also using spherical harmonic expansion, in which the equivalent depth is assigned (for each zonal wave number) and the frequency of the normal modes is obtained as the eigenvalues.This is different from the eigenvalue problem for tidal waves in which the wave frequencies and zonal wave number are specified and eigenvalues are obtained and used to compute equivalent depths, just as stated in the original eigenvalue problem (Eq.1).
The collocation method was first applied to compute Hough functions by Boyd (1976).His latitudinal basis functions replace associated Legendre functions by cosine functions of colatitude ϕ multiplied by a "parity factor" which is sin(ϕ) for odd zonal wave number s and the constant one for even zonal wave numbers.The parity factor is explained in Appendix C. In addition, the modified latitudinal variable is often used to analyze and solve differential equations in spherical geometry.The reason is that trigonometric functions are replaced by powers of µ, simplifying almost everything.And denoting the Chebyshev polynomials by T n (x), Chebyshev's famous identity shows that T n (µ) = T n (cos(ϕ)) = cos(nϕ), n = 0, 1, . ... Thus a Fourier cosine series in colatitude is, with the same coefficients, also a Chebyshev polynomial series in µ.Boyd (1976) and Orszag (1974) listed several advantages of Chebyshev polynomial collocation over spherical harmonic Galerkin approximations.First, cosines/Chebyshev polynomials are much simpler than associated Legendre functions, which are different for each different zonal wave number s.Second, collocation, which evaluates and interpolates, is much easier to program than the Galerkin method, which integrates.These advantages make it much easier to apply the Chebyshev collocation method than the spherical harmonic Galerkin method.See also Hesthaven et al. (2007, Chapter 3) for a discussion of advantages of Fourier collocation methods over the Fourier-Galerkin methods.
In this paper we compare the solution of the eigenvalue problem for the Laplace tidal operator using two numerical methods, the normalized ALP expansion method and the Chebyshev collocation method.Both methods are not widely used, but both have some advantages over the commonly used unnormalized ALP expansion.For the first method we note some details of numerical implementation as the denominators in some terms of matrix entries can become zero.For the second method a compact MATLAB code is provided to facilitate its use.We also discuss other related issues and show that there is no accuracy penalty in using the Chebyshev collocation method.

Computation of Hough functions
In this section, we compare two methods for computing Hough functions: one using the normalized ALP expansion, the other using the Chebyshev collocation method.

Computation of Hough functions using normalized associated Legendre polynomial expansion
The first method uses the expansion in terms of normalized ALPs (e.g., Groves, 1981).To solve the Laplace tidal equation, first expand in terms of the unnormalized associated Legendre polynomials P s r : Substituting into the Laplace tidal equation Eq. ( 1), one obtains where and These equations were first given by Hough (1898); see also Lindzen and Chapman (1969).The normalized associated Legendre polynomials P r,s are defined in terms of the unnormalized associated Legendre polynomials P s r by Expanding in terms of the normalized associated Legendre polynomials P r,s , we have (Dikii, 1965;Groves, 1981) where .
Equation ( 10) can be written in a matrix form for the coefficients vector x = [a s , a s+1 , a s+2 , a s+3 , . ..]T as the matrix eigenvalue problem F 0 x = λx, with matrix F 0 defined as Or it may be written as, respectively, [a s , a s+2 , . ..]T for symmetric modes, with matrix F 1 defined as and F 2 x 2 = λ 2 x 2 , x 2 = [a s+1 , a s+3 , . ..]T for anti-symmetric modes, with matrix F 2 defined as These are real symmetric matrices, and the eigenvalue problem can be solved accurately using the Jacobi methods (e.g., Golub and Van Loan, 1996, Chapter 8).The computed eigenvectors are the expansion coefficients.
A few remarks on unnormalized versus normalized ALP expansion are in order here.The unnormalized polynomials (not just ALPs, but Legendre and Chebyshev and Hermite polynomials too) have survived because the canonical unnormalized forms have polynomial coefficients that are integers or rational numbers.This is convenient for many applications, such as when using exact arithmetic in computer algebra.Note that this property carries over to the Galerkin matrix elements for the Hough differential equation, which are rational functions of r and s in Eq. ( 6).Also, for some purposes it is very convenient to use polynomials which are all 1 at µ = 1, as true for unnormalized Chebyshev and Legendre polynomials.The bad news is that unnormalized polynomials generate bigger roundoff errors in all calculations, not just computing matrix eigenvalues.The Galerkin matrix element formulas are more complicated for normalized polynomials.As we noted above, a particular advantage of working with normalized ALPs is that the discretization matrix becomes a symmetric matrix.Spectral discretizations often generate a few inaccurate eigenvalues with nonzero imaginary parts, but the eigenvalues of a symmetric tridiagonal matrix are always real.
A note on numerical implementation is relevant here, since denominators of terms in M r can become zero.We found that Eq. (6b), instead of Eq. ( 11b), of M r should be used, even though the two forms are equivalent.In addition, we should set that last term of Eq. (6b) of M r to zero when it becomes a form of 0/0.Thus, to compute the (s = 2, σ = 1) modes or SW2 (semidiurnal, westward propagating, zonal wave number 2) modes, we should set the last term of Eq. (6b) to zero when r = s = 2.
The Fortran 90 source code of the Jacobi eigenvalue algorithm implemented by Burkardt (2013) can be used to solve the two symmetric matrix eigenvalue problems.It can actually, for the (s = 1, σ = 0.5) modes or DW1 (diurnal, westward propagating, zonal wave number 1) tide, compute the one infinite eigenvalue with P 2,1 as the eigenmode, "the most important odd mode" (Lindzen and Chapman, 1969, p. 151) since P 2,1 ∝ sin φ cos φ.So in this way we will not miss any important eigenvalue or eigenfunction; see Sect. 3 for a discussion on the "missing" modes for the solar diurnal modes and the completeness of Hough functions.When using MAT-LAB, we can set any inf matrix entry to realmax and then www.geosci-model-dev.net/9/1477/2016/Geosci.Model Dev., 9, 1477-1488, 2016 use the MATLAB function eig to solve the matrix eigenvalue problem.It is also preferable to compute eigenvalues for symmetric and anti-symmetric modes separately, especially when there are interior singularities, e.g., for the DW1 tide.A MATLAB implementation is shown in Appendix B1.
Using the method of expansions in the normalized associated Legendre polynomials, truncated at r max = 60 on 94 Gaussian quadrature points, we compute eigenvalues and eigenfunctions for several important solar tides.We use solar day instead of sidereal day in our computations.The first several equatorial symmetric and anti-symmetric modes for DW1 are shown in Fig. 1.The first several equatorial symmetric and anti-symmetric modes for SW2 of scalar fields are shown in Fig. 2a-b.The first several equatorial symmetric and anti-symmetric modes for (s = 3, σ = 1.5) modes or TW3 (terdiurnal, westward propagating, zonal wave number 3) for temperature field are shown in Fig. 3.For completeness, a method of computing Hough functions for the horizontal wind components by Groves (1981) (with correction) is presented in Appendix A.

Computation of Hough functions using Chebyshev collocation method
The Chebyshev collocation method was first used by Boyd (1976) to solve the Laplace tidal problem.Expand in terms of the Chebyshev polynomials T n (µ): which includes a parity factor sin ϕ for the odd zonal wave number s (Orszag, 1974;Boyd, 1978), where ϕ is colatitude, ϕ = π/2−φ.See Appendix C for an explanation of the parity factor.The Chebyshev collocation points can be defined in different ways.When the interior or "roots" points are used, they are defined as (e.g., Boyd, 2001, p. 571) where N is total number of collocation points.By using the differential matrices, it is straightforward to apply the Chebyshev collocation methods to any differential operators.Discussion on property of Chebyshev polynomials and collocation method can be found in Boyd (2001) and Trefethen (2000).A MATLAB implementation is shown in Appendix B2.Parity requirement is discussed in Orszag (1974).To quote from Orszag (1974), If parity requirements are violated, then differentiability is lost (at the boundaries, i.e., at the poles), possibly resulting in slow convergence of series expansions and associated Gibbs' phenomena.It is important that assumed spectral representations not impose an incorrect symmetry on a solution if infinite-order accurate results are desired (see also Boyd, 1978).To show how accuracy is affected by the parity factor, we compare the eigenfunction expansion coefficients b n computed with or without the parity factor in Fig. 4. For both  terdiurnal and pentadiurnal tides, when the parity factor is removed, only limited lower-order algebraic convergence rates are achieved: fourth order for terdiurnal and seventh order for pentadiurnal.When the parity factor is included, spectral or exponential convergence is restored.Thus including the parity factor improves the accuracy dramatically, so solutions are less affected by singularities when they exist.It is important to include the parity factor when computing eigen-values and eigenfunctions for DW1 (s = 1, σ = 0.5) modes (see Sect. 2.3).A theoretical justification for the parity factor is given in Appendix C.
The MATLAB code listed in Appendix B2 includes a parity factor for the odd zonal wave number.It also computes Hough modes for horizontal wind components.The computed eigenvalue in this case is just (negative) γ , and from Eq. ( 3) we can compute the corresponding equivalent depths  h.Hough functions are simply the computed eigenvectors (with different normalization factors that are irrelevant) when Chebyshev differential matrices are used.So the eigenvalue and eigenvector problem we solve can be viewed as a direct discretization of the original operator eigenvalue problem (Eq.1).

Comparison of the two methods
Table 1 compares the number of good eigenvalues that can be obtained using the two methods.The "good" eigenvalue is defined as one whose relative error is less than 10 −6 , where λ is the eigenvalue computed at high truncation N = 160.This definition is somewhat arbitrary but is useful for comparisons.It shows that for DW1 about 60 % of the computed eigenvalues are good using the normalized ALP expansion method and about 50 % of the computed eigenvalues are good using the Chebyshev collocation method; for SW2 a little over 50 % of the computed eigenvalues are good using both methods; and for TW3 the number of good eigenvalues is about 75 % for both methods.We note that for DW1 only about 15 % of the computed eigenvalues are good without parity factor, contrasted to 50 % with parity factor.This again illustrates the importance of preserving correct parity.Considering the "unusual difficulties" in solving the eigenvalue problem of the Laplace tidal equation using general numerical methods, as remarked by Bailey et al. (1991), it is remarkable that the Chebyshev collocation method with a parity factor for odd zonal wave number can be used so successfully in solving the eigenvalue problem of the Laplace tidal equation.

A remark on the completeness of Hough functions
Although the completeness of Hough functions for zonal wave number s and period T = (s + 1)/2 days was questioned earlier by Lindzen (1965), completeness was later proved by Holl (1970) with further analysis by Homer (1992).Giwa (1974) proved by direct computation that, for zonal wave number s and period T = (s + 1)/2 days, Hough functions for tidal oscillations are the same as the associated Legendre polynomials P s s+1 , and Hough functions form a complete set of orthogonal functions.
One advantage in using the normalized associated Legendre polynomials as basis functions, as shown in Sect.2.1, is that the eigenvalue problem becomes an eigenvalue problem for two real symmetric matrices: one for symmetric modes and one for anti-symmetric modes.The spectral theory of (Hermitian) symmetric matrices tells us that these real symmetric matrices have "a complete set of orthogonal eigenvectors, and that the corresponding eigenvalues are real" (e.g., Lax, 2002, Chapter 28).Thus this approach in a heuristic way shows the completeness of Hough functions.

Summary and conclusions
In this paper, we briefly survey the numerical methods for computing eigenvalues and eigenvectors for the Laplace tidal operator.In particular we compare two numerical methods: the normalized ALP expansion and Chebyshev collocation.The normalized ALP expansion method leads to two symmetric matrices which can be solved very accurately.It also has an advantage in providing another conceptual understanding for the completeness of eigenfunctions (Hough functions) of the Laplace tidal operator.We also note some details on numerical implementation and provide a MAT-LAB code.
The Chebyshev collocation method was first used by Boyd (1976) for computing the eigenvalues for the Laplace tidal problem.Here we compared this method with the ALP expansion and found that both are producing comparable results.The Chebyshev collocation method uses Fourier cosine series in colatitude as the basis functions and is relatively easy to work with.A compact MATLAB code is provided to facilitate the use of the Chebyshev collocation method for the Laplace tidal problem.
The Chebyshev polynomial expansion method is merely a Fourier cosine expansion method in disguise (Boyd, 2001).In using the Chebyshev collocation method, it is important to include a parity factor in Chebyshev polynomial expansion for odd zonal wave number modes, as illustrated in the paper.

Figure 1 .
Figure 1.The first few symmetric and anti-symmetric Hough modes for DW1 (s = 1, σ = 0.5) of scalar fields, computed using the normalized associated Legendre polynomial (ALP) expansions.Panels (a, b) are for symmetric modes; (c, d) are for anti-symmetric modes.The labels are [−1] for the first negative mode with largest negative eigenvalue, [+1] for the first positive mode with largest positive eigenvalue, and [0] for the so-called missing mode with zero eigenvalue or infinite equivalent depth.

Figure 2 .
Figure 2. The first few symmetric and anti-symmetric Hough modes for SW2 (s = 2, σ = 1), computed using the normalized ALP expansions.The left panels are symmetric modes, and the right panels are anti-symmetric modes, except panels (e, f) which are reversed.Panels (a, b) are for the scalar fields, (c, d) for the zonal wind component, and (e, f) for the meridional wind component.The labels are conventional.

Figure 3 .
Figure 3.The first few symmetric and anti-symmetric Hough modes for TW3 (s = 3, σ = 1.5) of scalar fields, computed using the normalized ALP expansions.The left panels are symmetric modes, and the right panels are anti-symmetric modes.

Figure 4 .
Figure 4.The absolute value of the expansion coefficients b n in Eq. (15), truncated at N = 150.The left panels are for the terdiurnal tides, s = 3, σ = 1.5, for eigenfunction with eigenvalue γ = 17.2098:(a) without parity factor, (b) with parity factor.The right panels are for pentadiurnal tides s = 5, σ = 2.5, for eigenfunction with eigenvalue γ = 22.9721: (c) without parity factor, (d) with parity factor.An empirical fitting curve is also shown in red dash.

Figure C1 .
Figure C1.Schematic isolines for Fourier terms a s (ϕ) cos(sλ) for various zonal wave numbers s, shown in a polar projection.Positivevalued isolines are solid black, while negative-valued isolines are red dashed.The thick yellow line segments depict a part of a meridian.For odd wave numbers (upper panels), the yellow lines connect solid black contours to red dashed isolines -the function changes sign along the meridian.

Figure C2 .
Figure C2.Schematic of the behavior of a s (ϕ) cos(sλ) along a meridian.If a s (0) = 0, the Fourier term will have a jump discontinuity across the pole (thick black curve) when longitude jumps by π .