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

数学の具体的な計算にPythonを使って、数学もPythonも同時に学んでしまいましょう。

本記事は、以前の記事

pianofisica.hatenablog.com

で学んだ内容の（数学的にもPython的にも）発展編です。

今回はPython（SymPy）を使ってBernoulli多項式（ベルヌーイ多項式）を学ぼうと思います。

Bernoulli多項式（ベルヌーイ多項式）
- 定義
- 性質
- 母関数

Bernoulli多項式（ベルヌーイ多項式）

数学では、有用な性質を持っていることから特別な名前がつけられている多項式がたくさんあります。

今回はその中でもBernoulli多項式と呼ばれるものを取り上げたいと思います。

Bernoulli多項式は、スイスの数学者 Jacob Bernoulli（1654 – 1705）にちなんだ多項式です。

以下で具体的にみていくように、二項係数とBernoulli数（ベルヌーイ数）を組み合わせた多項式で、

整数の冪乗和とも深く関係しています。

例えば、高校の数学でも出てくる、自然数 $1$ から $N$ までの和

$\quad\displaystyle{S_1(N)=\sum_{n=1}^Nn=1+2+3+\dots+N}$

は、まだ子供だった J. C. F. Gauss（1777 - 1855）が $1$ から $50$ の場合について、

その答えを即座に求めて先生を驚かせたという逸話で彩られている有名な公式

$\quad\displaystyle{S_1(N)=\frac{N(N+1)}{2}}$

で与えられます。では、自然数 $1$ から $N$ までの二乗の和

$\quad\displaystyle{S_2(N)=\sum_{n=1}^Nn^2=1^2+2^2+3^2+\dots+N^2}$

はどのような公式になったかというと（高校数学の公式を覚えているでしょうか？）

$\quad\displaystyle{S_2(N)=\frac{N(N+1)(2N+1)}{6}}$

です（数学的帰納法などで示すことができます）。

さて、これらの公式で右辺を $N$ の多項式だと思うことにしましょう。

この多項式は何らかの系統的な手続きで与えられるでしょうか？

実はこれらの多項式は、Bernoulli多項式 $B_n$ と

$\quad\displaystyle{S_n(N)=\int_1^{N+1}B_n(x)\,dx}$

によって関係付けられます。さらに言えば

$\quad\displaystyle{k^n=\int_k^{k+1}B_n(x)\,dx}$

という関係が成り立ち、積分の区間と和の性質により、上の公式が導かれることがわかります。

以下では、Bernoulli多項式 $B_n$ がどのように与えられるのか、

冪乗和との関係がどのようにして導かれるのか、詳しく見ていくことにしましょう。

Pythonで具体的な計算を入力してその結果をみながら読み進めてもらうことで、

理解をより深めてもらうための一助になればと思います。

定義

Bernoulli多項式の定義式は

$\quad\displaystyle{B_n(x)=\sum_{k=0}^n\binom{n}{k}b_{n-k}x^k}$

です。ここで $\binom{n}{k}$ は二項係数

$\quad\displaystyle{\binom{n}{k}={}_nC_k=\frac{n!}{k!(n-k)!}}$

そして $b_{n}$ はBernoulli数です。

冒頭にも述べましたが、本記事は以前の記事

pianofisica.hatenablog.com

の発展編です。Bernoulli数に関してはそこで取り上げているので、

あわせて読んでいただけたらと思います。

定義式から具体的にいくつかBernoulli多項式を求めてみます：

import sympy as sp
sp.init_printing()
k,n=sp.symbols('k,n',integer=True)
sp.var('x')
def B(n,x):
   return sp.summation(sp.binomial(n,k)*sp.bernoulli(n-k)*x**k,(k,0,n))
B_list = []
for n in range(7):
    b = B(n,x)
    B_list.append(b)
B_list

から

