# A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions

Abstract: Spherical harmonic expansions form partial sums of fully normalised associated Legendre functions (ALFs). However, when evaluated increasingly close to the poles, the ultra-high degree and order (e.g. 2700) ALFs range over thousands of orders of magnitude. This causes existing recursion techniques for computing values of individual ALFs and their derivatives to fail. A common solution in geodesy is to evaluate these expansions using Clenshaw's method, which does not compute individual ALFs or their derivatives. Straightforward numerical principles govern the stability of this technique. Elementary algebra is employed to illustrate how these principles are implemented in Clenshaw's method. It is also demonstrated how existing recursion algorithms for computing ALFs and their first derivatives are easily modified to incorporate these same numerical principles. These modified recursions yield scaled ALFs and first derivatives, which can then be combined using Horner's scheme to compute partial sums, complete to degree and order 2700, for all latitudes (except at the poles for first derivatives). This exceeds any previously published result. Numerical tests suggest that this new approach is at least as precise and efficient as Clenshaw's method. However, the principal strength of the new techniques lies in their simplicity of formulation and implementation, since this quality should simplify the task of extending the approach to other uses, such as spherical harmonic analysis.

## Summary (5 min read)

### 1 Introduction

- Current geodetic practice is witnessing an increase in the construction and use of ultrahigh degree spherical harmonic expansions of the geopotential or topography.
- These have taken the form of simpleeffects models (eg. Featherstone, 1999; Nov´ak et al., 2001) for which synthetic geopotential coefficients up to degree and order2700 and2160, respectively, were produced without reference to a mass distribution.
- Hybrids of source and effects models also exist.
- Haagmans (2000) combines empirically determined coefficients with synthetic ones derived from numerical integration over isostatically compensated source masses to degree and order2160.
- The numerical means for including the necessary ALFs constitutes the principal chal- lenge to evaluating ultra-high degree spherical harmonic expansions.

### 1.1 Spherical Harmonic Expansions

- The general notationS (d), (d)m , X (d) m andP (d) nm( ) is used whenever a textual or mathematical reference applies to both the undifferentiated quantities and the first derivatives simultaneously.
- Underflows in the computation of anyP (d) nm( ) excludes the corresponding (matching degree and order) coefficients from contributing toS(d), whereas an overflow in the computation of one or moreP (d) nm( ) prevents any result forS (d) from being achieved at all.
- Thus, IEEE double precision only permits a maximum range of 620 orders of magnitude within which to compute and store the requiredP (d) nm( ) values.
- A numerically more stable alternative to these standard recursion relations is the Clenshaw (1955) method (cf. Tscherning and Poder, 1982; Gleason, 1985; Deakin, 1998).
- These principles are quite simple, both in concept and in application.

### 2 Forward Recursions for the Calculation of ALFs

- The most direct approach for evaluatingS (d) (Eq. 1) employs a recursive algorithm to computeP nm( ).
- These values ofP (d) nm( ) are multiplied by the correspondingEnm terms to yield the required series values ofX (d) m (Eq. 3), which subsequently yield (d)m (Eq. 2) and henceS(d) (Eq. 1).
- UnlikeP nm( ) cosm sinm , P (d6=0) nm ( ) cosm sinm is not normalised in 9 the correct usage of the word since its average squared value integrated over the unit sphere is not unity.
- This paper focuses solely on the computation of fully normalised ALFs and their first derivatives.
- No numerical tests were conducted for the above quasinormalisations.

### 2.1 Standard Forward Column Methods

- The most popular recursive algorithm used for computingP nm( ) in geodesy can be obtained by fully normalising, for example, Magnus et al. (1966, Eq. 4.3.3(2)).
- Note that in Fig. 2, the degree increases in rows down, the order increases in columns to the right, and the diagonal elements of the matrix are the sectoral values.
- This nomenclature will be employed throughout the paper.

### 2.2 Standard Forward Row Methods

