Python处理第一类切比雪夫多项式

第一类切比雪夫多项式简介

Chebyshev多项式是一种非常重要的正交多项式在逼近理论中有重要应用第一类切比雪夫多项式的根可用于多项式插值对弥补龙格现象有很大的帮助。其表达形式为

$$T_n=\cos (n\arccos x)$$

其中$n$为第一类切比雪夫多项式的阶数在Python中提供了Chebyshev类类构造函数为

chebyshev.Chebyshev(coef, domain=None, window=None, symbol='x')

其中coef为系数列表$a_0,a_1,\cdots ,a_n$表示生成

$$\sum _{i=0}^n{a}_i{T}_i(x)$$

domain表示$x$的定义域window表示缩放系数x为自变量符号。

from numpy.polynomial.chebyshev import Chebyshev
t3 = Chebyshev(coef=[4,3,2,1])
print(t3)
# 输出为4.0 + 3.0 T_1(x) + 2.0 T_2(x) + 1.0 T_3(x)

为了对第一类切比雪夫有个直观的认识可以绘制一下不同阶数的第一类切比雪夫的函数

import numpy as np
import matplotlib.pyplot as plt
for i in range(5):
    c = np.zeros(i+1)
    c[i] = 1
    t = Chebyshev(coef=c, domain=(-5,5))
    xs, ys = t.linspace()
    plt.plot(xs, ys, label=str(i))

plt.legend()
plt.show()

其中h.linspace表示在定义域范围内对多项式进行采样有一个参数n表示在定义域范围内等间隔生成n组$x,y$默认为100。

得图如下

求导和积分

Chebyshev支持简单的符号计算比如可通过deriv(n)求多项式的n阶导数；通过integ(n)可求n阶积分示例如下

>>> t3.deriv(1)
Chebyshev([6., 8., 6.], domain=[-1.,  1.], window=[-1.,  1.])
>>> t3.deriv(2)
Chebyshev([ 8., 24.], domain=[-1.,  1.], window=[-1.,  1.])
>>> t3.integ(1)
Chebyshev([0.375     , 3.        , 0.5       , 0.33333333, 0.125     ], domain=[-1.,  1.], window=[-1.,  1.])

求根和反演

roots可用于求根而fromroot可根据根来生成切比雪夫多项式

>>> rs = t3.roots()
>>> print(rs)
[-1.29715651+0.j          0.14857825-0.60281258j  0.14857825+0.60281258j]
>>> tNew = t3.fromroots(rs)
>>> print(tNew)
(0.9999999999999998+0j) + (0.749999999999999+0j) T_1(x) +
(0.49999999999999967+0j) T_2(x) + (0.25+0j) T_3(x)

可以发现roots和fromroots并非对称的关系。

拟合

Chebyshev类中同样提供了拟合函数fit定义为

Chebyshev.fit(x, y, deg, domain=None, rcond=None, full=False, w=None, window=None, symbol='x')

其中domain, window, symbol不必赘述其中x,y为待拟合多项式；deg为多项式的阶数。rcond表示截止误差。full为False时只返回拟合系数否则还返回拟合的标准差等。

>>> t3 = Chebyshev(coef=[4,3,2,1])
>>> xs, ys = p3.linspace()
>>> t3_3 = t3.fit(xs, ys, 3)
>>> print(t3_3)
3.9999999999999987 + 3.0000000000000018 T_1(x) +
1.999999999999999 T_2(x) + 1.0000000000000013 T_3(x)
>>> t3_4 = t3.fit(xs, ys, 4)
>>> print(t3_4)
4.0 + 3.0000000000000004 T_1(x) + 2.0000000000000013 T_2(x) +
1.0000000000000007 T_3(x) + 2.9674017981841603e-15 T_4(x)

可见其拟合效果还是不错的。

degree返回多项式的最高项次数cutdeg可以对多项式的次数做阶段例如

>>> t3.degree()
3
>>> t3.cutdeg(2)
Chebyshev([4., 3., 2.], domain=[-1.,  1.], window=[-1.,  1.])