Return to SymPy でフーリエ解析

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']

init_printing();

mathtext.fontset の設定

In [2]:
# グラフを描くためではなくデフォルト設定のため
import matplotlib.pyplot as plt

plt.rcParams['mathtext.fontset'] = 'cm'

関数の定義

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

In [3]:
# -1 <= x <= 1 で定義された関数 f0(x)
def f0(x):
    return 1 - abs(x)

# f0(x) を任意周期 2 の周期関数にする小道具
def xcyc(x):
    return x - 2*floor((x+1)/2)

# 周期 2 の周期関数 f(x)
def f(x):
    return f0(xcyc(x))

関数のグラフ

区間 $-1 \le x \le 1$ で定義された関数 $f_0(x)$ と，周期 $2$ の周期関数 $f(x)$ をそれぞれグラフにしてみる。

In [4]:
plot(f0(x), (x, -1, 1), xlim=(-5, 5), ylabel="$f_0(x)$");
In [5]:
# plot(f(x), (x, -5, 5)) とするところだが
# plot(f, (x, -5, 5)) とすると警告が出ない。

plot(f, (x, -5, 5), xlim=(-5, 5), ylabel="$f(x)$");

フーリエ級数の定義

定義どおりにフーリエ級数展開を定義すると以下のようになる。

In [6]:
var('n', integer = True)

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):
    sums = a(0)/2 + summation(
                        a(i)*cos(i*pi*x) + 
                        b(i)*sin(i*pi*x), (i, 1, n))
    return sums

SymPy は優秀で，絶対値 abs(x) を含む積分でも特に問題なく積分値を返してくれる。たとえば，a(3) や Fourier(3) は…

In [7]:
a(3)
Out[7]:
$\displaystyle \frac{4}{9 \pi^{2}}$
In [8]:
Fourier(3)
Out[8]:
$\displaystyle \frac{4 \cos{\left(\pi x \right)}}{\pi^{2}} + \frac{4 \cos{\left(3 \pi x \right)}}{9 \pi^{2}} + \frac{1}{2}$

フーリエ係数

フーリエ係数 $a_n, b_n$ は…

In [9]:
a(n)
Out[9]:
$\displaystyle \begin{cases} – \frac{2 \left(-1\right)^{n}}{\pi^{2} n^{2}} + \frac{2}{\pi^{2} n^{2}} & \text{for}\: n \neq 0 \\1 & \text{otherwise} \end{cases}$
In [10]:
b(n)
Out[10]:
$\displaystyle 0$

$a_n$ は $n$ の偶奇によってゼロになる場合があるので，$n=10$ まで表示してみる。

In [11]:
for n in range(11):
    print('a(%2d) = ' % n, a(n))

a( 0) =  1
a( 1) =  4/pi**2
a( 2) =  0
a( 3) =  4/(9*pi**2)
a( 4) =  0
a( 5) =  4/(25*pi**2)
a( 6) =  0
a( 7) =  4/(49*pi**2)
a( 8) =  0
a( 9) =  4/(81*pi**2)
a(10) =  0

フーリエ級数のグラフ

$n = 3, 5, 7$ のフーリエ級数を求め，グラフにしてみる。

In [12]:
f3 = Fourier(3)
f5 = Fourier(5)
f7 = Fourier(7)

plot(
     (f, "$f(x)$"), 
     (f3, "$n=3$"), 
     (f5, "$n=5$"), 
     (f7, "$n=7$"), (x, -5, 5), {"lw":1}, 
     xlim = (-5, 6.5), ylabel = False);

$n=15$ くらいだと，かなり $f(x)$ に近い感じがあらわれている。

In [13]:
f15 = Fourier(15)

plot(
     (f, "$f(x)$"), 
     (f15, "$n=15$"), (x, -5, 5), {"lw":1},
     xlim = (-5, 6.5), ylabel = False);

参考：fourier_series() を使う場合

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

In [14]:
s = fourier_series(f0(x), (x, -1, 1))

ためしに，フーリエ級数の最初から（ゼロでない）5項を出力してみます。

In [15]:
s.truncate(n = 5)
Out[15]:
$\displaystyle \frac{4 \cos{\left(\pi x \right)}}{\pi^{2}} + \frac{4 \cos{\left(3 \pi x \right)}}{9 \pi^{2}} + \frac{4 \cos{\left(5 \pi x \right)}}{25 \pi^{2}} + \frac{4 \cos{\left(7 \pi x \right)}}{49 \pi^{2}} + \frac{1}{2}$

我々が定義した Fourier(n) では n=7 に相当します。（n が偶数の場合は a(n) がゼロであったことに注意。）

In [16]:
Fourier(7)
Out[16]:
$\displaystyle \frac{4 \cos{\left(\pi x \right)}}{\pi^{2}} + \frac{4 \cos{\left(3 \pi x \right)}}{9 \pi^{2}} + \frac{4 \cos{\left(5 \pi x \right)}}{25 \pi^{2}} + \frac{4 \cos{\left(7 \pi x \right)}}{49 \pi^{2}} + \frac{1}{2}$