- The next approach is termed thestandard forward row recursion (Fig. 3), and appears to be rarely used in geodesy.
- As with the standard forward column recursion (Section 2.1), the sectoralPmm( ) serve as seed values for the forward row recursion, and can be computed using Eq. (13).
- The standard forward row recursion computes nonsectoralP nm( ) of constantn (a ‘row’ in Fig. 3) and sequentially decreasingm (to the left (ie., ‘forward’) from the diagonal in Fig. 3).
- Using the same argument to that introduced for the forward column recursion, the nonexistant value ofP n;n+1( ) required to computeP n;n 1( ) in Eq. (18) may be disregarded because the corresponding recursion coefficient,hn;n 1, is always zero.
- Note that, to computeP nm( ) using the forward row recursion, Eq. (18) uses the cor- responding sectoral values of the samen, rather than the samem, as seed values.

### 2.3 Numerical Problems with the Standard Forward Methods

- Even when applied in IEEE double precision, both the standard forward column (Eq. 11) and standard forward row (Eq. 18) recursions will underflow forM > 1900 in the colatitude range 20Æ 160Æ.
- The numerical instability of both these forward recursions is noted in the geodetic literature (eg. Gleason, 1985) and elsewhere (eg. Libbrecht, 1985).
- 0 (ie., towards the poles) and asm increases.
- Accordingly, the high degree and order values ofPmm( ) will exceed the range of magnitudes capable being stored in IEEE double precision, thereby resulting in an underflow.

### 2.4 Other Normalisations and the Edmonds Recursion

- Belikov (1991) and Belikov and Taybatorov (1991) present a suite of recursive algorithms for computing the quantitieŝP (d)nm, whereP̂ (d) nm are related to un-normalisedP (d) nm according to the modified normalisation P̂nm( ) = 2 m n! (n+m)! Pnm( ) (25) However, this approach is also subject to numerical limitations.
- The description of the test 14 results which support this claim indicate that these computations were only performed at the equator, although the point is not clear.
- Nevertheless, Risbo (1996) includes a Fortran 77 subroutine for this recursion.
- Therefore, although it may be possible to incorporate the new approaches presented here into the Risbo (1996) approach, it is clear that the Edmonds recursion alone does not solve the numerical problems encountered towards the poles as described in Section 1.2.

### 2.5 The Modified Forward Row Method

- A simple, yet effective, method by which this problem of underflowingP mm( ) may be avoided is to eliminate theum term from the recursion process in Eq. (13).
- This will be demonstrated numerically in Section 4. 16 2.6 Modified Forward Column Method Values ofP (d) nm( ) um may also be computed using what will be termed thefirst modified forward column recursion.
- Thus, for the remainder of this paper, no partial sumsS(1) will be computed at the poles.

### 3.2 Reverse Column Method

- Results from timing tests presented in Section 4.3 show the reverse column methods, described below, to be highly inefficient in comparison with the other approaches presented in this paper for evaluating the required partial sums.
- Thus, the reverse column methods are used here to highlight the basic similarities and differences between these two approaches.
- To compute any value ofPnm( ) Pmm( ) , the second modified forward column recursion (eg. 34) aggregates the necessary(almt) andblm recursive terms in the sequence of increasing degreel (sequentially down each column in Fig. 2).
- The effect of using the recursion in Eq. (46) in this way is to sequentially aggregate the (almt) andblm recursive terms, in the sequence of decreasingl, until the recursion terminates at the computation ofsmm = Pnm( ) Pmm( ) .

### 3.3 Standard Clenshaw Methods

- The standard Clenshaw methods, summarised below, closely resemble the reverse column recursions (Section 3.2).
- The Clenshaw (1955) approach, which was formulated originally to evaluate partial sums of Chebyshev polynomials, was adapted for use in geodesy by Gulick (1970) to compute partial sums ofP (d) nm( ).
- Section 3.3.1 introduces a sim- ple implementation of the Clenshaw (1955) approach, whilst Section 3.3.2 presents the implementation that is used more commonly in geodesy (cf. Gleason, 1985).

