使用便利类#

多项式包提供的便利类有:

名字

提供

Polynomial

幂级数

Chebyshev

切比雪夫级数

Legendre

勒让德级数

Laguerre

Laguerre 系列

Hermite

埃尔米特级数

HermiteE

HermiteE 系列

在这种情况下,级数是相应的多项式基函数乘以系数的有限和.例如,幂级数看起来像

\[p(x) = 1 + 2x + 3x^2\]

并且有系数 \([1, 2, 3]\).具有相同系数的切比雪夫级数看起来像

\[p(x) = 1 T_0(x) + 2 T_1(x) + 3 T_2(x)\]

更一般地

\[p(x) = \sum_{i=0}^n c_i T_i(x)\]

在这种情况下,:math:T_n 是度数为 \(n\) 的切比雪夫函数,但同样可以是其他类别的基础函数.所有类别的约定是系数 \(c[i]\) 与度数为 i 的基础函数相对应.

所有的类都是不可变的,并且具有相同的方法,特别是它们实现了 Python 的数值运算符 +, -, *, //, %, divmod, **, ==, 和 !=.由于浮点舍入误差,最后两个可能会有点问题.我们现在使用 NumPy 版本 1.7.0 快速演示各种操作.

基础#

首先我们需要一个多项式类和一个多项式实例来操作.这些类可以直接从多项式包中导入,或者从相关类型的模块中导入.这里我们从包中导入并使用常规的 Polynomial 类,因为它的熟悉度:

