切比雪夫插值介绍与应用

  在文章拉格朗日插值多项式的原理介绍及其应用中，笔者介绍了拉格朗日插值多项式及其应用，在文章插值多项式的龙格现象的介绍与模拟，我们发现插值多项式会在端点附近发生较大的扭动，即龙格现象。而切比雪夫插值是一种特定最优的点间距选取方式，其可以尽量减少插值误差。

插值误差

  在数值分析中，有如下定理：
  (定理1 插值误差公式)假设P(x)是n-1或者更低阶的插值多项式，其拟合n个点$(x_1,y_1),\cdots ,(x_n,y_n)$，则插值误差是：

$$f(x)-P(x)=\frac{1}{n!}(x-x_1)(x-x_2)\cdots (x-x_n)f^{(n)}(c)$$

其中c在最小和最大的n+1个数字$x,x_1,\cdots ,x_n$之间。

  从中我们可以发现，如果n确定，$x_1,\cdots ,x_n$是自由选取的n个点，那么插值误差的大小将取决于$|(x-x_1)(x-x_2)\cdots (x-x_n)|$和$f^{(n)}(c)$，而$f^{(n)}(c)$视函数f(x)而定，因此我们重点讨论$|(x-x_1)(x-x_2)\cdots (x-x_n)|$，我们将x的取值范围限定在区间[-1, 1]，有如下定理：

  （定理2）选择实数$-1\le x_1,\cdots ,x_n\le 1$，使得$\underset{-1\le x\le 1}{\max }|(x-x_1)\cdots (x-x_n)|$尽可能小，则$x_i=\cos \frac{(2i-1)\pi }{2n},i=1,\cdots ,n$，对应的最小值是$1/2^{n-1}$，实际上，通过$(x-x_1)\cdots (x-x_n)=\frac{1}{2^{n-1}}{T}_n(x)$可以得到极小值，其中$T_n(x)$表示n阶切比雪夫多项式。

  我们将上述定理中的$x_i=\cos \frac{(2i-1)\pi }{2n}$记为$x_i=\cos \frac{odd\pi }{2n}$，其中odd表示1到2n之间的奇数，将这样的根记为切比雪夫根。我们将使用切比雪夫根作为基点的插值多项式叫作切比雪夫插值多项式。

切比雪夫多项式

  定义n阶切比雪夫多项式 $T_n(x)=\cos (n\arccos x)$，其中$-1\le x\le 1$。对于切比雪夫多项式，有如下归纳公式：

$$T_{n+1}(x)=2x{T}_n(x)-T_{n-1}(x)$$

由此，根据数学归纳法，不难证明，$T_n(x)$是关于x的多项式，事实上，前几个切比雪夫多项式如下：

$$\begin{gathered} T_0(x)=1 \\ {T}_1(x)=x \\ {T}_2(x)=2x^2-1 \\ {T}_3(x)=4x^3-3x \end{gathered}$$

根据切比雪夫多项式的一些基本性质，我们可以证明定理2，但这里省略，如有兴趣，可以查阅文章最后的参考文献。
  在定理2中，当x的取值区间从[-1, 1]变成任意闭区间[a, b]时，我们只需对区间做伸缩变换，则有：

  (定理3 切比雪夫插值节点)在区间[a, b]， $x=\frac{a+b}{2}+\frac{b-a}{2}\cos \frac{(2i-1)\pi }{2n},i=1,2,\cdots ,n$，不等式$|(x-x_1)\cdots (x-x_n)|\le \frac{(\frac{b-a}{2}{)}^n}{2^{n-1}}$在区间[a, b]上恒成立。

应用

应用1 观察切比雪夫点的插值多项式的拟合结果

  我们对函数$f(x)=\frac{1}{1+12x^2}$在区间[-1, 1]上，我们选取n个切比雪夫根及其对应的函数值，对这些样本点进行多项式插值。
  Python实现程序如下：

# -*- coding: utf-8 -*-
# @Time : 2023/3/8 23:40
# @Author : Jclian91
# @File : chebyshev_ploy.py
# @Place : Xuhui, Shanghai
import math
import matplotlib.pyplot as plt


# sample function
# 函数f(x)=1/(1+12*x**2)
def sample_func(x):
    return 1 / (1 + 12 * x ** 2)


# get chebyshev points from sample function with interval [-1, 1]
def get_chebyshev_points(n):
    # n: number of chebyshev points
    step = 2 / (n-1)
    x_values = [math.cos((2*i+1)*math.pi/(2*n)) for i in range(n)]
    y_values = [sample_func(x) for x in x_values]
    return x_values, y_values


# get basic lagrange polynomial unit
def get_lagrange_polynomial_unit(x_values, k, x):
    # x_values: values of x in list x_values
    # k: kth lagrange polynomial unit
    # x: variable in kth lagrange polynomial unit
    poly_unit = 1
    for i in range(len(x_values)):
        if i != k:
            poly_unit *= (x-x_values[i])/(x_values[k]-x_values[i])
    return poly_unit


