Return to SymPy で理工系の数学C

SymPy でフーリエ解析

必要なパッケージを import します。

In [1]:
from sympy.abc import * 
from sympy import *

# SymPy Plotting Backends (SPB): グラフを描く際に利用
from spb import *

# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']

フーリエ級数

  1. 区間 $-\pi < x < \pi$ で定義された関数 $f(x)$ は，それがどんな関数であっても…
  2. 区間外では，周期 $ 2\pi $ の周期関数とみなして
  3. 三角関数 $\cos, \ \sin$ の重ね合わせて表すことができる！

つまり，以下のように書けるということ。

$$ f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \bigl( a_n \cos n x + b_n \sin nx \bigr) $$

ここで，フーリエ係数 $a_n, b_n$ は

$$a_n = \frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \cos nx \, dx $$$$b_n = \frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \sin nx \, dx $$

例題

区間 $[-\pi: \pi]$ で定義された関数 $f(x) = x^2$ のフーリエ級数展開。

周期 $2 \pi$ の周期関数として，

  • $[-3\pi:-\pi]$ では $f(x) = (x+2 \pi)^2$…
  • $[-\pi:\pi]$ では $f(x) = x^2$,
  • $[\pi:3\pi]$ では $f(x) = (x-2 \pi)^2$…

のようにすればいい。

まず，1周期分の $f_0(x) \equiv x^2$ を定義して，$[-\pi:\pi]$ の区間でグラフを描く。

In [2]:
def f0(x):
    return x**2
In [3]:
plot(f0(x), (x, -pi, pi));

$[-\pi:\pi]$ の区間外では，周期 $2 \pi$ の周期関数として，

  • $[-3\pi:-\pi]$ では $f(x) = f_0(x+2 \pi)$…
  • $[-\pi:\pi]$ では $f(x) = f_0(x)$,
  • $[\pi:3\pi]$ では $f(x) = f_0(x-2 \pi)$…

をグラフに描く。

In [4]:
plot((f0(x+2*pi), (x, -3*pi, -pi)), 
     (f0(x), (x, -pi, pi)), 
     (f0(x-2*pi), (x, pi, 3*pi)), legend = False);

周期毎に別の関数を描くのは面倒なので，切り捨てる関数 floor() を使って周期 $2\pi$ の関数 $f(x)$ を定義します。

In [5]:
def f(x):
    return f0(x-2*pi*floor((x+pi)/(2*pi)))
In [6]:
plot(f(x), (x, -3*pi, 3*pi));

