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

Abstract

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 (eg. 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 (1955) method, which does not compute individual ALFs or their derivatives. Straightforward numerical principles govern the stability of this technique. This paper employs elementary algebra to illustrate how these principles are implemented in Clenshaw's method. It also demonstrates 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.

High-frequency analysis of Earth gravity field models based on terrestrial gravity and GPS/levelling data: a case study in Greece

“…The numerical stability and accuracy of the normalized associated Legendre functions presented in Eq. 1 is a critical computing problem (Bosch 2000, Claessens 2005, Montenbruck and Gill 2000, especially in the case of very high degree and order (Wittwer et al 2008, Holmes andFeatherstone 2002). The numerical stability is treated here by implementing appropriate recurrence relations into our source code that are recommended for this purpose (Montenbruck and Gill 2000).…”

Discrete Spherical Harmonic Transforms: Numerical Preconditioning and Optimization

“…This strategy of biasing the EXPONENT is independent of the colatitude θ which is very important in this context. The literature on computing very high degree spherical harmonics contains numerous strategies such as using the Clenshaw summation approach [13], preconditioning the variable with 10 280 sinθ in SHTools [19], modifying the recursion to avoid high powers of sinθ [14], and others. These modifications have been reported to achieve degrees around 2700-2800 in synthesis computations [13,14] and in synthesis and analysis [19].…”

“…The literature on computing very high degree spherical harmonics contains numerous strategies such as using the Clenshaw summation approach [13], preconditioning the variable with 10 280 sinθ in SHTools [19], modifying the recursion to avoid high powers of sinθ [14], and others. These modifications have been reported to achieve degrees around 2700-2800 in synthesis computations [13,14] and in synthesis and analysis [19]. Notice that in synthesis only computations, it is only necessary to avoid underflows and set the corresponding incremental contributions to zero while in synthesis and analysis computations, the numerical recovery of the input spectral coefficients is expected in addition to avoiding possible underflows.…”

“…With proper normalization such as the geodetic one used in the following computations, one can achieve degrees and orders to around 1800 in REAL*8 or double precision arithmetic [4] and over 3600 in REAL*16 or quadruple precision arithmetic [5,6]. With latitudebased preconditioning, [13,14,19] achieved about 2700. With Clenshaw's approach, similar results can be achieved for appropriately decreasing spectral coefficients such as with the EGM06 model coefficients of degree and order 2190 [7,13].…”

“…With latitudebased preconditioning, [13,14,19] achieved about 2700. With Clenshaw's approach, similar results can be achieved for appropriately decreasing spectral coefficients such as with the EGM06 model coefficients of degree and order 2190 [7,13].…”