# get lagrange polynomial
def get_lagrange_polynomial(x_values, y_values, x):
    poly = 0
    for i, y in enumerate(y_values):
        poly += y * get_lagrange_polynomial_unit(x_values, i, x)
    return poly


# plot curves with matplotlib
def plot_function(n):
    # plot lagrange polynomial with n sample points from sample function
    sample_x_values, sample_y_values = get_chebyshev_points(n)
    sample_points_number = 500
    x_list = [-1 + i * 2 / (sample_points_number-1) for i in range(sample_points_number)]
    original_y_list = [sample_func(x) for x in x_list]
    y_list = [get_lagrange_polynomial(sample_x_values, sample_y_values, x)for x in x_list]
    plt.plot(x_list, original_y_list, label='f(x)=1/(1+12*x**2)')
    plt.plot(x_list, y_list, label='lagrange polynomial')
    plt.title(f'Runge phenomenon with {n} chebyshev points in function f(x)=1/(1+12*x**2)')
    plt.legend()
    plt.show()
    # plt.savefig(f"{n}_basic_points.png")


if __name__ == '__main__':
    n_points = 10
    plot_function(n_points)

当n=5时，图像如下：

当n=10时，图像如下：

当n=15时，图像如下：

当n=25时，图像如下：

从中我们可以发现，当n越大，插值多项式越逼近原函数，拟合误差越小，并且我们避免了插值多项式的龙格现象，这是因为我们选取的区间点是切比雪夫点。

应用2 数值逼近

  设计程序，使得在区间$[0,\frac{\pi }{2}]$上$sin(x)$的数值精确到小数点后10位。
  在区间$[0,\frac{\pi }{2}]$上，取n个切比雪夫根，对于插值多项式$P_{n-1}(x)$在区间$[0,\frac{\pi }{2}]$上的最大插值误差有：

$$|\sin x-P_{n-1}(x)|=\frac{1}{n!}|(x-x_1)\cdots (x-x_n))||f^{(n)}(c)|\le \frac{(\frac{\frac{\pi }{2}-0}{2}{)}^n}{n!2^{n-1}}\cdot 1$$

经计算，当n=10时，误差界满足要求，可以精确到小数点后10位，此时切比雪夫根应为$\frac{\pi }{4}+\frac{\pi }{4}\cos (odd\pi /20)$.
  Python实现程序如下：

# -*- coding: utf-8 -*-
# @Time : 2023/3/9 10:37
# @Author : Jclian91
# @File : similar_to_sine_function.py
# @Place : Xuhui, Shanghai
# 利用切比雪夫多项式逼近sin(x)函数，区间为[0, pi/2]
import math

PI = math.pi    # pi in math


# sample function
# 函数f(x)=sin(x)
def sample_func(x):
    return math.sin(x)


# get chebyshev points from sample function
def get_chebyshev_points(n):
    # n: number of chebyshev points
    x_values = [PI/4 + PI/4 * math.cos((2*i+1)*PI/(2*n)) for i in range(n)]
    y_values = [sample_func(x) for x in x_values]
    return x_values, y_values


# get basic lagrange polynomial unit
def get_lagrange_polynomial_unit(x_values, k, x):
    # x_values: values of x in list x_values
    # k: kth lagrange polynomial unit
    # x: variable in kth lagrange polynomial unit
    poly_unit = 1
    for i in range(len(x_values)):
        if i != k:
            poly_unit *= (x-x_values[i])/(x_values[k]-x_values[i])
    return poly_unit


# get lagrange polynomial
def get_lagrange_polynomial(x_values, y_values, x):
    poly = 0
    for i, y in enumerate(y_values):
        poly += y * get_lagrange_polynomial_unit(x_values, i, x)
    return poly


# get similar value to sin(x) function
def get_similar_value(x):
    # x: real number in interval [0, pi/2]
    true_value = sample_func(x)
    n = 10      # chebyshev basic points number, with error < 10^-10
    x_list, y_list = get_chebyshev_points(n)
    similar_value = get_lagrange_polynomial(x_list, y_list, x)
    error = true_value - similar_value
    return f"{true_value}\t{similar_value}\t{error}"


if __name__ == '__main__':
    for x0 in [0, 0.25, 0.5, 0.75, 1, 1.25, 1.5]:
        d = get_similar_value(x0)
        print(d)

输出结果如下：

x值	sinx值	逼近值	误差
0	0.0	3.104458877467575e-11	-3.104458877467575e-11
0.25	0.24740395925452294	0.24740395924349076	1.1032175173397718e-11
0.5	0.479425538604203	0.4794255386316042	-2.7401192426168564e-11
0.75	0.6816387600233341	0.681638759993194	3.0140112627918825e-11
1	0.8414709848078965	0.8414709848397693	-3.1872837702451307e-11
1.25	0.9489846193555862	0.9489846193206646	3.492162115037445e-11
1.5	0.9974949866040544	0.9974949865890702	1.4984236074155888e-11

参考文献

1. 数值分析（原书第2版） (美)Timothy Sauer著，裴玉茹 马赓宇译， 机械工业出版社