ベキ級数展開
ある関数 $f(x)$ を $x$ のベキ級数（整級数）
\[a_0+a_1(x-x_0)^1+a_2(x-x_0)^2+\cdots\]
で表すことを 関数の 級数展開 といいます。無限に項をとることによって、理論上は関数と級数展開は厳密に一致しますが、実用のために有限項で打ち切っても、適切な数の項をとることによって、ベキ級数は $x=x_0$ の近くで関数の良い近似を与えます (ただし、$x$ が $x_0$ から離れるほどその精度は下がっていきます)。たとえば、$\sin x$ のマクローリン級数
\[\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\: \cdots\]
などは級数展開の代表的な例の１つです。
sympy.series()
SymPy パッケージ の sympy.series() を使うと任意の関数をベキ級数展開できます。
以下のコードは Jupyter notebook の使用を前提に記述されています。
最初に必要な関数をまとめてインポートしておきます。
# SYMPY_SERIES
# In[1]
from sympy import init_printing, var, series, poly, cos, sin
# 数式をLaTeX表示
init_printing()
# シンボルxを定義
var('x')$\sin x$ を $x=0$ の近くで級数展開してみましょう。
# In[2]
# sinxをx=0の周りで級数展開
series(sin(x), x)デフォルト設定では、$x^6$ 以降の項はすべて O(x**6) にまとめられています。オプション引数 n でこの設定を変更することができます。たとえば、O(x**7) まで級数を表示したい場合は次のように記述します。
# In[3]
# sinxをx=0の周りで級数展開
series(sin(x), x, n=8)オプション引数 x0 で展開の中心点を指定できます (デフォルトは 0)。$\sin x$ を $x=1$ のまわりで展開してみます。
# In[4]
# sinxをx=1の周りで級数展開
series(sin(x), x, x0=1)\sin{\left (1 \right )} &+ \left(x - 1\right) \cos{\left (1 \right )} - \frac{\left(x - 1\right)^{2} \sin{\left (1 \right )}}{2} - \frac{\left(x - 1\right)^{3} \cos{\left (1 \right )}}{6}\\[6pt]
&+ \frac{\left(x - 1\right)^{4} \sin{\left (1 \right )}}{24} + \frac{\left(x - 1\right)^{5} \cos{\left (1 \right )}}{120} + O\left(\left(x - 1\right)^{6}; x\rightarrow 1\right)\end{align*}\]
sympy.series() で得た式を sympy.poly() で多項式クラス (sympy.polys.polytools.Polyクラス) のオブジェクトに変換しておくと、色々なメソッドが使えるようになって便利です。$\cos x$ を級数展開して、coeffs()メソッドで $x^n$ の係数をすべて取り出してみます。
# In[5]
# cosxをx=0の周りで級数展開
f = series(cos(x), x, n=10)
# 級数展開された多項式fを表示
display(f)
# fを多項式オブジェクトに変換
g = poly(f)
# 級数の係数をすべて取得
c = g.coeffs()
c.reverse()
# 級数の係数を表示
display(c)
母関数の級数展開
級数 $a_n$ が関数 $p(x)$ を級数展開したときの $x^n$ の係数となっているとき、$p(x)$ は級数 $a_n$ の母関数であるといいます（「級数を生み出す関数」という意味です）。たとえば、前の $2$ つの項を足して次の項をつくるフィボナッチ数列
\[1,\ 1,\ 2,\ 3,\ 5,\ 8,\ 13,\ 21,\ 34,\ \cdots\]
の母関数は
\[p(x)=\frac{1}{1-x-x^2}\]
で与えられることが知られています。$p(x)$ を sympy.series() で展開して確認してみましょう。
# In[6]
# フィボナッチ数列の母関数
p = 1/(1 - x - x**2)
# フィボナッチ数列の母関数を級数展開
f = series(p, x, n=10)
display(f)係数は確かにフィボナッチ数列となっていることがわかります。
# In[7]
# fを多項式オブジェクトに変換
g = poly(f)
# フィボナッチ級数のすべての係数を取得
c = g.coeffs()
c.reverse()
display(c)