### 3.3.1 The first Clenshaw method

- The simplest implementation of the Clenshaw (1955) technique, herein termed thefirst Clenshaw method, uses the recursions in Eqs. (46) and (47) to compute, directly, the intermediate sumsX (d) m Pmm( ) , without evaluating individual values ofP (d) nm( ) Pmm( ) .
- As in the reverse column recursion (Section 3.2), the recursive process terminates at the computation ofsmm (on the diagonal in Fig. 7), except that, in this case, the sectoralsmm = Xm Pmm( ) .
- This study is confined to the computation ofS (1).
- These seed values allow the recursive computation of all_slm , of constantm and sequentially decreasingl, from _snm to _smm .

### 3.3.2 The second Clenshaw method

- This allows the recursive computation of all slm , of constantm (a column in Fig. 7), and sequentially increasingl (upwards and towards the diagonal in Fig. 7), fromsnm to smm , wheresmm = Xm Pmm( )qm .
- The standard approach (cf. Gleason, 1985; Deakin, 1998) is to combine these values using the implementation of Horner’s scheme given in Eq. (50).

### 4.1 Viable Methods

- These algorithms are summarised in Table 1.
- As mentioned in Section 2.6, the first (MFC-1) and second (MFC-2) modified forward column recursions are, essentially, a single method.
- This leaves, to this point, five separate methods for computing spherical harmonic expansions.
- The tests of numerical accuracy will use analytic solutions for the sum of the square ofP (d) nm( ), to compare the modified forward row and modified forward column algorithms only.

### 4.2 Relative Numerical Precision

- The first step in comparing different methods for computing the partial sumsS (d) in Eq. (1) is to choose some appropriate values forEnm and .
- The relative precision (RP) for each was calculated using RP = s(d) s(d)s(d) (55) wheres(d) is the value of Eq. (54) for the summation computed in precision by each respective method ands(d) is the result for the same sum computed using the second Clenshaw method in precision.
- That is, all combinations of the five algorithms with the three implementations of Horner’s scheme showed that the choice of implementation of Horner’s scheme was irrelevant to the observed precision of the algorithm.
- One interesting feature is that the relative precision signatures obtained from all of the column methods for the northern latitudes are almost identical, whilst the signatures for the southern latitudes are not (cf. Figs. 8 through 10).
- This phenomenon is not observed in the accuracy tests for the two modified forward algorithms (Section 4.4).

### 4.3 Numerical Efficiency

- The five methods that successfully computedS (d), for M = 2700 and for integer values of to the poles, were tested for their relative numerical efficiency.
- The square roots and inverted square roots required to construct the recursion coefficients were computed once by each algorithm and then stored for multiple use in the synthesis 32 subroutines.
- The CPU times for the reverse column algorithm are excessively large and so were extrapolated from the computation times for a single parallel.
- All computations were performed, once again, on aSun Ultra 10 (333MHz) workstation that uses a virtual or ‘swapped’ RAM configuration, which is slower than actual RAM.
- It should be noted that these CPU times, in addition to showing the relative efficiency of each approach, are also functions of the computer architecture, compiler and programming language employed, as well as the programmer’s implementation of these algorithms.

### 4.4 Accuracy

- As mentioned, the numerical evaluations presented in Section 4.2 are tests of precision only, since both the tested algorithms and the ‘control’ algorithm computed in IEEE extended double precision may contain shared systematic errors.
- These could be due to any one of compiler, computer architecture or programming errors, for example.
- This necessitates that the squaring ofP (d) nm( ) um 10 280, as well as their combination using Horner’s scheme, be performed in IEEE extended double precision.
- The modified forward column algorithm becomes increasingly less accurate than the modified forward row algorithm as the computation approaches the poles.
- These results support those obtained for these two algorithms in the precision tests (Section 4.2).

### 5 Summary, Conclusion and Recommendation

