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 . 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 :
where is a sequence of functions that satisfy the linear recurrence relation
where the coefficients and are known in advance.
The algorithm is most useful when are functions that are complicated to compute directly, but and are particularly simple. In the most common applications, does not depend on , and is a constant that depends on neither nor .
To perform the summation for given series of coefficients , compute the values by the "reverse" recurrence formula:
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 . After computing and ,
the desired sum can be expressed in terms of them and the simplest functions and :
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
The functions are simply
\phi_0(x) &= 1, \\
\phi_k(x) &= x^k = x\phi_{k-1}(x)
\end{align}
and are produced by the recurrence coefficients and .
In this case, the recurrence formula to compute the sum is
and, in this case, the sum is simply
which is exactly the usual Horner's method.
=Special case for Chebyshev series=
Consider a truncated Chebyshev series
The coefficients in the recursion relation for the Chebyshev polynomials are
with the initial conditions
Thus, the recurrence is
and the final results are
An equivalent expression for the sum is given by
=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
Leaving off the initial term, the remainder is a summation of the appropriate form. There is no leading term because .
The recurrence relation for is
making the coefficients in the recursion relation
and the evaluation of the series is given by
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 , so the end of the recurrence is simply ; the term is added separately:
Note that the algorithm requires only the evaluation of two trigonometric quantities and .
=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
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
\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 is the average
value of and the second element is the average slope.
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 in the recurrence relation, and .
The standard Clenshaw algorithm can now be applied to yield
\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 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 and to simultaneously compute and the derivative , provided that, in evaluating and , we take .
See also
- Horner scheme to evaluate polynomials in monomial form
- De Casteljau's algorithm to evaluate polynomials in Bézier form