Python（SymPy）でFourier級数展開する

数学の具体的な計算にPython（SymPy）を使って、数学もPython（SymPy）も同時に学んでしまいましょう。今回はPython（SymPy）を使って関数のFourier（フーリエ）級数展開を見てみたいと思います。具体例として取り上げた関数について、そのFourier級数展開から円周率の級数表式が得られますが、おまけとしてこれをPythonで数値的に確かめてみます。

三角関数の直交性
Fourier級数展開
Fourier級数展開の例

三角関数の直交性

関数

$\qquad\displaystyle{\frac{1}{\sqrt{2\pi}},~\frac{\cos x}{\sqrt{\pi}},~\frac{\cos 2x}{\sqrt{\pi}},~\frac{\cos 3x}{\sqrt{\pi}},\, \dots,~\frac{\sin x}{\sqrt{\pi}},~\frac{\sin 2x}{\sqrt{\pi}},~\frac{\sin 3x}{\sqrt{\pi}},\,\dots}$

は $-\pi\leq x\leq \pi$ の変域で正規化された直交関数系をなします。実際に計算してみましょう：

import sympy as sp

x = sp.symbols('x')
def c(m, x): 
   if m == 0: return 1/sp.sqrt(2*sp.pi)
   else: return sp.cos(m*x)/sp.sqrt(sp.pi)
def s(m, x): return sp.sin(m*x)/sp.sqrt(sp.pi)

def cc(m,n): return sp.integrate(c(m,x)*c(n,x), (x, -sp.pi, sp.pi))
def cs(m,n): return sp.integrate(c(m,x)*s(n,x), (x, -sp.pi, sp.pi))
def ss(m,n): return sp.integrate(s(m,x)*s(n,x), (x, -sp.pi, sp.pi))

として

CC = [[],[],[],[],[]]
for m in range(5):
  for n in range(5):
    CC[m].append(cc(m,n))
print(CC)

より

$\qquad\displaystyle{\int^\pi_{-\pi} dx \frac{\cos mx}{\sqrt{\pi}} \frac{\cos nx}{\sqrt{\pi}}=\delta_{m,n} \qquad (m\neq0)}$

がわかります。同様に

SS = [[],[],[],[],[]]
for m in range(5):
  for n in range(5):
    SS[m].append(ss(m,n))
print(SS)

から

$\qquad\displaystyle{\int^\pi_{-\pi} dx \frac{\sin mx}{\sqrt{\pi}} \frac{\sin nx}{\sqrt{\pi}}=\delta_{m,n} \qquad (m\neq0)}$

CS = [[],[],[],[],[]]
for m in range(5):
  for n in range(5):
    CS[m].append(cs(m,n))
print(CS)

から

$\qquad\displaystyle{\int^\pi_{-\pi} dx \frac{\sin mx}{\sqrt{\pi}} \frac{\cos nx}{\sqrt{\pi}}=0}$

がわかります。以上より、うえで与えた関数系の正規直交性がわかります。

Fourier級数展開

変域 $-\pi\leq x\leq \pi$ で定義された関数 $f(x)$ の正規化された三角関数による展開式

$\qquad \displaystyle{\begin{aligned}f(x)&=A_0\frac{1}{\sqrt{2\pi}}+A_1\frac{\cos x}{\sqrt{\pi}}+B_1\frac{\sin x}{\sqrt{\pi}}+A_2\frac{\cos 2x}{\sqrt{\pi}}+B_2\frac{\sin 2x}{\sqrt{\pi}}+\dots\\ &=A_0\frac{1}{\sqrt{2\pi}}+\sum_{n=1}^\infty \left(A_n\frac{\cos nx}{\sqrt{\pi}}+B_n\frac{\sin nx}{\sqrt{\pi}}\right)\end{aligned}}$

をFourier級数展開といいます。

関数の正規直交性から、展開係数は

$\qquad \displaystyle{A_0=\int^\pi_{-\pi} dx \, \frac{1}{\sqrt{2\pi}} \, f(x),\quad A_n=\int^\pi_{-\pi} dx \, \frac{\cos nx}{\sqrt{\pi}} \,f(x), \quad B_n=\int^\pi_{-\pi} dx \, \frac{\sin nx}{\sqrt{\pi}}\, f(x)}$

によって求めることができます。

Fourier級数展開の例

