Clenshaw algorithm

In numerical analysis, the Clenshaw algorithm, also called Clenshaw summation, 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.

In full generality, the Clenshaw algorithm computes the weighted sum of a finite series of functions ϕ k ( x ) :

S ( x ) = ∑ k = 0 n a k ϕ k ( x )

where ϕ k , k = 0 , 1 , … is a sequence of functions that satisfy the linear recurrence relation

ϕ k + 1 ( x ) = α k ( x ) ϕ k ( x ) + β k ( x ) ϕ k − 1 ( x ) ,

where the coefficients α k ( x ) and β k ( x ) are known in advance.

The algorithm is most useful when ϕ k ( x ) are functions that are complicated to compute directly, but α k ( x ) and β k ( x ) are particularly simple. In the most common applications, α ( x ) does not depend on k , and β is a constant that depends on neither x nor k .

To perform the summation for given series of coefficients a 0 , … , 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 + α k ( x ) b k + 1 ( x ) + β k + 1 ( x ) b k + 2 ( x ) .

Note that this computation makes no direct reference to the functions ϕ 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 ϕ 0 ( x ) and ϕ 1 ( x ) :

S ( x ) = ϕ 0 ( x ) a 0 + ϕ 1 ( x ) b 1 ( x ) + β 1 ( x ) ϕ 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 ) = ∑ k = 0 n a k x k .

The functions are simply

ϕ 0 ( x ) = 1 , ϕ k ( x ) = x k = x ϕ k − 1 ( x )

and are produced by the recurrence coefficients α ( x ) = x and β = 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 ) + ⋯ + a n T n ( x ) .

The coefficients in the recursion relation for the Chebyshev polynomials are

α ( x ) = 2 x , β = − 1 ,

with the initial conditions

T 0 ( x ) = 1 , T 1 ( x ) = x .

Thus, the recurrence is

b k ( x ) = a k + 2 x 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 ) = 2 a 0 + 2 x b 1 ( x ) − b 2 ( x ) ,

(note the doubled a₀ coefficient) followed by

p n ( x ) = 1 2 [ 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 ( θ ) = C 0 θ + C 1 sin ⁡ θ + C 2 sin ⁡ 2 θ + ⋯ + C n sin ⁡ n θ .

Leaving off the initial C 0 θ term, the remainder is a summation of the appropriate form. There is no leading term because ϕ 0 ( θ ) = sin ⁡ 0 θ = sin ⁡ 0 = 0 .

The recurrence relation for sin ⁡ k θ is

sin ⁡ ( k + 1 ) θ = 2 cos ⁡ θ sin ⁡ k θ − sin ⁡ ( k − 1 ) θ ,

making the coefficients in the recursion relation

α k ( θ ) = 2 cos ⁡ θ , β k = − 1.

and the evaluation of the series is given by

b n + 1 ( θ ) = b n + 2 ( θ ) = 0 , b k ( θ ) = C k + 2 cos ⁡ θ b k + 1 ( θ ) − b k + 2 ( θ ) , f o r n ≥ k ≥ 1.

The final step is made particularly simple because ϕ 0 ( θ ) = sin ⁡ 0 = 0 , so the end of the recurrence is simply b 1 ( θ ) sin ⁡ ( θ ) ; the C 0 θ term is added separately:

m ( θ ) = C 0 θ + b 1 ( θ ) sin ⁡ θ .

Note that the algorithm requires only the evaluation of two trigonometric quantities cos ⁡ θ and sin ⁡ θ .

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 ( θ 1 ) − m ( θ 2 ) = ∑ k = 1 n 2 C k sin ⁡ ( 1 2 k ( θ 1 − θ 2 ) ) cos ⁡ ( 1 2 k ( θ 1 + θ 2 ) ) .

Clenshaw summation can be applied in this case provided we simultaneously compute m ( θ 1 ) + m ( θ 2 ) and perform a matrix summation,

M ( θ 1 , θ 2 ) = [ ( m ( θ 1 ) + m ( θ 2 ) ) / 2 ( m ( θ 1 ) − m ( θ 2 ) ) / ( θ 1 − θ 2 ) ] = C 0 [ μ 1 ] + ∑ k = 1 n C k F k ( θ 1 , θ 2 ) ,

where

δ = 1 2 ( θ 1 − θ 2 ) , μ = 1 2 ( θ 1 + θ 2 ) , F k ( θ 1 , θ 2 ) = [ cos ⁡ k δ sin ⁡ k μ sin ⁡ k δ δ cos ⁡ k μ ] .

The first element of M ( θ 1 , θ 2 ) is the average value of m and the second element is the average slope. F k ( θ 1 , θ 2 ) satisfies the recurrence relation

F k + 1 ( θ 1 , θ 2 ) = A ( θ 1 , θ 2 ) F k ( θ 1 , θ 2 ) − F k − 1 ( θ 1 , θ 2 ) ,

where

A ( θ 1 , θ 2 ) = 2 [ cos ⁡ δ cos ⁡ μ − δ sin ⁡ δ sin ⁡ μ − sin ⁡ δ δ sin ⁡ μ cos ⁡ δ cos ⁡ μ ]

takes the place of α in the recurrence relation, and β = − 1 . The standard Clenshaw algorithm can now be applied to yield

B n + 1 = B n + 2 = 0 , B k = C k I + A B k + 1 − B k + 2 , f o r n ≥ k ≥ 1 , M ( θ 1 , θ 2 ) = C 0 [ μ 1 ] + B 1 F 1 ( θ 1 , θ 2 ) ,

where B k are 2×2 matrices. Finally we have

m ( θ 1 ) − m ( θ 2 ) θ 1 − θ 2 = M 2 ( θ 1 , θ 2 ) .

This technique can be used in the limit θ 2 = θ 1 = μ and δ = 0 to simultaneously compute m ( μ ) and the derivative d m ( μ ) / d μ , provided that, in evaluating F 1 and A , we take lim δ → 0 ( sin ⁡ δ ) / δ = 1 .

References

Clenshaw algorithm Wikipedia

(Text) CC BY-SA

Contents