Using the Convenience Classes

由多项式包提供的方便类是:

名称 提供
多项式 电源系列
切比雪夫 切比雪夫系列
Legendre Legendre系列
Laguerre Laguerre系列
Hermite Hermite系列
HermiteE HermiteE系列

在此上下文中的系列是相应的多项式基函数乘以系数的有限和。例如,功率系列看起来像

并具有系数[1, 2, 3]具有相同系数的切比雪夫系列看起来像

其中在这种情况下,T_n是度数n的切比雪夫函数,但可以很容易地是任何其他类的基函数。所有类的约定是系数c[i]与度i的基函数一起。

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

Basics

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

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

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

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

打印多项式产生一个较短的形式,而没有域和窗口:

>>> print p
poly([ 1.  2.  3.])

当我们得到拟合时,我们将处理域和窗口,目前我们忽略它们并且通过基本代数和算术运算。

加减:

>>> p + p
Polynomial([ 2.,  4.,  6.], [-1.,  1.], [-1.,  1.])
>>> p - p
Polynomial([ 0.], [-1.,  1.], [-1.,  1.])

乘法:

>>> p * p
Polynomial([  1.,   4.,  10.,  12.,   9.], [-1.,  1.], [-1.,  1.])

权力:

>>> p**2
Polynomial([  1.,   4.,  10.,  12.,   9.], [-1.,  1.], [-1.,  1.])

师:

floor division,'//'是多项式类的除法运算符,在这方面,多项式被视为整数。对于Python版本在某些时候,它将被弃用:

>>> p // P([-1, 1])
Polynomial([ 5.,  3.], [-1.,  1.], [-1.,  1.])

余:

>>> p % P([-1, 1])
Polynomial([ 6.], [-1.,  1.], [-1.,  1.])

Divmod:

>>> quo, rem = divmod(p, P([-1, 1]))
>>> quo
Polynomial([ 5.,  3.], [-1.,  1.], [-1.,  1.])
>>> rem
Polynomial([ 6.], [-1.,  1.], [-1.,  1.])

评价:

>>> 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代替本身,导出扩展后的新的4次多项式。如果多项式被认为是函数,则这是函数的组合:

>>> p(p)
Polynomial([  6.,  16.,  36.,  36.,  27.], [-1.,  1.], [-1.,  1.])

根:

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

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

>>> p + [1, 2, 3]
Polynomial([ 2.,  4.,  6.], [-1.,  1.], [-1.,  1.])
>>> [1, 2, 3] * p
Polynomial([  1.,   4.,  10.,  12.,   9.], [-1.,  1.], [-1.,  1.])
>>> p / 2
Polynomial([ 0.5,  1. ,  1.5], [-1.,  1.], [-1.,  1.])

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

>>> 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], [-1.,  1.], [-1.,  1.])

其给出了Chebyshev形式的多项式p这是因为T_1(x) = x而用x代替x不会改变原始多项式。然而,所有的乘法和除法将使用切比雪夫系列,因此结果的类型。

Calculus

多项式实例可以进行积分和微分。

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

第一个示例集成一次p,第二个示例集成两次。默认情况下,积分的下限和积分常数为0,但两者都可以指定。

>>> p.integ(lbnd=-1)
Polynomial([-1.,  2.,  3.], [-1.,  1.], [-1.,  1.])
>>> p.integ(lbnd=-1, k=1)
Polynomial([ 0.,  2.,  3.], [-1.,  1.], [-1.,  1.])

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

>>> p = P([1, 2, 3])
>>> p.deriv(1)
Polynomial([ 2.,  6.], [-1.,  1.], [-1.,  1.])
>>> p.deriv(2)
Polynomial([ 6.], [-1.,  1.], [-1.,  1.])

Other Polynomial Constructors

通过指定系数构造多项式仅仅是获得多项式实例的一种方式,它们也可以通过指定它们的根,通过从其他多项式类型的转换和通过最小二乘法拟合来创建。拟合在其自己的章节中讨论,其他方法如下所示:

>>> 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.], [-1.,  1.], [-1.,  1.])
>>> p.convert(kind=T)
Chebyshev([ -9.  ,  11.75,  -3.  ,   0.25], [-1.,  1.], [-1.,  1.])

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

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

在numpy版本> = 1.7.0中,还可以使用基础cast类方法。cast方法的工作方式类似于convert方法,而基础方法返回给定度的基础多项式:

>>> P.basis(3)
Polynomial([ 0.,  0.,  0.,  1.], [-1.,  1.], [-1.,  1.])
>>> T.cast(p)
Chebyshev([ -9.  ,  11.75,  -3.  ,   0.25], [-1.,  1.], [-1.,  1.])

类型之间的转换可能很有用,但建议用于日常使用。从切比雪夫系列50度到相同程度的多项式系列的数值精度的损失可使数值评估的结果基本上是随机的。

Fitting

拟合是原因,窗口属性是方便类的一部分。为了说明该问题,下面绘制直到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="$T_%d$"%i)
...
>>> plt.legend(loc="upper left")
<matplotlib.legend.Legend object at 0x3b3ee10>
>>> plt.show()

源代码pngpdf

../_images/routines-polynomials-classes-1.png

In the range -1 <= x <= 1 they are nice, equiripple functions lying between +/- 1. The same plots over the range -2 <= x <= 2 look very different:

>>> 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="$T_%d$"%i)
...
>>> plt.legend(loc="lower right")
<matplotlib.legend.Legend object at 0x3b3ee10>
>>> plt.show()

源代码pngpdf

../_images/routines-polynomials-classes-2.png

可以看出,“好”部分已经缩小到无关紧要。在使用Chebyshev多项式拟合时,我们要使用x在-1和1之间的区域,这是窗口指定的。但是,要拟合的数据不太可能具有该时间间隔中的所有数据点,因此我们使用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')
[<matplotlib.lines.Line2D object at 0x2136c10>]
>>> xx, yy = p.linspace()
>>> plt.plot(xx, yy, lw=2)
[<matplotlib.lines.Line2D object at 0x1cf2890>]
>>> p.domain
array([ 0.        ,  6.28318531])
>>> p.window
array([-1.,  1.])
>>> plt.show()

源代码pngpdf

../_images/routines-polynomials-classes-3.png