$\quad\displaystyle{ \begin{array}{l} B_{0}(x)=1 \\ B_{1}(x)=x-\frac{1}{2} \\ B_{2}(x)=x^{2}-x+\frac{1}{6} \\ B_{3}(x)=x^{3}-\frac{3}{2} x^{2}+\frac{1}{2} x \\ B_{4}(x)=x^{4}-2 x^{3}+x^{2}-\frac{1}{30}\\ \quad\vdots \end{array} }$

です。特に引数で $x=0$ とおくと定義式の和の初項のみが値を持って

$\quad\displaystyle{ B_n(0)=b_n }$

より、Bernoulli数を与えます：

Ber_list = []
for n in range(15):
    ber = B(n,0)
    Ber_list.append(ber)
Ber_list

$\quad\displaystyle{ b_0=1 \qquad b_1=-\frac12 \qquad b_2=\frac16 \qquad b_3=0 \qquad b_4=-\frac{1}{30} \quad \dots }$

性質

Bernoulli多項式の特筆すべき性質は（既に述べたように）

$\quad\displaystyle{\int_{k}^{k+1} B_{n}(x) d x=k^{n}}$

となる点です。

あるいは、この性質が成り立つような多項式としてBernoulli多項式が定められる、ともいえます。

この性質が本当に成り立っているのか、具体的に計算して確認してみましょう：

B1int_list = []
for k in range(7):
    b1int = sp.integrate(B_list[1], (x, k, k+1) )
    B1int_list.append(b1int)
B1int_list

より

$\quad\displaystyle{ \int_{k}^{k+1} B_{1}(x) d x= \{0, \quad 1,\quad 2,\quad 3,\quad 4,\quad 5, \quad 6,\ \dots\}, \qquad (k=0,1,2,3,\dots) }$

B3int_list = []
for k in range(7):
    b3int = sp.integrate(B_list[3], (x, k, k+1) )
    B3int_list.append(b3int)
B3int_list

より

$\quad\displaystyle{ \int_{k}^{k+1} B_{3}(x) d x= \{0, \quad 1,\quad 8,\quad 27,\quad 64,\quad 125, \quad 216, \ \dots \}, \qquad (k=0,1,2,3,\dots) }$

また

B5int_list = []
for k in range(7):
    b5int = sp.integrate(B_list[5], (x, k, k+1) )
    B5int_list.append(b5int)
B5int_list

より

$\quad\displaystyle{ \int_{k}^{k+1} B_{5}(x) d x= \{0, \quad 1,\quad 32,\quad 243,\quad 1024,\quad 3125, \quad 7776, \ \dots \}, \quad (k=0,1,2,3,\dots) }$

などと、確かに成り立っていることがわかります。

よって

$\quad\displaystyle{ \begin{aligned} S_n(N) =\sum_{k=1}^N k^{n} &=\sum_{k=1}^N \int_{k}^{k+1} B_{n}(x) d x\\ &=\left(\int_{1}^{2}+\int_{2}^{3}+\dots+\int_{N-1}^{N}+\int_{N}^{N+1}\right) B_{n}(x) d x =\int_{1}^{N+1} B_{n}(x) d x \end{aligned} }$

であることがわかります。実際

Bint_list = []
N=sp.symbols('N',integer=True)
for n in range(7):
    bint = sp.factor(sp.integrate(B_list[n], (x, 1, N+1) ))
    Bint_list.append(bint)
Bint_list

により

$\quad\displaystyle{ \begin{array}{l} S_{0}(N)=N \\ S_{1}(N)=\frac{N(N+1)}{2} \\ S_{2}(N)=\frac{N(N+1)(2 N+1)}{6} \\ S_{3}(N)=\frac{N^{2}(N+1)^{2}}{4}\\ S_{4}(N)=\frac{N(N+1)(2 N+1)\left(3 N^{2}+3 N-1\right)}{30}\\ \quad \vdots \end{array} }$

といった公式を得ることができ、本記事の冒頭に述べた公式を正しく与えています。

母関数

