Clenshaw algorithm

In numerical analysis, the Clenshaw algorithm is a recursive method to evaluate a linear combination of Chebyshev polynomials. 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.

Clenshaw algorithm

Suppose that $\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}$ and $\beta _{k}$ are known in advance. Note that in the most common applications, $\alpha (x)$ does not depend on $k$, and $\beta$ is a constant that depends on neither $x$ nor $k$.

Our goal is to evaluate a weighted sum of these functions

$S(x)=\sum _{k=0}^{n}a_{k}\phi _{k}(x)$

Given the coefficients $a_{0},\ldots ,a_{n}$, compute the values $b_{k}(x)$ 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).$

The linear combination of the $\phi _{k}$ satisfies:

$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 for more information and stability analyses.

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

$\phi _{0}(x)=1,$

$\phi _{k}(x)=x^{k}=x\phi _{k-1}(x)$

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}+xb_{k+1}(x)$

and, in this case, the sum is simply

$S(x)=a_{0}+xb_{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 sum is

$p_{n}(x)=a_{0}+xb_{1}(x)-b_{2}(x).$

One way to evaluate this is to continue the recurrence one more step, and compute

$b_{0}(x)=2a_{0}+2xb_{1}(x)-b_{2}(x),$

(note the doubled a₀ coefficient) followed by

$p_{n}(x)=b_{0}(x)-xb_{1}(x)-a_{0}={\frac {1}{2}}\left[b_{0}(x)-b_{2}(x)\right].$

Geodetic applications

Clenshaw's algorithm is extensively used in geodetic applications where it is usually referred to as Clenshaw summation. A simple application is summing the trigonometric series to compute the meridian arc. 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

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

Failed to parse (syntax error): {\displaystyle b_k(\theta) = C_k + 2\cos \theta \,b_{k+1}(\theta) - b_{k+2}(\theta),\quad\text{for $n\ge k \ge 1$}. }

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$.

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})=\sum _{k=1}^{n}2C_{k}\sin {\bigl (}{\textstyle {\frac {1}{2}}}k(\theta _{1}-\theta _{2}){\bigr )}\cos {\bigl (}{\textstyle {\frac {1}{2}}}k(\theta _{1}+\theta _{2}){\bigr )}.$

Clenshaw summation can be applied in this case 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

${\mathsf {F}}_{k}(\theta _{1},\theta _{2})={\begin{bmatrix}\cos k\delta \sin k\mu \\{\frac {\sin k\delta }{\delta }}\cos k\mu \end{bmatrix}},$

$\delta =\textstyle {\frac {1}{2}}(\theta _{1}-\theta _{2})$, and $\mu =\textstyle {\frac {1}{2}}(\theta _{1}+\theta _{2})$. 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 \\-{\frac {\sin \delta }{\delta }}\sin \mu &\cos \delta \cos \mu \end{bmatrix}}.$

The standard Clenshaw algorithm can now be applied to yield

Failed to parse (unknown function "\begin{aligned}"): {\displaystyle \begin{aligned} \mathsf B_{n+1} &= \mathsf B_{n+2} = \mathsf 0, \\ \mathsf B_k &= C_k \mathsf I + \mathsf A \mathsf B_{k+1} - \mathsf B_{k+2}, \qquad\text{for $n\ge k \ge 1$},\\ \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{aligned} }

where ${\mathsf {B}}_{k}$ are $2\times 2$ matrices. Finally we have

${\frac {m(\theta _{1})-m(\theta _{2})}{\theta _{1}-\theta _{2}}}={\mathsf {M}}_{2}(\theta _{1},\theta _{2}).$

In the limit $\theta _{2}\rightarrow \theta _{1}$, this relation correctly evaluates the derivative $dm(\theta _{1})/d\theta _{1}$, providing that, in evaluating ${\mathsf {F}}_{k}$ and ${\mathsf {A}}$, we take $\lim _{\delta \rightarrow 0}(\sin k\delta )/\delta =k$.

See also

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