『Python数値計算ノート』ではアフィリエイトプログラムを利用して商品を紹介しています。

【SymPy】級数展開

ベキ級数展開

ある関数 f(x)x のベキ級数(整級数)
 a0+a1(xx0)1+a2(xx0)2+
で表すことを 関数の 級数展開 といいます。無限に項をとることによって、理論上は関数と級数展開は厳密に一致しますが、実用のために有限項で打ち切っても、適切な数の項をとることによって、ベキ級数は x=x0 の近くで関数の良い近似を与えます (ただし、xx0 から離れるほどその精度は下がっていきます)。たとえば、sinx のマクローリン級数
 sinx=xx33!+x55!x77!+
などは級数展開の代表的な例の1つです。

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

sinxx=0 の近くで級数展開してみましょう。

# In[2]

# sinxをx=0の周りで級数展開
series(sin(x), x)
xx36+x5120+O(x6)

デフォルト設定では、x6 以降の項はすべて O(x**6) にまとめられています。オプション引数 n でこの設定を変更することができます。たとえば、O(x**7) まで級数を表示したい場合は次のように記述します。

# In[3]

# sinxをx=0の周りで級数展開
series(sin(x), x, n=8)
xx36+x5120x75040+O(x8)

オプション引数 x0 で展開の中心点を指定できます (デフォルトは 0)。sinxx=1 のまわりで展開してみます。

# In[4]

# sinxをx=1の周りで級数展開
series(sin(x), x, x0=1)
sin(1)+(x1)cos(1)(x1)2sin(1)2(x1)3cos(1)6+(x1)4sin(1)24+(x1)5cos(1)120+O((x1)6;x1)

sympy.series() で得た式を sympy.poly() で多項式クラス (sympy.polys.polytools.Polyクラス) のオブジェクトに変換しておくと、色々なメソッドが使えるようになって便利です。cosx を級数展開して、coeffs()メソッドで xn の係数をすべて取り出してみます。

# 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)
1x22+x424x6720+x840320+O(x10)[1,1,12,124,1720,140320]

母関数の級数展開

級数 an が関数 p(x) を級数展開したときの xn の係数となっているとき、p(x) は級数 an の母関数であるといいます(「級数を生み出す関数」という意味です)。たとえば、前の 2 つの項を足して次の項をつくるフィボナッチ数列
 1, 1, 2, 3, 5, 8, 13, 21, 34, 
の母関数は
 p(x)=11xx2
で与えられることが知られています。p(x) を sympy.series() で展開して確認してみましょう。

# In[6]

# フィボナッチ数列の母関数
p = 1/(1 - x - x**2)

# フィボナッチ数列の母関数を級数展開
f = series(p, x, n=10)

display(f)
1+x+2x2+3x3+5x4+8x5+13x6+21x7+34x8+55x9+O(x10)

係数は確かにフィボナッチ数列となっていることがわかります。

# In[7]

# fを多項式オブジェクトに変換
g = poly(f)

# フィボナッチ級数のすべての係数を取得
c = g.coeffs()
c.reverse()

display(c)
[1,1,1,2,3,5,8,13,21,34,55]

コメント

  1. tiida より:

    お尋ねしますが,
    オーダーの記号の設定をすることは可能でしょうか?
    例えば,上記であれば,
    sin xのx=1のまわりのテイラー展開で最後の項が
    O((x-1)^6;x→1)
    となっていますが,これを
    O((x-1)^5)
    等に変更できるのでしょうか?できればOも\mathcal{O}ではなく小文字のoに変更できればなおよいのですが.

    • BlogCat より:

      sympy.sereis() はオプション引数 n で展開打ち切り次数を設定できます。
      デフォルトで n=6 に設定されていますが、たとえば
      series(sin(x), x, x0=1, n=5)
      とすれば、最後の項が O((x-1)^5) となります。
      O を小文字にする方法はわかりません。申し訳ないです。