orthogonal polynomial interpolation

Concept

basically in order to use orthogonal polynomial for interpolation, we setup the orthogonal polynomial coefficient matrix with selected polynomials and solve it for coefficients.

$$A=\left[\begin{array}{cccc}{P}_0(x_0)& P_1(x_0)& \cdots & P_{n-1}(x_0)\\ P_0(x_1)& P_1(x_1)& \cdots & P_{n-1}(x_1)\\ ⋮& ⋮& \ddots & ⋮\\ P_0(x_{n-1})& P_1(x_{n-1})& \cdots & P_{n-1}(x_{n-1})\end{array}\right]$$ $$A\vec{c}=\vec{b}$$

Or we can use

$$p_{n-1}(x)=\sum _{j=1}{c}_j{P}_j(x)$$

where $P$ is the selected polynomial. Then

$$c_j=\frac{(f,P_j{)}_w}{(P_j,P_j{)}_w}=\frac{{\int }_{-1}^{1}w(x)f(x)P_j(x)dx}{{\int }_{-1}^{1}w(x)P_j^2(x)dx}=\frac{{\int }_{-1}^{1}w(x)f(x)P_j(x)dx}{||P_j|{|}^2}$$

where $(f,P_j{)}_w$ is weighted inner-product.

interpolation error

$$f-p=\frac{f^{(n)}(\xi )}{n!}q(x)$$

where

$$q(x)=(x-x_1)(x-x_2)\dots (x-x_n)$$

Crude Error

assume $t_1<t_2<\cdots <t_n$ we have

$$\max |f(t)-p_{n-1}(t)|\le \frac{M{h}^n}{4n}$$

where $h={\max }_i\Delta t_i$ and $M$ is constant with $|f^{(n)}|\le M$.

Code

import numpy as np
 
def T(j, x):
    """Compute Chebyshev polynomial T_j(x) = cos(j * arccos(x))"""
    return np.cos(j * np.arccos(x))
 
def chebyshev_interp(f, n, t):
    """
    Use Chebyshev polynomials T_j(x) and Gauss-Chebyshev quadrature
    to interpolate function f on interval [-1, 1].
    
    f: function to interpolate
    n: number of basis functions
    t: array of evaluation points
    """
    # chebyshev_nodes
    x_nodes = np.empty(n)
    for i in range(n):
        x_nodes[i] = np.cos((2 * i + 1) * np.pi / (2 * n))
    
    f_vals = np.empty(n)
    for i in range(n):
        f_vals[i] = f(x_nodes[i])
    
    coeffs = np.empty(n)
    for j in range(n):
        num = 0
        denom = 0
        for i in range(n):
            Tij = T(j, x_nodes[i])
            num += f_vals[i] * Tij
            denom += Tij * Tij
        num *= np.pi / n
        denom *= np.pi / n
        coeffs[j] = num / denom
    
    res = np.empty_like(t)
    for i in range(len(t)):
        s = 0
        for j in range(n):
            s += coeffs[j] * T(j, t[i])
        res[i] = s
    
    return res