使用便利類別#

多項式套件提供的便利類別有

名稱

提供

Polynomial

冪級數

Chebyshev

切比雪夫級數

Legendre

勒讓德級數

Laguerre

拉蓋爾級數

Hermite

埃爾米特級數

HermiteE

埃爾米特E級數

此處的級數是相應多項式基底函數乘以係數的有限總和。例如,冪級數看起來像

\[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)\]

其中在此例中 \(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')

次方

>>> 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

但不同類型可以用於替換。事實上,這就是多項式類別之間針對類型、定義域和視窗轉換的方式

>>> 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 方法也可以轉換定義域和視窗

>>> 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 中,basiscast 類別方法也可用。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 次的切比雪夫級數傳遞到相同次數的多項式級數會導致數值評估結果的數值精度損失,基本上是隨機的。

擬合#

擬合是 domainwindow 屬性成為便利類別一部分的原因。為了說明問題,下面繪製了直到 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