Clenshaw algorithm explained

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

Clenshaw algorithm

In full generality, the Clenshaw algorithm computes the weighted sum of a finite series of functions

\phi _{k(x)} :$$S(x)\; =\; \backslash sum\_^n\; a\_k\; \backslash phi\_k(x)$$where \phi _{k,} k =0 , 1 , \ldots is a sequence of functions that satisfy the linear recurrence relation$$\backslash phi\_(x)\; =\; \backslash alpha\_k(x)\backslash ,\backslash phi\_k(x)\; +\; \backslash beta\_k(x)\backslash ,\backslash phi\_(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:

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)\; =\; \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 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)\; =\; \backslash sum\_^n\; a\_k\; x^k.$$The functions are simplyand 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\_(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)\; +\; \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\_(x)\; -\; b\_(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)\; =\; a\_0\; +\; 2xb\_1(x)\; -\; b\_2(x),$$(note the doubled a₀ coefficient) followed by$$p\_n(x)\; =\; \backslash tfrac\; \backslash left[a\_0+b\_0(x)\; -\; b\_2(x)\backslash right].$$

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(\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\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$$\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 byThe 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(\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

\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(\backslash theta\_1)-m(\backslash theta\_2)\; =\; C\_0(\backslash theta\_1-\backslash theta\_2)\; +\; \backslash sum\_^n\; 2\; C\_k\; \backslash sin\backslash bigl(k(\backslash theta\_1-\backslash theta\_2)\backslash bigr)\; \backslash cos\backslash bigl(k(\backslash theta\_1+\backslash theta\_2)\backslash bigr).$$Clenshaw summation can be applied in this case^[2] provided we simultaneously compute

m(\theta _{1)+m(\theta} _{2)} and perform a matrix summation,$$\backslash mathsf\; M(\backslash theta\_1,\backslash theta\_2)\; =\; \backslash begin\; (m(\backslash theta\_1)\; +\; m(\backslash theta\_2))\; /\; 2\backslash \backslash \; (m(\backslash theta\_1)\; -\; m(\backslash theta\_2))\; /\; (\backslash theta\_1\; -\; \backslash theta\_2)\; \backslash end\; =\; C\_0\; \backslash begin\; \backslash mu\; \backslash \backslash \; 1\; \backslash end\; +\; \backslash sum\_^n\; C\_k\; \backslash mathsf\; F\_k(\backslash theta\_1,\backslash theta\_2),$$whereThe first element of M(\theta _{1,\theta} _{2)} is the averagevalue of m and the second element is the average slope.F _{k(\theta} _{1,\theta} _{2)} satisfies the recurrencerelation$$\backslash mathsf\; F\_(\backslash theta\_1,\backslash theta\_2)\; =\; \backslash mathsf\; A(\backslash theta\_1,\backslash theta\_2)\; \backslash mathsf\; F\_k(\backslash theta\_1,\backslash theta\_2)\; -\; \backslash mathsf\; F\_(\backslash theta\_1,\backslash theta\_2),$$wheretakes the place of \alpha in the recurrence relation, and \beta =-1 .The standard Clenshaw algorithm can now be applied to yieldwhere B _{k} are 2×2 matrices. Finally we have$$\backslash frac\; =\; \backslash mathsf\; M\_2(\backslash theta\_1,\; \backslash 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 F _{1} and A , we take \lim _{\delta} (\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

Notes and References

1. Clenshaw . C. W.. A note on the summation of Chebyshev series. 10.1090/S0025-5718-1955-0071856-0. Mathematical Tables and Other Aids to Computation. 0025-5718. 9. 51. 118. July 1955 . 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)} .
2. Karney. C. F. F.. 2024. 68. The area of rhumb polygons. Stud. Geophys. Geod.. 10.1007/s11200-024-0709-z. free. Appendix B.