用Python进行数值分析 - 插值(2)

import numpy as np
import matplotlib.pyplot as plt

# 切比雪夫插值多项式专题
# numpy API 文档：
# https://numpy.org/doc/stable/reference/generated/numpy.polynomial.chebyshev.Chebyshev.fit.html
from numpy.polynomial.chebyshev import Chebyshev

# 切比雪夫多项式是基于牛顿差商多项式基础之上的(可以理解为牛顿差商的进化版)，它能有效解决插值中的龙格现象
# 其一般表达式为：T(n) = cos[n*arccos(x)]
# 什么是切比雪夫的根？--> 需要了解点切比雪夫理论
# (这个牛人有很多理论，不要搞迷糊了，关于多项式的可以参考：https://www.renrendoc.com/paper/196219602.html）
# 使用切比雪夫根作为基点的插值多项式就叫做切比雪夫插值多项式

# 第一类切比雪夫n阶多项式
# x取值范围: [-1, 1]
def chebyShevFunc(n, x):
    return np.cos(n*np.arccos(x))

# 第二类切比雪夫n阶多项式
# x取值范围: [-1, 1]
def chebyShevFunc2(n, x):
    """
    parameters:
    n: 从0开始,代表n阶切比雪夫多项式
    x: 任意标量实数
    return:
    第二类切比雪夫递推计算结果
    """
    if n == 0:
        return 1
    if n == 1:
        return x
    return 2*x*chebyShevFunc2(n-1, x) - chebyShevFunc2(n-2, x)


# 使用第一类切比雪夫n阶多项式拟合函数
def chebyShevPolynomialDemo1():
    # 通过已知点使用切比雪夫拟合多项式函数表达式
    # y = x**3 + 2
    xElements = np.array([-1, -0.5, 0, 0.5, 1])
    yElements = np.array([1, 1.875, 2, 2.125, 3])
    result1 = Chebyshev.fit(xElements, yElements, deg = 3, domain=None, full=False)
    print('general result: ', result1)
    # 得到表达式：2.0 + 0.75*T_1(x) - (5.93439169e-17)*T_2(x) + 0.25*T_3(x)
    result2 = Chebyshev.fit(xElements, yElements, deg = 3, domain=None, full=True)
    print('detailed result: ', result2)

    # 用图来画一画
    samplePoints = 500  # 样本点都选取500个
    xdots = [-1 + i * 2 / (samplePoints-1) for i in range(samplePoints)] 
    # 如果想区别显示出来，可以将上面的xdots稍作调整，例如，放大到 [-2, 2]区间上去
    # xdots = [-2 + i * 4 / (samplePoints-1) for i in range(samplePoints)]
    ydots = [(x**3 + 2) for x in xdots]
    fitYDots = [(2.0 + 0.75*chebyShevFunc(1, x) - (5.93439169e-17)*chebyShevFunc(2, x) + 0.25*chebyShevFunc(3, x)) for x in xdots]
    plt.plot(xdots, ydots, label='actual function: f(x)=x**3+2', color='red')
    plt.plot(xdots, fitYDots, label='chebyshev fitting function', color='blue')
    plt.title(f'fitting function and actual function demo1')
    plt.legend()
    plt.show()
    # 通过图像，可以看到在该区间完美拟合

# 使用第一类切比雪夫n阶多项式拟合函数
def chebyShevPolynomialDemo2():
    # 通过已知点使用切比雪夫拟合多项式函数表达式
    # y = x**2 - 2
    xElements = np.array([-1, -0.5, 0, 0.5, 1])
    yElements = np.array([-1, -1.75, -2, -1.75, -1])
    result1 = Chebyshev.fit(xElements, yElements, deg=2, domain=None, full=False)
    print('general result: ', result1)
    # 得到表达式：-1.5 + (2.10650008e-16)*T_1(x) + 0.5 * T_2(x)
    result2 = Chebyshev.fit(xElements, yElements, deg=2, domain=None, full=True)
    print('detailed result: ', result2)
    
    # 用图来画一画
    samplePoints = 500  # 样本点都选取500个
    xdots = [-1 + i * 2 / (samplePoints-1) for i in range(samplePoints)]
    ydots = [(x**2 - 2) for x in xdots]
    fitYDots = [(-1.5 + (2.10650008e-16)*chebyShevFunc(1, x) + 0.5*chebyShevFunc(2, x)) for x in xdots]
    plt.plot(xdots, ydots, label='actual function: f(x)=x**2-2', color='red')
    plt.plot(xdots, fitYDots, label='chebyshev fitting function', color='blue')
    plt.title(f'fitting function and actual function demo2')
    plt.legend()
    plt.show()
    # 通过图像，可以看到在该区间完美


# 在更大的区间范围[-2, 2]上如何选取切比雪夫插值点去拟合呢？
# 选择方法：假如给定区间[a, b], 切比雪夫插值点为：(a+b)/2 + ((b-a)/2)*pi/2*n
def chebyShevPolynomialDemo3():
    # 通过已知点使用切比雪夫拟合多项式函数表达式
    # y = x**2 - 2
    xElements = np.array([(0.2)*np.pi, (0.6)*np.pi, np.pi, 1.4*np.pi, 1.8*np.pi])
    #计算系数的过程按照牛顿差商的法则
    coeff = [((0.2)*np.pi)**2-2]
    coeff.append(0.8*np.pi)
    coeff.append(1)
    coeff.append(0) # 意味着可以省略掉后面的多项式
    coeff.append(0) # 意味着可以省略掉后面的多项式
    #意味着表达式为：[(np.pi/5)**2-2] + (0.8*np.pi)*(x-np.pi/5) + (x-np.pi/5)(x-0.6*np.pi)
    
    # 用图来画一画
    samplePoints = 500  # 样本点都选取500个
    xdots = [-2 + i * 4 / (samplePoints-1) for i in range(samplePoints)]
    ydots = [(x**2 - 2) for x in xdots]
    fitYDots = [((np.pi/5)**2-2 + (0.8*np.pi)*(x-np.pi/5) + (x-np.pi/5)*(x-0.6*np.pi)) for x in xdots]
    plt.plot(xdots, ydots, label='actual function: f(x)=x**2-2', color='red')
    plt.plot(xdots, fitYDots, label='chebyshev fitting function', color='blue')
    plt.title(f'fitting function and actual function demo3')
    plt.legend()
    plt.show()
    # 通过图像，可以看到在该区间完美


if __name__ == "__main__":
    # chebyShevPolynomialDemo1()
    # chebyShevPolynomialDemo2()
    chebyShevPolynomialDemo3()