Clenshaw algorithm#Geodetic applications

In numerical analysis, the Clenshaw algorithm, also called Clenshaw summation, is a recursive method to evaluate a linear combination of Chebyshev polynomials.{{Cite journal | last1 = Clenshaw | first1 = C. W.| title = A note on the summation of Chebyshev series| url = https://www.ams.org/journals/mcom/1955-09-051/S0025-5718-1955-0071856-0/| doi = 10.1090/S0025-5718-1955-0071856-0| journal = Mathematical Tables and Other Aids to Computation| issn = 0025-5718| volume = 9| issue = 51| page = 118| date=July 1955 | doi-access = free}} Note that this paper is written in terms of the Shifted Chebyshev polynomials of the first kind T^*_n(x) = T_n(2x-1). The method was published by Charles William Clenshaw in 1955. It is a generalization of Horner's method for evaluating a linear combination of monomials.

It generalizes to more than just Chebyshev polynomials; it applies to any class of functions that can be defined by a three-term recurrence relation.{{Citation |last1=Press |first1=WH |last2=Teukolsky |first2=SA |last3=Vetterling |first3=WT |last4=Flannery |first4=BP |year=2007 |title=Numerical Recipes: The Art of Scientific Computing |edition=3rd |publisher = Cambridge University Press |publication-place=New York |isbn=978-0-521-88068-8 |chapter=Section 5.4.2. Clenshaw's Recurrence Formula |chapter-url=http://apps.nrbook.com/empanel/index.html?pg=222}}

Clenshaw algorithm

In full generality, the Clenshaw algorithm computes the weighted sum of a finite series of functions \phi_k(x):

S(x) = \sum_{k=0}^n a_k \phi_k(x)

where \phi_k,\; k=0, 1, \ldots is a sequence of functions that satisfy the linear recurrence relation

\phi_{k+1}(x) = \alpha_k(x)\,\phi_k(x) + \beta_k(x)\,\phi_{k-1}(x),

where the coefficients \alpha_k(x) and \beta_k(x) are known in advance.

The algorithm is most useful when \phi_k(x) are functions that are complicated to compute directly, but \alpha_k(x) and \beta_k(x) are particularly simple. In the most common applications, \alpha(x) does not depend on k, and \beta is a constant that depends on neither x nor k.

To perform the summation for given series of coefficients a_0, \ldots, a_n, compute the values b_k(x) by the "reverse" recurrence formula:

\begin{align}

b_{n+1}(x) &= b_{n+2}(x) = 0, \\

b_k(x) &= a_k + \alpha_k(x)\,b_{k+1}(x) + \beta_{k+1}(x)\,b_{k+2}(x).

\end{align}

Note that this computation makes no direct reference to the functions \phi_k(x). After computing b_2(x) and b_1(x),

the desired sum can be expressed in terms of them and the simplest functions \phi_0(x) and \phi_1(x):

S(x) = \phi_0(x)\,a_0 + \phi_1(x)\,b_1(x) + \beta_1(x)\,\phi_0(x)\,b_2(x).

See Fox and Parker{{Citation |first1=Leslie |last1=Fox |first2=Ian B. |last2=Parker |title=Chebyshev Polynomials in Numerical Analysis |publisher=Oxford University Press |year=1968 |isbn=0-19-859614-6}} for more information and stability analyses.

Examples

=Horner as a special case of Clenshaw=

A particularly simple case occurs when evaluating a polynomial of the form

S(x) = \sum_{k=0}^n a_k x^k.

The functions are simply

\begin{align}

\phi_0(x) &= 1, \\

\phi_k(x) &= x^k = x\phi_{k-1}(x)

\end{align}

and are produced by the recurrence coefficients \alpha(x) = x and \beta = 0.

In this case, the recurrence formula to compute the sum is

b_k(x) = a_k + x b_{k+1}(x)

and, in this case, the sum is simply

S(x) = a_0 + x b_1(x) = b_0(x),

which is exactly the usual Horner's method.

=Special case for Chebyshev series=

Consider a truncated Chebyshev series

p_n(x) = a_0 + a_1 T_1(x) + a_2 T_2(x) + \cdots + a_n T_n(x).

The coefficients in the recursion relation for the Chebyshev polynomials are

\alpha(x) = 2x, \quad \beta = -1,

with the initial conditions

T_0(x) = 1, \quad T_1(x) = x.

Thus, the recurrence is

b_k(x) = a_k + 2xb_{k+1}(x) - b_{k+2}(x)

and the final results are

b_0(x) = a_0 + 2xb_1(x) - b_2(x),

p_n(x) = \tfrac{1}{2} \left[a_0+b_0(x) - b_2(x)\right].

An equivalent expression for the sum is given by

p_n(x) = a_0 + xb_1(x) - b_2(x).

=Meridian arc length on the ellipsoid=

Clenshaw summation is extensively used in geodetic applications.{{Citation

| last1=Tscherning

| first1=C. C.

| last2=Poder

| first2=K.

| year=1982

| title=Some Geodetic applications of Clenshaw Summation

| journal=Bolletino di Geodesia e Scienze Affini

| volume=41

| number=4

| pages=349–375

| url=http://cct.gfy.ku.dk/publ_cct/cct80.pdf

| access-date=2012-08-02

| archive-url=https://web.archive.org/web/20070612091533/http://cct.gfy.ku.dk/publ_cct/cct80.pdf

| archive-date=2007-06-12

| url-status=dead

}} A simple application is summing the trigonometric series to compute the meridian arc distance on the surface of an ellipsoid. These have the form

m(\theta) = C_0\,\theta + C_1\sin \theta + C_2\sin 2\theta + \cdots + C_n\sin n\theta.

Leaving off the initial C_0\,\theta term, the remainder is a summation of the appropriate form. There is no leading term because \phi_0(\theta) = \sin 0\theta = \sin 0 = 0.

The recurrence relation for \sin k\theta is

\sin (k+1)\theta = 2 \cos\theta \sin k\theta - \sin (k-1)\theta,

making the coefficients in the recursion relation

\alpha_k(\theta) = 2\cos\theta, \quad \beta_k = -1.

and the evaluation of the series is given by

\begin{align}

b_{n+1}(\theta) &= b_{n+2}(\theta) = 0, \\

b_k(\theta) &= C_k + 2\cos \theta \,b_{k+1}(\theta) - b_{k+2}(\theta),\quad\mathrm{for\ } n\ge k \ge 1.

\end{align}

The final step is made particularly simple because \phi_0(\theta) = \sin 0 = 0, so the end of the recurrence is simply b_1(\theta)\sin(\theta); the C_0\,\theta term is added separately:

m(\theta) = C_0\,\theta + b_1(\theta)\sin \theta.

Note that the algorithm requires only the evaluation of two trigonometric quantities \cos \theta and \sin \theta.

=Difference in meridian arc lengths=

Sometimes it necessary to compute the difference of two meridian arcs in a way that maintains high relative accuracy. This is accomplished by using trigonometric identities to write

m(\theta_1)-m(\theta_2) = C_0(\theta_1-\theta_2) + \sum_{k=1}^n 2 C_k

\sin\bigl({\textstyle\frac12}k(\theta_1-\theta_2)\bigr)

\cos\bigl({\textstyle\frac12}k(\theta_1+\theta_2)\bigr).

Clenshaw summation can be applied in this case{{ cite journal

| last = Karney

| first = C. F. F.

| year = 2024

| volume = 68

| number = 3-4

| pages = 99-120

| title = The area of rhumb polygons

| journal = Stud. Geophys. Geod.

| doi = 10.1007/s11200-024-0709-z

| doi-access = free

| postscript = Appendix B| arxiv = 2303.03219

}}

provided we simultaneously compute m(\theta_1)+m(\theta_2)

and perform a matrix summation,

\mathsf M(\theta_1,\theta_2) = \begin{bmatrix}

(m(\theta_1) + m(\theta_2)) / 2\\

(m(\theta_1) - m(\theta_2)) / (\theta_1 - \theta_2)

\end{bmatrix} =

C_0 \begin{bmatrix} \mu \\ 1 \end{bmatrix} +

\sum_{k=1}^n C_k \mathsf F_k(\theta_1,\theta_2),

where

\begin{align}

\delta &= \tfrac{1}{2}(\theta_1-\theta_2), \\[1ex]

\mu &= \tfrac{1}{2}(\theta_1+\theta_2), \\[1ex]

\mathsf F_k(\theta_1,\theta_2) &=

\begin{bmatrix}

\cos k \delta \sin k \mu \\

\dfrac{\sin k \delta}\delta \cos k \mu

\end{bmatrix}.

\end{align}

The first element of \mathsf M(\theta_1,\theta_2) is the average

value of m and the second element is the average slope.

\mathsf F_k(\theta_1,\theta_2) satisfies the recurrence

relation

\mathsf F_{k+1}(\theta_1,\theta_2) =

\mathsf A(\theta_1,\theta_2) \mathsf F_k(\theta_1,\theta_2) -

\mathsf F_{k-1}(\theta_1,\theta_2),

where

\mathsf A(\theta_1,\theta_2) = 2\begin{bmatrix}

\cos \delta \cos \mu & -\delta\sin \delta \sin \mu \\

- \displaystyle\frac{\sin \delta}\delta \sin \mu & \cos \delta \cos \mu

\end{bmatrix}

takes the place of \alpha in the recurrence relation, and \beta=-1.

The standard Clenshaw algorithm can now be applied to yield

\begin{align}

\mathsf B_{n+1} &= \mathsf B_{n+2} = \mathsf 0, \\[1ex]

\mathsf B_k &= C_k \mathsf I + \mathsf A \mathsf B_{k+1} -

\mathsf B_{k+2}, \qquad\mathrm{for\ } n\ge k \ge 1, \\[1ex]

\mathsf M(\theta_1,\theta_2) &=

C_0 \begin{bmatrix}\mu\\1\end{bmatrix} +

\mathsf B_1 \mathsf F_1(\theta_1,\theta_2),

\end{align}

where \mathsf B_k are 2×2 matrices. Finally we have

\frac{m(\theta_1) - m(\theta_2)}{\theta_1 - \theta_2} =

\mathsf M_2(\theta_1, \theta_2).

This technique can be used in the limit \theta_2 = \theta_1 = \mu and \delta = 0 to simultaneously compute m(\mu) and the derivative dm(\mu)/d\mu, provided that, in evaluating \mathsf F_1 and \mathsf A, we take \lim_{\delta \to 0} (\sin k \delta)/\delta = k.

See also

References

{{reflist|30em}}

{{DEFAULTSORT:Clenshaw Algorithm}}

Category:Numerical analysis