以上、具体的な計算から、Bernoulli多項式に関する公式が成り立つことを確認してみました。

しかしながら、それは数学的な証明ではないので、納得がいかないという人もいるかもしれません。

そこで、最後にBernoulli多項式の母関数

$\quad\displaystyle{ \frac{t\,e^{xt}}{e^t-1}=\sum_{n=0}^\infty B_n(x)\frac{t^n}{n!} }$

から、その性質を導いておこうと思います。まず、Bernoulli数 $b_k$ は

$\quad\displaystyle{ \frac{t}{e^t-1}=\sum_{k=0}^\infty b_k\frac{t^k}{k!} }$

で与えられます（むしろこの式でBernoulli数を定義しているといえます）。

母関数を２つの関数の積とみなし、それぞれをTaylor展開することによって

$\quad\displaystyle{ \begin{aligned} \frac{t\,e^{xt}}{e^t-1}=\frac{t}{e^t-1}\cdot e^{xt} &=\left(\sum_{k=0}^\infty b_k\frac{t^k}{k!}\right)\left(\sum_{l=0}^\infty \frac{x^lt^l}{l!}\right) =\sum_{k=0}^\infty\sum_{l=0}^\infty b_k\frac{x^l}{k!l!}t^{k+l}\\ &=\sum_{n=0}^\infty\sum_{k=0}^\infty b_kx^{n-k}\frac{n!}{k!(n-k)!} \frac{t^n}{n!} =\sum_{n=0}^\infty\left(\sum_{k=0}^\infty b_kx^{n-k}\binom{n}{k}\right)\frac{t^n}{n!}\\ &=\sum_{n=0}^\infty B_n(x)\frac{t^n}{n!} \end{aligned} }$

を得ます。実際

sp.var('t')
def G(t):
   return t*sp.exp(x*t)/(sp.exp(t)-1)
Gtayl_list = []
for n in range(7):
    g = sp.factorial(n-1)*sp.series(G(t), t, 0, n).coeff(t, n-1)
    Gtayl_list.append(g)
Gtayl_list

により、確かにBernoulli多項式の母関数になっていることが確かめられます。

母関数を積分することによって

$\quad\displaystyle{ \begin{aligned} \sum_{n=0}^\infty \frac{t^{n}}{n!}\int_k^{k+1}B_n(x) dx &=\int_k^{k+1} dx\frac{te^{xt}}{e^t-1} =\left[\frac{e^{xt}}{e^t-1}\right]_{x=k}^{k+1}\\ &=e^{kt}=\sum_{n=0}^\infty\frac{t^n}{n!}\,k^n \end{aligned} }$

より、公式

$\quad\displaystyle{ k^n=\int_{k}^{k+1} B_{n}(x) d x}$

を示すことができます。

また母関数の引数を $\ x\to x+1$ としたものと引き算して

$\quad\displaystyle{ \begin{aligned} \sum_{n=0}^\infty \left(B_n(x+1)-B_n(x)\right)\frac{t^n}{n!} &=\frac{t\,e^{(x+1)t}}{e^t-1}-\frac{t\,e^{xt}}{e^t-1}\\ &=t\,e^{xt}=\sum_{n=0}^\infty \frac{x^n\,t^{n+1}}{n!} \end{aligned} }$

により公式

$\quad\displaystyle{ B_n(x+1)-B_n(x)=nx^{n-1}}$

を得ます。実際

Bsub_list = []
for n in range(7):
    s = sp.factor(B_list[n].subs(x,x+1)-B_list[n])
    Bsub_list.append(s)
Bsub_list

により、確かに公式が成り立っています。

Bernoulli多項式は、その母関数も含めて、

とてもよくできた多項式だなという感動を共有できたら幸いです。

本記事で特に断りなく利用したSymPyの使い方については過去記事

pianofisica.hatenablog.com

が参考になるかと思います。あわせて読んでみてください。

キーワード：Bernoulli多項式、Python、SymPy