- This paper has shown that standard Clenshaw methods for evaluating high degree spherical harmonic expansions derive their stability from simple numerical principles.
- An efficient technique for the computation of the gravimetric quantities from geoptential earth models.
- 547–555. Libbrecht KG (1985) Practical considerations for the generation of large order spherical harmonics, also known as Geophys J Int 139.
- The modified forward column and modified forward row recursions are immediately more versatile than the standard Clenshaw methods, since they do not automatically combine quantities of the same orderm.

Did you find this useful? Give us your feedback

...read more

##### Citations

271 citations

### Cites methods from "A unified approach to the Clenshaw ..."

...If starting values used in the recursion are appropriately scaled, as is summarized by Holmes and Featherstone (2002), these can be computed to high accuracy up to a maximum spherical harmonic degree of about 2700....

[...]

207 citations

### Cites background or methods from "A unified approach to the Clenshaw ..."

...Wenzel single normalization (1998), or Horner’s schemes (in cos φ) on partial sums such as in Holmes and Featherstone (2002),...

[...]

...Figure 3a shows that Plm(x) becomes nonsignificant above some maximum order, a phenomenon used by Jekeli et al. (2007)—see below. On Fig. 3b it is clear that the Hlm’s increase regularly with l before becoming significant, then oscillate around a stable value; and the Plm’s behave similarly (starting at ∼cos mφ), which we will use in our algorithm. When using recursive relations on the Hlm(x), as it has become customary in geodesy, usual tricks, e.g. Wenzel single normalization (1998), or Horner’s schemes (in cos φ) on partial sums such as in Holmes and Featherstone (2002),...

[...]

...Figure 3a shows that Plm(x) becomes nonsignificant above some maximum order, a phenomenon used by Jekeli et al. (2007)—see below....

[...]

...When using recursive relations on the Hlm(x), as it has become customary in geodesy, usual tricks, e.g. Wenzel single normalization (1998), or Horner’s schemes (in cos ϕ) on partial sums such as in Holmes and Featherstone (2002), do not work at such level....

[...]

111 citations

110 citations

### Cites background or methods from "A unified approach to the Clenshaw ..."

...Holmes and Featherstone [2002] noted that the standard forward column method, applied at double precision, will produce numerical underflow for m > 1900; we therefore expect that our algorithm will produce reliable output for orders considerably higher than m = 72....

[...]

...The recursion relations given in equations (A6), (A7), (A8), and (A12) are adapted fromHolmes and Featherstone [2002, equations (13), (11), (15), and (21)]. The recursion coefficients are largely the same as those given in Holmes and Featherstone [2002]; we have relabeled their f as d, and named a new coefficient c....

[...]

108 citations

##### References

5,050 citations

4,377 citations

### "A unified approach to the Clenshaw ..." refers background in this paper

...Risbo (1996) claims that the Edmonds (1957) recursion for D-matricies can be used to compute fully normalised ALFs up to degree 200,000....

[...]

3,100 citations

818 citations

### "A unified approach to the Clenshaw ..." refers methods in this paper

...For example, all methods will report underflows (8 ) when EGM96 coefficients are employed for the lower degrees....

[...]

...Lemoine FG, Kenyon SC, Factor JK, Trimmer RG, Pavlis NK, Chinn DS, Cox CM, Klosko SM, Luthcke SB, Torrence MH, Wang YM, Williamson RG, Pavlis EC, Rapp RH, Olson TR (1998) The development of the joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) geopotential model EGM96, NASA/TP-1998206861, National Aeronautics and Space Administration, Washington....

[...]

...For Enm , one might use empirically generated coefficients such as EGM96 (Lemoine et al., 1998) and/or GPM98B (Wenzel, 1998) to compute the lower degreeEnm ....

[...]

...For Enm , one might use empirically generated coefficients such as EGM96 (Lemoine et al., 1998) and/or GPM98B (Wenzel, 1998) to compute the lower degree Enm ....

[...]