Clenshaw algorithm#Geodetic applications

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

$$S(x)\; =\; \backslash sum\_\{k=0\}^n\; a\_k\; \backslash phi\_k(x)$$

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

$$\backslash phi\_\{k+1\}(x)\; =\; \backslash alpha\_k(x)\backslash ,\backslash phi\_k(x)\; +\; \backslash beta\_k(x)\backslash ,\backslash phi\_\{k-1\}(x),$$

where the coefficients $\backslash alpha\_k(x)$ and $\backslash beta\_k(x)$ are known in advance.

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

To perform the summation for given series of coefficients $a\_0,\; \backslash ldots,\; a\_n$, compute the values $b\_k(x)$ by the "reverse" recurrence formula:

$$\backslash 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 $\backslash 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 $\backslash phi\_0(x)$ and $\backslash phi\_1(x)$:

$$S(x)\; =\; \backslash phi\_0(x)\backslash ,a\_0\; +\; \backslash phi\_1(x)\backslash ,b\_1(x)\; +\; \backslash beta\_1(x)\backslash ,\backslash phi\_0(x)\backslash ,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.

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

$$S(x)\; =\; \backslash sum\_\{k=0\}^n\; a\_k\; x^k.$$

The functions are simply

$$\backslash 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 $\backslash alpha(x)\; =\; x$ and $\backslash 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.

Consider a truncated Chebyshev series

$$p\_n(x)\; =\; a\_0\; +\; a\_1\; T\_1(x)\; +\; a\_2\; T\_2(x)\; +\; \backslash cdots\; +\; a\_n\; T\_n(x).$$

The coefficients in the recursion relation for the Chebyshev polynomials are

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

with the initial conditions

$$T\_0(x)\; =\; 1,\; \backslash 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)\; =\; \backslash tfrac\{1\}\{2\}\; \backslash left[a\_0+b\_0(x)\; -\; b\_2(x)\backslash right].$$

An equivalent expression for the sum is given by

$$p\_n(x)\; =\; a\_0\; +\; xb\_1(x)\; -\; b\_2(x).$$

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(\backslash theta)\; =\; C\_0\backslash ,\backslash theta\; +\; C\_1\backslash sin\; \backslash theta\; +\; C\_2\backslash sin\; 2\backslash theta\; +\; \backslash cdots\; +\; C\_n\backslash sin\; n\backslash theta.$$

Leaving off the initial $C\_0\backslash ,\backslash theta$ term, the remainder is a summation of the appropriate form. There is no leading term because $\backslash phi\_0(\backslash theta)\; =\; \backslash sin\; 0\backslash theta\; =\; \backslash sin\; 0\; =\; 0$.

The recurrence relation for $\backslash sin\; k\backslash theta$ is

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

making the coefficients in the recursion relation

$$\backslash alpha\_k(\backslash theta)\; =\; 2\backslash cos\backslash theta,\; \backslash quad\; \backslash beta\_k\; =\; -1.$$

and the evaluation of the series is given by

$$\backslash 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 $\backslash phi\_0(\backslash theta)\; =\; \backslash sin\; 0\; =\; 0$, so the end of the recurrence is simply $b\_1(\backslash theta)\backslash sin(\backslash theta)$; the $C\_0\backslash ,\backslash theta$ term is added separately:

$$m(\backslash theta)\; =\; C\_0\backslash ,\backslash theta\; +\; b\_1(\backslash theta)\backslash sin\; \backslash theta.$$

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

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(\backslash theta\_1)+m(\backslash 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

$$\backslash 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 $\backslash mathsf\; M(\backslash theta\_1,\backslash theta\_2)$ is the average

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

$\backslash mathsf\; F\_k(\backslash theta\_1,\backslash 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 $\backslash alpha$ in the recurrence relation, and $\backslash beta=-1$.

The standard Clenshaw algorithm can now be applied to yield

$$\backslash 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 $\backslash 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 $\backslash theta\_2\; =\; \backslash theta\_1\; =\; \backslash mu$ and $\backslash delta\; =\; 0$ to simultaneously compute $m(\backslash mu)$ and the derivative $dm(\backslash mu)/d\backslash mu$, provided that, in evaluating $\backslash mathsf\; F\_1$ and $\backslash mathsf\; A$, we take $\backslash lim\_\{\backslash delta\; \backslash to\; 0\}\; (\backslash sin\; k\; \backslash delta)/\backslash delta\; =\; k$.

• Horner scheme to evaluate polynomials in monomial form
• De Casteljau's algorithm to evaluate polynomials in Bézier form

{{reflist|30em}}

{{DEFAULTSORT:Clenshaw Algorithm}}

Category:Numerical analysis