/usr/local/lib/python3.8/dist-packages/spb/series.py:587: UserWarning: NumPy is unable to evaluate with complex numbers some of the functions included in this symbolic expression: [floor]. Hence, the evaluation will use real numbers. If you believe the resulting plot is incorrect, change the evaluation module by setting the `modules` keyword argument.
  warnings.warn("NumPy is unable to evaluate with complex "

少しオプションを設定して描く例：

以下の例では，

  1. 凡例を設定し，
  2. 横軸の表示範囲を設定し，
  3. 縦軸の表示範囲を設定し，
  4. $x$ 軸のフォーマットを $\pi$ の倍数になるように設定しています。
In [7]:
fpi = float(pi)
p = plot(f(x), (x, -3*pi, 3*pi), 
         "$f(x)$", legend=True, 
         xlim = (-3*pi, 3*pi), 
         ylim = (-2, 12), show = False)

ax = p.ax
ax.set_xticks([(i*fpi) for i in range(-3,4)])
xlab = ["-3 $\pi$", "-2 $\pi$", "-$\pi$", "0", 
        "$\pi$", "$2 \pi$", "3 $\pi$"]
ax.set_xticklabels(xlab);

フーリエ級数は，以下のように書けました。

$$ f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \bigl( a_n \cos n x + b_n \sin nx \bigr) $$

ここで，フーリエ係数 $a_n, b_n$ は

$$a_n = \frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \cos nx \, dx $$$$b_n = \frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \sin nx \, dx $$

SymPy にはフーリエ級数を扱う関数がありますが，簡単なので，積分と和をとる演算を以下のように定義してしまいます。

In [8]:
# フーリエ係数の定義，フーリエ級数の定義。

n = Symbol('n', integer = True, positive = True)

def a(n):
    return 1/pi*integrate(f0(x)*cos(n*x), (x, -pi, pi))
def b(n):
    return 1/pi*integrate(f0(x)*sin(n*x), (x, -pi, pi))

def Fourier(n, x):
    var('i')
    return a(0)/2 + summation(
      a(i)*cos(i*x) + b(i)*sin(i*x), (i, 1, n))

ちなみに，上のセルの定義で使われている関数 integrate() は積分を実行し， summation() は和をとり，説明は以下の通りです。

$\displaystyle \int_a^b f(x)\, dx = $ integrate(f(x), x, a, b);

$\displaystyle \sum_{i = 1}^n a_i = $ summation(a(i), (i, 1, n));

$f(x) = x^2$ は偶関数なので，奇関数部分のフーリエ係数である $b_n$ はゼロのはず。実際に調べてみると…

In [9]:
for i in range(1, 5):
    print("b(%d) = " % i,b(i))

b(1) =  0
b(2) =  0
b(3) =  0
b(4) =  0
In [10]:
# 一般の n については...
b(n)
Out[10]:
$\displaystyle 0$
In [11]:
a(0)
Out[11]:
$\displaystyle \frac{2 \pi^{2}}{3}$
In [12]:
for i in range(1, 5):
    print("a(%d) = " % i,a(i))

a(1) =  -4
a(2) =  1
a(3) =  -4/9
a(4) =  1/4
In [13]:
# 一般の n については...
a(n)
Out[13]:
$\displaystyle \frac{4 \left(-1\right)^{n}}{n^{2}}$
In [14]:
f1 = Fourier(1, x)
f1
Out[14]:
$\displaystyle – 4 \cos{\left(x \right)} + \frac{\pi^{2}}{3}$
In [15]:
f2 = Fourier(2, x)
f2
Out[15]:
$\displaystyle – 4 \cos{\left(x \right)} + \cos{\left(2 x \right)} + \frac{\pi^{2}}{3}$
In [16]:
f3 = Fourier(3, x)
f3
Out[16]:
$\displaystyle – 4 \cos{\left(x \right)} + \cos{\left(2 x \right)} – \frac{4 \cos{\left(3 x \right)}}{9} + \frac{\pi^{2}}{3}$
In [17]:
f4 = Fourier(4, x)
f4
Out[17]:
$\displaystyle – 4 \cos{\left(x \right)} + \cos{\left(2 x \right)} – \frac{4 \cos{\left(3 x \right)}}{9} + \frac{\cos{\left(4 x \right)}}{4} + \frac{\pi^{2}}{3}$

それぞれの次数までのフーリエ級数展開の式と，もとの周期関数 $f(x)$ を重ねてグラフにします。

In [18]:
p=plot(
     (f(x), (x, -3*pi, 3*pi), "f(x)"), 
     (f4, (x, -3*pi, 3*pi), "n=4", {"linewidth":0.7}),
     (f3, (x, -3*pi, 3*pi), "n=3", {"linewidth":0.7}),
     (f2, (x, -3*pi, 3*pi), "n=2", {"linewidth":0.7}),
     (f1, (x, -3*pi, 3*pi), "n=1", {"linewidth":0.7}),
     xlim = (-3*pi, 3*pi), ylim = (-2, 12), 
     size=(6,6), show = False
)

ax = p.ax
ax.set_xticks([(i*fpi) for i in range(-3,4)])
xlab = ["-3 $\pi$", "-2 $\pi$", "-$\pi$", "0", 
        "$\pi$", "$2 \pi$", "3 $\pi$"]
ax.set_xticklabels(xlab);

任意の周期をもつ関数のフーリエ級数展開

周期 $2\pi$ の決め打ちのフーリエ級数ではなく，任意の周期をもつ場合は，一般に周期を $2L$ として…

$$ f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \bigl( a_n \cos \frac{n\pi x}{L} + b_n \sin \frac{n\pi x}{L} \bigr) $$$$a_n = \frac{1}{L} \int_{-L}^{L} f(x) \cos \frac{n\pi x}{L} \, d{x} $$$$b_n = \frac{1}{L} \int_{-L}^{L} f(x) \sin \frac{n\pi x}{L} \, d{x} $$

例題

区間 $[-1:1]$ で定義された関数 $f(x) = x$ が，区間外では周期 $2$ の周期関数であるとして，$n=5$ までのフーリエ級数展開を求める。

上記の式で $L = 1$ とすればよい。

一般に，1区間1周期で $f_0(x)$ として定義された関数を，その区間外で周期 $2L$ の周期関数にするには，上記のような $x_{\rm cyc}$ を使えばよい。

In [19]:
def xcyc(x, L):
    return x - 2*L*floor((x+L)/(2*L))
In [20]:
def f0(x):
    return x
def f(x):
    return f0(xcyc(x, 1))
In [21]:
plot(f(x), (x, -5, 5));

In [22]:
# L = 1 決めうちで
def a(n):
    return integrate(f0(x)*cos(n*pi*x), (x, -1, 1))
def b(n):
    return integrate(f0(x)*sin(n*pi*x), (x, -1, 1))

def Fourier(n, x):
    var('i')
    return a(0)/2 + summation(
      a(i)*cos(i*pi*x) + b(i)*sin(i*pi*x), (i, 1, n))

$f_0(x) = x$ は奇関数なので，偶関数部分のフーリエ係数である $a_n$ はゼロのはず。実際に調べてみると…

In [23]:
# 任意の n について...
a(n)
Out[23]:
$\displaystyle 0$
In [24]:
b(1)
Out[24]:
$\displaystyle \frac{2}{\pi}$
In [25]:
b(2)
Out[25]:
$\displaystyle – \frac{1}{\pi}$
In [26]:
b(3)
Out[26]:
$\displaystyle \frac{2}{3 \pi}$
In [27]:
b(n)
Out[27]:
$\displaystyle – \frac{2 \left(-1\right)^{n}}{\pi n}$
In [28]:
f2 = Fourier(2, x)
f2
Out[28]:
$\displaystyle – \frac{2 \left(- \sin{\left(\pi x \right)} + \frac{\sin{\left(2 \pi x \right)}}{2}\right)}{\pi}$
In [29]:
f3 = Fourier(3, x)
f3
Out[29]:
$\displaystyle – \frac{2 \left(- \sin{\left(\pi x \right)} + \frac{\sin{\left(2 \pi x \right)}}{2} – \frac{\sin{\left(3 \pi x \right)}}{3}\right)}{\pi}$
In [30]:
f4 = Fourier(4, x)
f4
Out[30]:
$\displaystyle – \frac{2 \left(- \sin{\left(\pi x \right)} + \frac{\sin{\left(2 \pi x \right)}}{2} – \frac{\sin{\left(3 \pi x \right)}}{3} + \frac{\sin{\left(4 \pi x \right)}}{4}\right)}{\pi}$
In [31]:
plot(
     (f(x), (x, -5, 5), "f(x)"), 
     (f4, (x, -5, 5), "n=4", {"linewidth":0.7}),
     (f3, (x, -5, 5), "n=3", {"linewidth":0.7}),
     (f2, (x, -5, 5), "n=2", {"linewidth":0.7}),
     xlim = (-5, 5), ylim = (-1.3, 1.8), 
     size = (6, 5)
);

フーリエ積分・フーリエ変換

任意の周期 $2L$ をもつ関数の複素フーリエ級数展開を，非周期的現象にまで拡張したものが「フーリエ積分」であり，フーリエ係数の拡張が「フーリエ変換」。

周期性のない関数に対する（連続極限としての）フーリエ積分

$$ f(x) = \frac{1}{2\pi} \int_{-\infty}^{\infty}\ F(k)\ e^{i k x}\ dk, \quad F(k) \equiv \int_{-\infty}^{\infty} f(x) \ e^{-i k x} \ dx $$

例題

以下の周期性のない関数 $f(x)$ のフーリエ変換 $F(k)$ を求めよ。

$$ f(x) = \begin{cases}
\frac{1}{a} & (|x| \leq \frac{a}{2}) \\
0 & (|x|> \frac{a}{2})
\end{cases}$$

次に，以下の極限を求めよ。$$ \lim_{a \rightarrow 0} F(k) =?$$

もし $\displaystyle |x| \leq \frac{a}{2}$ なら $\displaystyle \frac{1}{a}$ を返し，それ以外なら $0$ を返す関数 $f(x)$ の定義。

if 文を使いそうになるが，この場合は Piecewise() を使うらしい。

In [32]:
def f(a, x):
    return Piecewise((1/a, abs(x) <= a/2), 
                     (0, abs(x) > a/2))
In [33]:
plot(f(2, x), (x, -5, 5), ylim = (-0.5, 1));

$$ F(k) = \int_{-\infty}^{\infty} f(x) \ e^{-i k x} \ dx = \int_{-\frac{a}{2}}^{\frac{a}{2}} \frac{1}{a}\ e^{-i k x} \ dx$$ということを使って，フーリエ変換の無限積分を有限区間の積分にします。

In [34]:
var('a k', nonzero = True)

F = integrate(1/a*exp(-I*k*x), (x, -a/2, a/2))
F.simplify()
Out[34]:
$\displaystyle \frac{2 \sin{\left(\frac{a k}{2} \right)}}{a k}$
In [35]:
Eq(Limit(F.simplify(), a, 0), 
   Limit(F.simplify(), a, 0).doit())
Out[35]:
$\displaystyle \lim_{a \to 0^+}\left(\frac{2 \sin{\left(\frac{a k}{2} \right)}}{a k}\right) = 1$

算数の問題として，

$$ \lim_{a \rightarrow 0} F(k) = \lim_{a \rightarrow 0} \frac{\sin\frac{ak}{2}}{\frac{ak}{2}} = 1$$

となることはわかった。では，フーリエ変換 $F(k)$ が $1$ となる $f(x)$ とはどのような関数であろうか。

以下のグラフから推察されるように，$a$ の値を小さくしていくと $f(x)$ は $x = 0$ の近くでのみ値をもつ関数になる。

$\displaystyle \int_{-\infty}^{\infty} f(x)\, dx =\int_{-\frac{a}{2}}^{\frac{a}{2}} \frac{1}{a} \,dx = 1$ という面積を保ちながら，$a \rightarrow 0$ で幅 $a$ がどんどん狭くなっていくため，高さ $\displaystyle \frac{1}{a}$ がどんどん高くなっていきます。

In [36]:
plot((f(0.1, x), (x, -5, 5), "a=0.1"), 
     (f(0.2, x), (x, -5, 5), "a=0.2"), 
     (f(0.5, x), (x, -5, 5), "a=0.5"), 
     (f(2, x), (x, -5, 5), "a=2"), legend = True);

この関数 $f(x)$ の $a \rightarrow 0$ の極限は，「デルタ関数」$\delta(x)$ と呼ばれます。定義は，（フーリエ変換 $F(k)$ が $1$ であるフーリエ積分であるから）

$$\delta(x) \equiv \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{i k x}\ dk$$

デルタ関数は以下のような性質を持ちます。

$$\delta(x) = 0, \ \mbox{if $x \neq 0$}$$$$\int_{-\infty}^{\infty} g(x) \delta(x – c)\, dx = g(a)$$

特に，

$$\int_{-\infty}^{\infty} \delta(x)\, dx = 1$$

参考：fourier_series() 関数を使う場合

SymPy にはフーリエ級数展開する関数があります。これを使う例。

例題

区間 $[-\pi: \pi]$ で定義された関数 $f(x) = x^2$ のフーリエ級数展開。

まず，1周期分の関数 $f_0(x) = x^2$ を定義し…

In [37]:
f0 = x**2
s = fourier_series(f0, (x, -pi, pi))

$n=2$ までのフーリエ級数展開は n = 3 として…

n = 2 ではなく，n = 3 とする理由は $a_0, a_1, a_2$ と3個あるから？… というよりは，truncate(n = 3) の n = 3 がゼロではない最初の3項を表示するという意味があるから。

In [38]:
f2 = s.truncate(n = 3)
f2
Out[38]:
$\displaystyle – 4 \cos{\left(x \right)} + \cos{\left(2 x \right)} + \frac{\pi^{2}}{3}$

同様にして…

In [39]:
f3 = s.truncate(n = 4)
f3
Out[39]:
$\displaystyle – 4 \cos{\left(x \right)} + \cos{\left(2 x \right)} – \frac{4 \cos{\left(3 x \right)}}{9} + \frac{\pi^{2}}{3}$
In [40]:
f4 = s.truncate(n = 5)
f4
Out[40]:
$\displaystyle – 4 \cos{\left(x \right)} + \cos{\left(2 x \right)} – \frac{4 \cos{\left(3 x \right)}}{9} + \frac{\cos{\left(4 x \right)}}{4} + \frac{\pi^{2}}{3}$

例題

区間 $[-1:1]$ で定義された関数 $f(x) = x$ が，区間外では周期 $2$ の周期関数であるとして，$n=5$ までのフーリエ級数展開を求める。

In [41]:
f0 = x
s = fourier_series(f0, (x, -1, 1))
In [42]:
# この場合は n = 5 で 5次まで表示
f5 = s.truncate(n = 5)
f5
Out[42]:
$\displaystyle \frac{2 \sin{\left(\pi x \right)}}{\pi} – \frac{\sin{\left(2 \pi x \right)}}{\pi} + \frac{2 \sin{\left(3 \pi x \right)}}{3 \pi} – \frac{\sin{\left(4 \pi x \right)}}{2 \pi} + \frac{2 \sin{\left(5 \pi x \right)}}{5 \pi}$

SymPy でフーリエ級数展開の例題を解く