チェビシェフ補間で多項式近似式を求める

cheby.py

# 参考: 『インターフェース』 CQ 出版, 2017 年 9 月号, p.109-111, 三上直樹
# ある函数 f(x) が下のような N 次の多項式 g(x) で表現できると假定し、そこそこ精度のある近似式を求める。
# f(x) ≈ g(x)
# g(x) ≈ c[0] * T[0] + c[1] * T[1] + ... + c[N] * T[N]
# T[N] はチェビシェフ多項式である。
# ここでは例として f(x) = Sin(Pi * x / 2) の近似式を具体的に計算してみる。区間は -1 ≦ x ≦ 1 とする。
# 区間が a ≦ x ≦ b である場合は次式 u(a, b) で変数を変換する。
# u = lambda a, b: (2 * x - (b + a)) / (b - a) 
import sympy as sym
import numpy as np

def cheby(f, N):
    x = sym.Symbol("x")

    k = np.arange(N+1)

    s = [2 / (N + 1), 2 * k + 1, np.pi / (2 * (N + 1)), (1 / (N + 1))]
    p = s[1] * s[2]    
    
    # 第 1 種チェビシェフ多項式を求める。
    T = [sym.chebyshevt_poly(n, x) for n in k]
    
    # 数列 r を求める。
    r = np.cos(p)
        
    # 数列 r を使って係数 c[0] を求める。
    fVals = f(r)
    c     = np.zeros(N+1)
    c[0]  = s[3] * np.sum(fVals)
        
    # 係数 c[1] ～ c[N] を求める。
    c[1:N+1] = [s[0] * np.sum(fVals * np.cos(p * n)) for n in np.arange(1, N+1)]

    # 係数 c とチェビシェフ多項式 T とで近似式 g(x) を組み立てる。
    g = np.sum([c[n] * T[n] for n in k])

    # 0 次 ～ N 次の係数を取り出す。
    return np.array([round(g.coeff(x, n), 8) for n in k])

########
# test #
########
if __name__ == "__main__":

    f = lambda x : np.sin(np.pi * x / 2) # この函数の近似式をチェビシェフ補間による多項式で求める。
    N = 5                                # N 次まで近似してみる。
    
    cApprx = cheby(f, 5)
    print(cApprx)