Values and Coefficients

The Two Faces of Polynomial Interpolants¶

Given $n+1$ data points $(x_j, f_j)$ where $x_j$ are Chebyshev points, the polynomial interpolant $p_n(x)$ is unique. But we can represent it in two equivalent ways:

Value Representation¶

Using Lagrange polynomials $\ell_j(x)$:

The polynomial is uniquely determined by the sampled function values $\{f_j\}$.

Coefficient Representation¶

Using Chebyshev polynomials $T_k(x)$:

The polynomial is uniquely determined by its Chebyshev coefficients $\{c_k\}$.

The Connection: A Linear Map¶

Evaluating the Chebyshev series at the data points gives:

This is a matrix-vector product:

Since $T_k(x_j) = T_k(\cos(j\pi/n)) = \cos(jk\pi/n)$, this is exactly the Discrete Cosine Transform matrix!

Two Representations¶

Consider a polynomial $p(x)$ of degree at most $n$.

Value Space¶

Store the values at the $n+1$ Chebyshev points:

where $x_k = \cos(k\pi/n)$ for $k = 0, 1, \ldots, n$.

Advantages:

• Direct sampling: just evaluate $f(x_k)$

• Differentiation via matrices: $p'(x_k) = (D\mathbf{f})_k$

• Nonlinear operations easy: $\sin(f)$ is just componentwise

Coefficient Space¶

Store the Chebyshev coefficients:

Advantages:

• Integration has a simple formula

• Coefficient decay reveals smoothness

• Truncation gives best polynomial approximation

The Discrete Cosine Transform¶

The key insight: under $x = \cos\theta$, Chebyshev polynomials become cosines:

At Chebyshev points $x_j = \cos(j\pi/n)$, we have $\theta_j = j\pi/n$, so:

This is exactly the DCT-I matrix!

Values to Coefficients¶

Given values $\mathbf{f}$, the coefficients are:

where $\sum''$ means the first and last terms are halved.

import scipy.fft as fft

def vals2coeffs(values):
    """Convert values at Chebyshev points to coefficients."""
    n = len(values) - 1
    if n == 0:
        return values.copy()

    # Use DCT-I
    coeffs = fft.dct(values[::-1], type=1, norm='forward')

    # Scale interior coefficients
    coeffs[1:n] *= 2.0
    return coeffs

Coefficients to Values¶

Given coefficients $\mathbf{c}$, the values are:

def coeffs2vals(coeffs):
    """Convert coefficients to values at Chebyshev points."""
    n = len(coeffs) - 1
    if n == 0:
        return coeffs.copy()

    # Undo scaling
    coeffs_scaled = coeffs.copy()
    coeffs_scaled[1:n] /= 2.0

    # Use inverse DCT-I
    values = fft.idct(coeffs_scaled, type=1, norm='forward')
    return values[::-1]

Example: $f(x) = x^3$¶

n = 4
x = np.cos(np.pi * np.arange(n) / (n-1))  # Chebyshev points
f = x**3                                    # Values

c = vals2coeffs(f)
# c = [0, 0.75, 0, 0.25]

The polynomial $x^3$ has the exact Chebyshev expansion:

Example: $f(x) = x^5$¶

Similarly:

The pattern: monomials have sparse Chebyshev representations (only odd or even terms appear).

The Clenshaw Algorithm¶

To evaluate $p(x) = \sum_{k=0}^{n} c_k T_k(x)$ without computing each $T_k$, use Clenshaw’s algorithm:

def clenshaw(x, coeffs):
    """Evaluate Chebyshev series at x using Clenshaw algorithm."""
    n = len(coeffs) - 1
    if n == 0:
        return coeffs[0] * np.ones_like(x)

    b_k1 = np.zeros_like(x)
    b_k2 = np.zeros_like(x)

    for k in range(n, 1, -1):
        b_k2, b_k1 = b_k1, coeffs[k] + 2*x*b_k1 - b_k2

    return coeffs[0] + x*b_k1 - b_k2

This is analogous to Horner’s method for monomial expansions.

Why Two Representations?¶

Different operations are natural in each space:

The Freedom to Choose¶

The key insight is that translation between representations is cheap ($O(n \log n)$), so we can work in whichever space is most convenient for each operation.

Example: Computing $\int_{-1}^{1} f(x)^2 \, dx$

Given the interpolant $p_n(x)$ for $f(x)$:

1. Square in value space: $(f_0^2, f_1^2, \ldots, f_n^2)$ — just componentwise!

2. Transform to coefficient space: Use DCT

3. Integrate in coefficient space: Use the closed-form Chebyshev integral formula

This hybrid approach is often optimal.