>>> from numpy.polynomial import Polynomial as P
>>> p = P([1,2,3])
>>> p
Polynomial([1., 2., 3.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

请注意,长版本的打印输出有三个部分.第一部分是系数,第二部分是域,第三部分是窗口:

>>> p.coef
array([1., 2., 3.])
>>> p.domain
array([-1.,  1.])
>>> p.window
array([-1.,  1.])

打印多项式会生成更熟悉格式的多项式表达式:

>>> print(p)
1.0 + 2.0·x + 3.0·x²

请注意,多项式的字符串表示默认使用Unicode字符（Windows除外）来表示幂和下标.也可以使用基于ASCII的表示（Windows默认）.多项式的字符串格式可以在包级别通过 set_default_printstyle 函数进行切换:

>>> np.polynomial.set_default_printstyle('ascii')
>>> print(p)
1.0 + 2.0 x + 3.0 x**2

或通过字符串格式化控制单个多项式实例:

>>> print(f"{p:unicode}")
1.0 + 2.0·x + 3.0·x²

当我们进行拟合时,我们将处理域和窗口,目前我们忽略它们并运行基本的代数和算术运算.

加法和减法:

>>> p + p
Polynomial([2., 4., 6.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
>>> p - p
Polynomial([0.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

乘法:

>>> p * p
Polynomial([ 1.,   4.,  10.,  12.,   9.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

Powers:

>>> p**2
Polynomial([ 1.,   4., 10., 12.,  9.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

Division:

地板除法,’//’,是多项式类的除法运算符,在这方面多项式被视为整数.对于Python版本< 3.x,’/’运算符映射到’//’,就像在Python中一样,对于更高版本,’/’仅适用于标量除法.在某些时候它将被弃用:

>>> p // P([-1, 1])
Polynomial([5.,  3.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

余数:

>>> p % P([-1, 1])
Polynomial([6.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

Divmod:

>>> quo, rem = divmod(p, P([-1, 1]))
>>> quo
Polynomial([5.,  3.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
>>> rem
Polynomial([6.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

评估:

>>> x = np.arange(5)
>>> p(x)
array([  1.,   6.,  17.,  34.,  57.])
>>> x = np.arange(6).reshape(3,2)
>>> p(x)
array([[ 1.,   6.],
       [17.,  34.],
       [57.,  86.]])

替换:

将多项式代入 x 并展开结果.这里我们将 p 代入其自身,导致展开后得到一个新的四次多项式.如果将多项式视为函数,这就是函数的复合:

>>> p(p)
Polynomial([ 6., 16., 36., 36., 27.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

根:

>>> p.roots()
array([-0.33333333-0.47140452j, -0.33333333+0.47140452j])

在算术运算中,显式使用多项式实例并不总是方便的,因此元组、列表、数组和标量会自动转换:

>>> p + [1, 2, 3]
Polynomial([2., 4., 6.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
>>> [1, 2, 3] * p
Polynomial([ 1.,  4., 10., 12.,  9.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
>>> p / 2
Polynomial([0.5, 1. , 1.5], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

在域、窗口或类上不同的多项式不能在算术中混合:

>>> from numpy.polynomial import Chebyshev as T
>>> p + P([1], domain=[0,1])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<string>", line 213, in __add__
TypeError: Domains differ
>>> p + P([1], window=[0,1])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<string>", line 215, in __add__
TypeError: Windows differ
>>> p + T([1])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<string>", line 211, in __add__
TypeError: Polynomial types differ

但不同的类型可以用于替换.事实上,这就是如何在类型、域和窗口之间转换多项式类的方式:

>>> p(T([0, 1]))
Chebyshev([2.5, 2. , 1.5], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

这给出了Chebyshev形式的`p`多项式.这是因为 \(T_1(x) = x\) 并且将 \(x\) 替换为 \(x\) 不会改变原始多项式.然而,所有的乘法和除法都将使用Chebyshev级数进行,因此结果的类型.

所有多项式实例都应是不可变的,因此增强操作（+=、-= 等）以及任何可能违反多项式实例不可变性的其他功能都故意未实现.

微积分#

多项式实例可以被积分和微分.:

>>> from numpy.polynomial import Polynomial as P
>>> p = P([2, 6])
>>> p.integ()
Polynomial([0., 2., 3.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
>>> p.integ(2)
Polynomial([0., 0., 1., 1.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

第一个示例将 p 积分一次,第二个示例将其积分两次.默认情况下,积分的下限和积分常数均为 0,但两者都可以指定.:

>>> p.integ(lbnd=-1)
Polynomial([-1.,  2.,  3.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
>>> p.integ(lbnd=-1, k=1)
Polynomial([0., 2., 3.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

在第一种情况下,积分的下限设置为 -1,积分常数为 0.在第二种情况下,积分常数也设置为 1.微分更简单,因为唯一的选项是多项式被微分的次数:

>>> p = P([1, 2, 3])
>>> p.deriv(1)
Polynomial([2., 6.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
>>> p.deriv(2)
Polynomial([6.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

其他多项式构造函数#

通过指定系数来构造多项式只是获得多项式实例的一种方法,它们也可以通过指定它们的根、从其他多项式类型转换以及通过最小二乘拟合来创建.拟合在其自己的部分中讨论,其他方法在下面演示:

>>> from numpy.polynomial import Polynomial as P
>>> from numpy.polynomial import Chebyshev as T
>>> p = P.fromroots([1, 2, 3])
>>> p
Polynomial([-6., 11., -6.,  1.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
>>> p.convert(kind=T)
Chebyshev([-9.  , 11.75, -3.  ,  0.25], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

convert 方法也可以转换域和窗口:

>>> p.convert(kind=T, domain=[0, 1])
Chebyshev([-2.4375 ,  2.96875, -0.5625 ,  0.03125], domain=[0.,  1.], window=[-1.,  1.], symbol='x')
>>> p.convert(kind=P, domain=[0, 1])
Polynomial([-1.875,  2.875, -1.125,  0.125], domain=[0.,  1.], window=[-1.,  1.], symbol='x')

在 numpy 版本 >= 1.7.0 中,`basis` 和 cast 类方法也是可用的.cast 方法的工作方式类似于 convert 方法,而 basis 方法返回给定次数的基多项式:

>>> P.basis(3)
Polynomial([0., 0., 0., 1.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
>>> T.cast(p)
Chebyshev([-9.  , 11.75, -3. ,  0.25], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')

类型之间的转换可能是有用的,但*不*推荐常规使用.从50阶切比雪夫级数转换为相同阶的多项式级数时,数值精度的损失会使数值评估的结果基本上变得随机.

拟合#

拟合是 domain 和 window 属性成为便捷类一部分的原因.为了说明这个问题,下面绘制了直到5次的多项式切比雪夫值.

>>> import matplotlib.pyplot as plt
>>> from numpy.polynomial import Chebyshev as T
>>> x = np.linspace(-1, 1, 100)
>>> for i in range(6):
...     ax = plt.plot(x, T.basis(i)(x), lw=2, label=f"$T_{i}$")
...
>>> plt.legend(loc="upper left")
>>> plt.show()
../_images/routines-polynomials-classes-1.png

在范围 -1 <= x <= 1 内,它们是很好的等波纹函数,位于 +/- 1 之间.在范围 -2 <= x <= 2 内的相同图看起来非常不同:

>>> import matplotlib.pyplot as plt
>>> from numpy.polynomial import Chebyshev as T
>>> x = np.linspace(-2, 2, 100)
>>> for i in range(6):
...     ax = plt.plot(x, T.basis(i)(x), lw=2, label=f"$T_{i}$")
...
>>> plt.legend(loc="lower right")
>>> plt.show()
../_images/routines-polynomials-classes-2.png

可以看出,”好”的部分已经缩减到微不足道.在使用切比雪夫多项式进行拟合时,我们希望使用 x 在 -1 到 1 之间的区域,这就是 window 指定的内容.然而,要拟合的数据不太可能所有数据点都在这个区间内,所以我们使用 domain 来指定数据点所在的区间.当进行拟合时,首先通过线性变换将域映射到窗口,然后使用映射后的数据点进行通常的最小二乘拟合.拟合的窗口和域是返回系列的一部分,并在计算值、导数等时自动使用.如果它们没有在调用中指定,拟合例程将使用默认窗口和包含所有数据点的最小域.下面是一个对噪声正弦曲线进行拟合的示例.

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from numpy.polynomial import Chebyshev as T
>>> np.random.seed(11)
>>> x = np.linspace(0, 2*np.pi, 20)
>>> y = np.sin(x) + np.random.normal(scale=.1, size=x.shape)
>>> p = T.fit(x, y, 5)
>>> plt.plot(x, y, 'o')
>>> xx, yy = p.linspace()
>>> plt.plot(xx, yy, lw=2)
>>> p.domain
array([0.        ,  6.28318531])
>>> p.window
array([-1.,  1.])
>>> plt.show()
../_images/routines-polynomials-classes-3.png
On this page