ここでFourier級数の具体例を見てみましょう。

例は寺澤寛一『自然科学者のための数学概論』（岩波書店）を参照しました。

<br />

具体例：f(x)=x

$\qquad \displaystyle{f(x)=x \quad (-\pi\leq x\leq \pi)}$

sp.init_printing()
def f(x): return x
def a(n): return sp.integrate(c(n,x)*f(x), (x, -sp.pi, sp.pi))
def b(n): return sp.integrate(s(n,x)*f(x), (x, -sp.pi, sp.pi))
A = []
B = []
for n in range(10):
  A.append(a(n))
  B.append(b(n))
[A,B]

より

$\qquad\displaystyle{A_0=0, \quad A_n=0, \quad B_n=(-1)^{n+1}\sqrt{\pi}\,\frac{2}{n} \qquad (n\geq1)}$

と求まり、そのFourier級数展開が

$\qquad\displaystyle{x=2\left\{\frac{\sin x}{1}-\frac{\sin 2x}{2}+\frac{\sin 3x}{3}-\frac{\sin 4x}{4}+\dots\right\}}$

で与えられることがわかります。

実は、Fourier級数展開を求めるコマンドはPython（SymPy）に実装されています：

ff = sp.fourier_series(f(x), (x, -sp.pi, sp.pi)).truncate(8)  
# truncate( ) で表示する展開の次数を指定
ff

グラフの描画

$\ n=1,3,7$ までの展開がどれくらい関数を近似しているのかを実感するために、グラフを描いてみます：

import matplotlib.pyplot as plt
import numpy as np

D = 50
xmin = -np.pi
xmax = np.pi

def Ff(n, x):
   return sp.fourier_series(f(x), (x, -sp.pi, sp.pi)).truncate(n)  
p = np.linspace( xmin, xmax, D)
def FF(n,k): return float(Ff(n, x).subs(x, p[k]))
Q1 = []
Q3 = []
Q7 = []
for k in range(D):
   Q1.append(FF(1,k))
   Q3.append(FF(3,k))
   Q7.append(FF(7,k))
q = p
q1 = np.array(Q1)
q3 = np.array(Q3)
q7 = np.array(Q7)

plt.plot(p ,q)
plt.plot(p ,q1)
plt.plot(p ,q3)
plt.plot(p ,q7)
plt.show()

から得られるグラフが以下です：

f:id:pianofisica:20201005150546p:plain

青色の線がもとの関数、オレンジ色、緑色、赤色の順で、展開のより高次の振動モードを含めた場合になっています。

円周率（Leibniz／Eulerの級数）

さて、上で得られたFourier級数展開で $x=\pi/2$ とおくことで等式

$\qquad\displaystyle{\frac{\pi}{4}=1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\dots}$

がわかります。この級数はLeibnizの級数、またはEulerの級数と呼ばれます。

この級数を数値的に検証してみましょう：

k = sp.var('k', integer = True)
def p(n): return 4+4*sp.summation((-1)**k/(2*k+1), (k, 1, n))

とすることで、級数を $n$ 項まで足したときの値が数値的に求められます。

例えば1000項を足したときには

float(p(1000))

より、級数の値は

$\qquad 3.14259165\dots$

となります。円周率の値は

float(sp.pi)

より

$\qquad \pi=3.141592653589793\dots$

ですから、1000＝ $10^3$ 項の和で３桁程度までは数値的に再現できていることがわかります。

1万項の和だと

float(p(10000))

より（少し計算に時間がかかりますが）

$\qquad 3.14169264\dots$

となって精度が１桁ほど上がります。しかし計算量が多いようで、出力に少し時間がかかりました（もちろん、使っているコンピュータの性能によりますが）。あるいはリストを使えばより早く計算できます：

N = 10000
P = [4]
def q(k): return 4*(-1)**k/(2*k+1) + P[k-1]
for k in range(1,N):
  P.append(q(k))
P[N-1]

項の数を $10^{7}$ にとってみます：

N = 10**7
P = [4]
def q(k): return 4*(-1)**k/(2*k+1) + P[k-1]
for k in range(1,N):
  P.append(q(k))
P[N-1]

から（少し計算に時間がかかりますが）

$\qquad 3.14159255\dots$

と、徐々に精度が上がっていっていることがわかります。

キーワード：Python、SymPy、Fourier級数、直交関数