Chebyshev polynomials: orthogonality (secret superpower!)
Chebyshev polynomials: integrals/derivative
exercise: approximating sin/cos

Chebyshev polynomials (½)

Chebyshev polynomials are defined on the interval 𝑥 ∈ [−1, +1]:

𝑇₀(𝑥) = 1     𝑇₁(𝑥) = 𝑥     𝑇ₙ(𝑥) = 2𝑥 𝑇ₙ₋₁(𝑥) − 𝑇ₙ₋₂(𝑥)

𝑇ₙ(𝑥) = 𝑐𝑜𝑠(𝑛 𝑎𝑐𝑜𝑠(𝑥))

These are orthogonal:

                    ________       / 0    if n != m
∫₋₁⁺¹ 𝑇ₙ(𝑥)𝑇ₘ(𝑥) / √ 1 − 𝑥²  𝑑𝑥 = {  π    if n = m = 0
                                   \ π/2  if n = m != 0

                          / 0   if i != j
∑ ₖ₌₀ᴺ⁻¹ 𝑇ᵢ(𝑥ₖ) 𝑇ⱼ(𝑥ₖ) = {  𝑁   if i = j = 0
                          \ 𝑁/2 if i = j != 0

𝑥ₖ = 𝑐𝑜𝑠((2𝑘 + 1) π / (2𝑁))     𝑘 = 0, 1, ... , 𝑁 − 1

This is your secret new superpower! Use it to approximate 𝑒𝑥𝑝(−𝑥²). How quickly do coefficients of the Chebyshev series vanish?

Chebyshev polynomials (2/2)

So far, we have used Chebyshev polynomials of the first kind, 𝑇ₙ(𝑥). There are also Chebyshev polynomials of the second kind, 𝑈ₙ(𝑥):

𝑇₀(𝑥) = 1     𝑇₁(𝑥) = 𝑥      𝑇ₙ(𝑥) = 2𝑥 𝑇ₙ₋₁(𝑥) − 𝑇ₙ₋₂(𝑥)
𝑈₀(𝑥) = 1     𝑈₁(𝑥) = 2𝑥     𝑈ₙ(𝑥) = 2𝑥 𝑈ₙ₋₁(𝑥) − 𝑈ₙ₋₂(𝑥)

𝑇ₙ(𝑥) = 𝑐𝑜𝑠(𝑛 𝑎𝑐𝑜𝑠(𝑥))       𝑈ₙ(𝑥) = 𝑠𝑖𝑛((𝑛 + 1) 𝑎𝑐𝑜𝑠(𝑥)) / 𝑠𝑖𝑛(𝑎𝑐𝑜𝑠(𝑥))

There are similar orthogonality relations for 𝑈ₙ(𝑥) - check a reference. Moreover, there are relations for derivatives and integrals:

                                     (𝑛 + 1) 𝑇ₙ₊₁(𝑥) − 𝑥𝑈ₙ(𝑥)) 
𝜕𝑇ₙ(𝑥)/𝜕𝑥 = 𝑛 𝑈ₙ₋₁(𝑥)    𝜕𝑈ₙ(𝑥)/𝜕𝑥 = -------------------------
                                              (𝑥² − 1)

            𝑇ₙ₊₁(𝑥)     𝑇ₙ₋₁(𝑥)
∫𝑇ₙ(𝑥)𝑑𝑥 = --------- - ---------    ∫𝑈ₙ(𝑥)𝑑𝑥 = 𝑇ₙ₊₁(𝑥) / (𝑛 + 1)
           2 (𝑛 + 1)   2 (𝑛 − 1)

Chebyshev polynomials for integration

              𝑇ₙ₊₁(𝑥)     𝑇ₙ₋₁(𝑥)
∫ 𝑇ₙ(𝑥) 𝑑𝑥 = --------- - ---------
             2 (𝑛 + 1)   2 (𝑛 − 1)

Use your Chebyshev approximation of 𝑒𝑥𝑝(−𝑥²) to derive the integral of the function.

When this is used for numerical integration of functions, it’s also known as Clenshaw-Curtis quadrature.

Chebyshev polynomials for integration: code

You can find the code in C++ in day5-ex1.cxx.

If you look at the coefficients, every other coefficient is zero because 𝑒𝑥𝑝(−𝑥²) is an even function. One easy way to deal with this is to exploit this symmetry, and use something like 𝑒𝑥𝑝(−(2.5 + 2.5𝑥)²) for 𝑥 > 0 instead.

For high precision work, it’s useful to calculate with higher numerical precision. One could use the long double data type in C/C++, or use an arbitary precision floating point library or calculator (like calc, see day5-ex1.calc).

approximating sin/cos

In many applications, you need both sin(x) and cos(x) with 𝑥 ∈ [−π, π], and you need it for many different values of x. In the C/C++ standard library, there is the sincos function which calculates both simultaneously. Pick one or more of the following problems:

approximate sin(x) and cos(x) with Chebyshev polynomials
approximate sin(x) with a Chebyshev polynomial on the interval [0, π], get cos(x) through trignonometric identities
approximate sin(x) and cos(x) through argument doubling on the unit circle
approximate sin(x) on the interval [0, π/2]. Can you turn this into a fast sin/cos implementation?

What’s the accuracy and how does speed compare?

approximating sin/cos: code (½)

Example code can be found in day5-ex2.cxx.

In the code, we aim for 𝒪(10⁻⁷) precision (although not all algorithms reach it). When compiled with g++ (or clang++), results similar to these can be obtained:

routine	error sin/cos	time/call	time/call
		(g++) [ns]	(clang++) [ns]
libc sin/cos	0/0	20.4	32.7
Cheb. for 𝑥∈[−π,π]	3.6e-7/3.0e-7	33.7	3.6
Cheb. for sin for 𝑥∈[0,π]	3.0e-7/1.9e-4	50.1	5.4
8-fold arg. doubling, Taylor	8.0e-5/5.6e-5	18.7	3.3
Cheb. for 𝑥∈[0,π/2]	9.5e-7/1.3e-6	29.8	3.7

approximating sin/cos: code (2/2)