入力例で学ぶPython (SymPy) の使い方（基礎）

本記事ではプログラミング言語『Python』の数式処理ライブラリ『SymPy』の使い方を紹介します。

以前の記事ではPython (SymPy) の使い方のごく基本を紹介しました：

今回は、そこで学んだことの具体的な使用例としてBernoulli数を求めてみたいと思います。

Bernoulli数は関数

$\displaystyle{f(x)=\frac{x}{e^x-1}}$

を $x=0$ のまわりでTaylor展開した

$\displaystyle{f(x)=\sum_{n=0}^\infty f^{(n)}(0) \frac{x^n}{n!}}$

の展開係数として求められます：

$\displaystyle{B_{n}=f^{(n)}(0)}$

つまり

$\displaystyle{\frac{x}{e^x-1}=\sum_{n=0}^\infty B_{n}\frac{x^{n}}{n!}}$

です。

まずSymPyを読み込み、関数 $f(x)$ を定義します：

import sympy as sp
sp.init_printing()
sp.var('x')
f=x/(sp.exp(x)-1)

関数 $f(x)$ の $x=0$ のまわりの11次まで（10次以下）のTaylor展開が

sp.series(f,x, 0, 11)

で求められます。

Taylor展開の $k$ 次の係数は

b_list = []
for i in range(11):
    b = sp.series(f,x, 0, i+1).coeff(x, i)
    b_list.append(b)

によってリスト化することができます。

Beroulli数は上記関数のTaylor展開の $k$ 次の係数に $k!$ を掛けたものなので

上で得たリストの各要素を $k!$ 倍したリスト

B_list = []
for i in range(11):
    B = b_list[i]*sp.factorial(i)
    B_list.append(B)

を作成します。

最後に、このリストを表示してみます：

print(B_list)

すると数列

$\displaystyle{~1, ~-\frac{1}{2},~\frac{1}{6},~0,~-\frac{1}{30},~0,~\frac{1}{42},~0,~-\frac{1}{30},~0,~\frac{5}{66},~\dots}$

を得ます。これらがBernoulli数です。

実は、わざわざ上のようにして求めなくても、

SymPyにはBernoulli数がすでに組み込まれています：たとえば

sp.bernoulli(2)

などです。

そこで、Taylor展開の係数から求めたリスト B_list の各要素に対して、

それが組み込みのBernoulli数に等しいか？という真偽判定をさせてみます：

TF_list = []
for i in range(11):
    tf = (0==sp.bernoulli(i)-B_list[i])
    TF_list.append(tf)
print(TF_list)

11個すべての要素について一致することが確認できたと思います。

本記事で扱った内容を数学的にもPython的にも、もう少し発展させたのが以下の記事です：

本記事の続編として是非あわせて読んでみてください。

また、本記事で紹介した内容をMaximaで計算させているのが以下の記事です：

比較しながらあわせて読むと、MaximaとPythonの違いがわかるかと思います。

キーワード：Python、SymPy、Bernoulli数（ベルヌーイ数）