Clenshaw algorithm

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

Clenshaw algorithm

In full generality, the Clenshaw algorithm computes the weighted sum of a finite series of functions ${\varphi }_k(x)$: $$S(x)=\sum _{k=0}^n{a}_k{\varphi }_k(x)$$ where ${\varphi }_k,k=0,1,\dots$ is a sequence of functions that satisfy the linear recurrence relation $${\varphi }_{k+1}(x)={\alpha }_k(x){\varphi }_k(x)+{\beta }_k(x){\varphi }_{k-1}(x),$$ where the coefficients ${\alpha }_k(x)$ and ${\beta }_k(x)$ are known in advance.

The algorithm is most useful when ${\varphi }_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,\dots ,a_n$, compute the values $b_k(x)$ by the "reverse" recurrence formula: $$\begin{array}{cc}{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{array}$$

Note that this computation makes no direct reference to the functions ${\varphi }_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 ${\varphi }_0(x)$ and ${\varphi }_1(x)$: $$S(x)={\varphi }_0(x)a_0+{\varphi }_1(x)b_1(x)+{\beta }_1(x){\varphi }_0(x)b_2(x).$$

See Fox and Parker 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{array}{cc}{\varphi }_0(x)& =1,\\ {\varphi }_k(x)& =x^k=x{\varphi }_{k-1}(x)\end{array}$$ 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,\beta =-1,$$ with the initial conditions $$T_0(x)=1,T_1(x)=x.$$

Thus, the recurrence is $$b_k(x)=a_k+2x{b}_{k+1}(x)-b_{k+2}(x)$$ and the final sum is $$p_n(x)=a_0+x{b}_1(x)-b_2(x).$$

One way to evaluate this is to continue the recurrence one more step, and compute $$b_0(x)=a_0+2x{b}_1(x)-b_2(x),$$ (note the doubled a₀ coefficient) followed by $$p_n(x)=\frac{1}{2}[a_0+b_0(x)-b_2(x)].$$

Meridian arc length on the ellipsoid

Clenshaw summation is extensively used in geodetic applications. 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 ${\varphi }_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 ,{\beta }_k=-1.$$ and the evaluation of the series is given by $$\begin{array}{cc}{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 ),\mathrm{f}\mathrm{o}\mathrm{r}n\ge k\ge 1.\end{array}$$ The final step is made particularly simple because ${\varphi }_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}2C_k\sin (\frac{1}{2}k({\theta }_1-{\theta }_2))\cos (\frac{1}{2}k({\theta }_1+{\theta }_2)).$$ Clenshaw summation can be applied in this case provided we simultaneously compute $m({\theta }_1)+m({\theta }_2)$ and perform a matrix summation, $$𝖬({\theta }_1,{\theta }_2)=\left[\begin{array}{c}(m({\theta }_1)+m({\theta }_2))/2\\ (m({\theta }_1)-m({\theta }_2))/({\theta }_1-{\theta }_2)\end{array}\right]=C_0\left[\begin{array}{c}\mu \\ 1\end{array}\right]+\sum _{k=1}^n{C}_k{𝖥}_k({\theta }_1,{\theta }_2),$$ where $$\begin{array}{cc}\delta & =\frac{1}{2}({\theta }_1-{\theta }_2),\\ \mu & =\frac{1}{2}({\theta }_1+{\theta }_2),\\ {𝖥}_k({\theta }_1,{\theta }_2)& =\left[\begin{array}{c}\cos k\delta \sin k\mu \\ \frac{\sin k\delta }{\delta }\cos k\mu \end{array}\right].\end{array}$$ The first element of $𝖬({\theta }_1,{\theta }_2)$ is the average value of $m$ and the second element is the average slope. ${𝖥}_k({\theta }_1,{\theta }_2)$ satisfies the recurrence relation $${𝖥}_{k+1}({\theta }_1,{\theta }_2)=𝖠({\theta }_1,{\theta }_2){𝖥}_k({\theta }_1,{\theta }_2)-{𝖥}_{k-1}({\theta }_1,{\theta }_2),$$ where $$𝖠({\theta }_1,{\theta }_2)=2\left[\begin{array}{cc}\cos \delta \cos \mu & -\delta \sin \delta \sin \mu \\ -\frac{\sin \delta }{\delta }\sin \mu & \cos \delta \cos \mu \end{array}\right]$$ takes the place of $\alpha$ in the recurrence relation, and $\beta =-1$. The standard Clenshaw algorithm can now be applied to yield $$\begin{array}{cc}{𝖡}_{n+1}& ={𝖡}_{n+2}=\mathrm{𝟢},\\ {𝖡}_k& =C_k𝖨+𝖠{𝖡}_{k+1}-{𝖡}_{k+2},\mathrm{f}\mathrm{o}\mathrm{r}n\ge k\ge 1,\\ 𝖬({\theta }_1,{\theta }_2)& =C_0\left[\begin{array}{c}\mu \\ 1\end{array}\right]+{𝖡}_1{𝖥}_1({\theta }_1,{\theta }_2),\end{array}$$ where ${𝖡}_k$ are 2×2 matrices. Finally we have $$\frac{m({\theta }_1)-m({\theta }_2)}{{\theta }_1-{\theta }_2}={𝖬}_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 ${𝖥}_1$ and $𝖠$, we take $\underset{\delta \to 0}{\lim }(\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