使用 python 實現 chebyshev 多項式及內插法

不囉唆, 詳如以下代碼：

import numpy as np
def Cp(x, n): # 快速疊代法, 計算 1st kind Chebyshev 多項式
  if (n == 0):
      return 1
  if (n == 1):
      return x
  pk_1 = 1
  pn_1 = x
  k = 1
  while k < n :
    pk = pn_1 * x * 2  - pk_1
    pk_1 = pn_1
    pn_1 = pk
    k += 1      
  return pn_1

test_f = lambda x: np.exp(-x*x) # 測試函式
Ln    = 20
px    = [0.0] * Ln
py    = [0.0] * Ln
theta = [0.0] * Ln
for k in range(Ln) :
    theta[k] = np.pi * (k + 0.5) / Ln;# θₖ = np.pi * (k + 0.5) / Ln
    px[k] = np.cos(theta[k]);# pxₖ = cos(θₖ)
    py[k] = test_f(px[k]);# pyₖ = f(pxₖ) to be used in Chebyshev Interpolation  

def Ci(i): # Coefficient, bind with pyₖ, θₖ, Ln
  sum = 0
  for k in range(Ln) :
    sum += py[k] * np.cos(i * theta[k])
  return sum * 2 / Ln # 2/n * Σₙf(pxₖ)*Tₖ(pxₖ), k = 0, 2, ... n - 1

def Chebyshev_interpolation(t): # Chebyshev Interpolation Polynomials
    sum = Ci(0) / 2;# Σₙ Cₖ*Tₖ(x) - C₀/2 = Σₖ Cₖ*Tₖ(x) + C₀/2, k = 1, 2, ... n-1, C₀*T₀(x) = C₀
    k = 1
    while k < Ln :
        sum += Ci(k) * Cp(t, k)
        k += 1
    return sum

for k in range(Ln) :
    xk = px[k] + 0.1
    c  = test_f(xk) # 實際值
    h  = Chebyshev_interpolation(xk)# Chebyshev 合成值
    print(f"x={xk}, 內插={h}, 實際值={c}, 誤差 = {h-c}")