Extension to Infinite Dimensions¶

For a general Lipschitz continuous function $f(x)$, the Chebyshev series is infinite:

This “infinite vector” is called a quasivector in the literature. While working with infinitely many coefficients requires care (functional analysis!), the practical reality is that:

The Chebfun Philosophy¶

Modern software like Chebfun (MATLAB) and ApproxFun (Julia) represent functions as their Chebyshev coefficients, automatically:

1. Sampling the function at Chebyshev points

2. Computing coefficients via DCT

3. Adaptively choosing $n$ until coefficients decay to machine precision

This gives “numerical functions” that can be manipulated like symbolic functions but with guaranteed accuracy.

# Conceptual chebfun-style workflow
def chebfun(f, tol=1e-14):
    """Create adaptive Chebyshev approximation."""
    for n in [16, 32, 64, 128, 256, 512, 1024]:
        x = chebpts(n)
        vals = f(x)
        coeffs = vals2coeffs(vals)

        # Check if coefficients have decayed
        if np.max(np.abs(coeffs[-3:])) < tol * np.max(np.abs(coeffs)):
            return coeffs  # Converged!

    raise ValueError("Function too complex or non-smooth")

Complexity Summary¶

The DCT/FFT connection makes transforming between representations cheap, so we can use whichever is more convenient for each operation.

Aliasing: The Sampling Pitfall¶

When we sample a function at $N$ discrete points, we cannot distinguish high-frequency components from low-frequency ones. This is aliasing—the bane of spectral methods.

The Nyquist Limit¶

For polynomial interpolation at $N+1$ points: we can represent polynomials of degree at most $N$. Anything higher gets aliased to lower frequencies.

Aliasing in Action¶

Consider sampling $\cos(10\theta)$ at $N = 8$ points ($\theta_j = j\pi/8$):

N = 8
theta = np.pi * np.arange(N) / N
f_high = np.cos(10 * theta)  # Frequency 10 > N/2 = 4

# These samples are IDENTICAL to cos(6*theta)!
f_low = np.cos(6 * theta)
np.allclose(f_high, f_low)  # True!

The frequency-10 wave “folds back” and appears as frequency 6. We can’t tell them apart from samples alone.

Why Aliasing Matters¶

Aliasing causes silent errors:

1. No warning: The computed coefficients look reasonable

2. Wrong answer: High-frequency content corrupts low-frequency coefficients

3. Hard to detect: Everything “works” until you compare with a finer grid

The 2/3 Rule for Nonlinear Terms¶

When computing products like $u \cdot v$ in spectral methods, aliasing becomes critical.

The problem: If $u$ and $v$ each have frequencies up to $N/2$, then $u \cdot v$ has frequencies up to $N$—exceeding our resolution!

def dealias_product(u_hat, v_hat):
    """Compute u*v with 2/3 dealiasing."""
    N = len(u_hat)
    M = 3 * N // 2  # Padding factor

    # Zero-pad to finer grid
    u_pad = np.zeros(M, dtype=complex)
    v_pad = np.zeros(M, dtype=complex)
    u_pad[:N//2] = u_hat[:N//2]
    u_pad[-(N//2):] = u_hat[-(N//2):]
    v_pad[:N//2] = v_hat[:N//2]
    v_pad[-(N//2):] = v_hat[-(N//2):]

    # Multiply in physical space
    uv = np.fft.ifft(u_pad) * np.fft.ifft(v_pad)

    # Transform back and truncate
    uv_hat = np.fft.fft(uv)
    result = np.zeros(N, dtype=complex)
    result[:N//2] = uv_hat[:N//2]
    result[-(N//2):] = uv_hat[-(N//2):]
    return result * (M / N)  # Normalization

Practical Guidelines¶

Detecting Aliasing¶

How to know if your solution is aliased:

1. Coefficient plateau: If $|c_k|$ doesn’t decay to machine precision, you may need more points

2. Resolution test: Double $N$ and compare—if answer changes significantly, you were underresolved

3. Energy in high modes: If significant energy near $k = N/2$, aliasing is likely

def check_resolution(coeffs, tol=1e-10):
    """Warn if coefficients suggest underresolution."""
    if np.max(np.abs(coeffs[-5:])) > tol * np.max(np.abs(coeffs)):
        print("Warning: coefficients not fully resolved—aliasing possible!")

The Moral¶

“Aliasing is the price of discretization.”

When sampling continuous functions at discrete points, information is lost. Nyquist tells us exactly how much bandwidth we can capture. For spectral methods to work reliably:

• Linear problems: Ensure solution is smooth enough for the grid

• Nonlinear problems: Dealias religiously

• Unknown problems: Always check convergence