【NumPy】多項式の計算

1次元多項式

1次元多項式オブジェクト

　numpy.poly1d() は 1次元多項式オブジェクト (numpy.lib.polynomial.poly1dクラスのインスタンス) を生成します。

np.poly1d(c_or_r, r=False, variable=None)

　c_or_r には係数または根を渡します。第 2 引数 r が False ならば c (係数) 、True の場合は r (根) が渡されたことになります。第 3 引数 variable には変数として用いる記号を渡します (デフォルトは x です)。
　
　たとえば、$P(x)=x^2+2x+3$ という多項式を作る場合は、第 1 引数に各項の係数のリスト [1, 2, 3] を渡します。

# NUMPY_POLY1D

# In[1]

import numpy as np

# p1 = x**2 + 2*x + 3
p1 = np.poly1d([1, 2, 3])

print(p1)
#    2
# 1 x + 2 x + 3

　print() を使うと二行で表示されます。少しわかりにくいのですが、上の行にある 2 という数字が $x$ の指数を表しています。Jupyter Notebook の display()関数を使うと poly1d([1, 2, 3]) のように表示されます。
　
　以下、1次元多項式オブジェクトを p と略記します。
　p.o または p.order で多項式の次数を取得できます。

# In[2]

# p1 = x**2 + 2*x + 3 の次数を取得
print(p1.order)

# 2

　p.c, p.coef, p.coefficients で各項の係数を取得できます。

# In[3]

# p1 = x**2 + 2*x + 3 の各項の係数を取得
print(p1.coef)
# [1 2 3]

　p.variable は変数として使用されている記号です。

# In[4]

# p1 = x**2 + 2*x + 3 の変数を取得
print(p1.variable)

# x

　p.r で「多項式 = 0 」の根（解）を得ることができます。

# In[5]

# p1 = x**2 + 2*x + 3 = 0 の根を取得
print(p1.r)

# [-1.+1.41421356j -1.-1.41421356j]

　numpy.poly1d() の第 2 引数 r を True にすると、第 1 引数に渡した配列は根であるとみなされます。たとえば、[1, 2] を渡すと (x - 1) * (x - 2) という多項式が生成されます。

# In[6]

# p2 = (x-1)*(x-2)
p2 = np.poly1d([1, 2], r = True)

print(p2)
#    2
# 1 x - 3 x + 2

　第 3 引数 (variable) で多項式の変数として使用する記号を指定することができます。

# In[7]

# p3 = (z-1)*(z-2)*(z-3)
p3 = np.poly1d([1, 2, 3], r = True, variable = "z")

print(p3)
#    3     2
# 1 z - 6 z + 11 z - 6

多項式の演算と代入

　算術演算子によって、数値と多項式、あるいは多項式オブジェクト同士で加減乗除を行なうことができます。

# NUMPY_POLY1D_CALCULATION

# In[1]

import numpy as np

# p = x + 1
p = np.poly1d([1, 1])

# q = 3*x**2 + 5*x + 4
q = np.poly1d([3, 5, 4])

print(p + q + 10)
#    2
# 3 x + 6 x + 15

　多項式同士の除算は、(商, 余り) のタプルが返ります。
　たとえば $P(x)=x+1,\ Q(x)=3x^2+5x+4$ のとき、
　
\[\frac{Q(x)}{P(x)}=3x+2+\frac{2}{x+1}\]
なので、商は $3x+2$, 余りは $2$ となりますが、多項式オブジェクトはこれを

(poly1d([3., 2.]), poly1d([2.]))

のように表します。

# In[2]

# P = x + 1, q = 3*x**2 + 5*x + 4

print(q / p)
# (poly1d([3., 2.]), poly1d([2.]))

　多項式に数値を代入することもできます。

# In[3]

# q = x**2 + 3*x + 1 に x=5 を代入
print(q(5))

# 41

　配列を代入することもできます。

# In[4]

# q = x**2 + 3*x + 1 に x=[1,2,3] を代入
print(q(np.arange(1, 4)))

# [ 5 11 19]

　多項式に多項式を代入することも可能です。

# In[5]

# q=x**2+3*x+1にp=2*x+1を代入
print(q(p))

#    2
# 4 x + 10 x + 5

多項式の微分と積分

　deriv()メソッド で多項式を微分することができます。

# NUMPY_POLY1D_DERIV

# In[1]

import numpy as np

# y = 2x**3 + 5x**2 + 3*x + 7
y = np.poly1d([2, 5, 3, 7])

# yの1階導関数
d_y = y.deriv()

print(d_y)
#    2
# 6 x + 10 x + 3

　p.deriv(k) は k階微分を取得します。

# In[2]

# y = 2x**3 + 5x**2 + 3*x + 7

# yの2階導関数
d2_y = y.deriv(2)

print(d2_y)
# 12 x + 10

　integ()メソッドは多項式を積分します。

# NUMPY_POLY1D_INTEG

# y = 2x**3 + 5x**2 + 3*x + 7 を積分
i_y = y.integ()

print(i_y)
#      4         3       2
# 0.5 x + 1.667 x + 1.5 x + 7 x

　引数 k で積分定数を指定することができます。

# In[4]

# y = 2x**3 + 5x**2 + 3*x + 7 を積分(積分定数5)
i_y = y.integ(k = 5)

print(i_y)
#      4         3       2
# 0.5 x + 1.667 x + 1.5 x + 7 x + 5