使用便捷类#
polynomial 包提供的便捷类有:
在这种上下文中,级数是相应的多项式基函数乘以系数的有限和.例如,幂级数看起来像
且有系数 \([1, 2, 3]\) .具有相同系数的切比雪夫级数如下所示
更一般地
其中在此情况下, \(T_n\) 是 \(n\) 次切比雪夫函数,但同样可以是任何其他类的基函数.所有类的约定是系数 \(c[i]\) 与 i 次基函数相关联.
所有类都是不可变的,并且具有相同的方法,特别是,它们实现了 Python 数值运算符 +, -, , //, %, divmod, \, == 和 !=.由于浮点舍入误差,最后两个可能会有些问题.现在,我们将使用 NumPy 1.7.0 版本快速演示各种操作.
基础知识#
首先,我们需要一个多项式类和一个多项式实例来使用.这些类可以直接从 polynomial 包或相关类型的模块中导入.这里我们从包中导入,并使用传统的 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')
幂:
>>> p**2
Polynomial([ 1., 4., 10., 12., 9.], domain=[-1., 1.], window=[-1., 1.], symbol='x')
除法:
向下取整除法“//”是多项式类的除法运算符,在这种情况下,多项式被视为整数.对于 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 替换为自身,从而在展开后得到一个新的 4 次多项式.如果将多项式视为函数,这就是函数的组合:
>>> 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])
显式使用 Polynomial 实例并不总是方便,因此元组,列表,数组和标量会在算术运算中自动转换:
>>> 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
但是不同的类型可以用于替换.事实上,这就是在类型,域和窗口转换中如何在 Polynomial 类之间进行转换的方式:
>>> p(T([0, 1]))
Chebyshev([2.5, 2. , 1.5], domain=[-1., 1.], window=[-1., 1.], symbol='x')
这将给出切比雪夫形式的多项式 p .这是可行的,因为 \(T_1(x) = x\) 并且用 \(x\) 替换 \(x\) 不会改变原始多项式.但是,所有乘法和除法都将使用切比雪夫级数完成,因此结果的类型也是切比雪夫级数.
所有多项式实例都应是不可变的,因此,增广运算( += , -= 等)和任何其他会违反多项式实例不变性的功能都有意未实现.
微积分#
可以对多项式实例进行积分和微分.:
>>> 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 方法也可以转换 domain 和 window:
>>> 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()
在 -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()
可以看出,“好”的部分已经缩小到无关紧要的地步.在使用切比雪夫多项式进行拟合时,我们希望使用 x 在 -1 和 1 之间的区域,这就是 window 指定的内容.但是,要拟合的数据不太可能将其所有数据点都位于该区间内,因此我们使用 domain 来指定数据点所在的区间.完成拟合后,首先通过线性变换将 domain 映射到 window,然后使用映射的数据点完成通常的最小二乘拟合.拟合的 window 和 domain 是返回的序列的一部分,并在计算值,导数等时自动使用.如果在调用中未指定它们,则拟合例程将使用默认的 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()
