pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

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

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

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

pianofisica.hatenablog.com

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



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

です。

Bernoulli数をSymPyで求める

まず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的にも、もう少し発展させたのが以下の記事です:

pianofisica.hatenablog.com

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


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

pianofisica.hatenablog.com

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


キーワードPython、SymPy、Bernoulli数(ベルヌーイ数)