Matplotlib でフーリエ級数展開の説明用の図を描く

フーリエ級数展開の説明用の図を Matplotlib で描く。以下のページで使っているので。

In [1]:
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np

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

plt.rcParams['font.family'] = 'serif'
plt.rcParams['mathtext.fontset'] = 'cm'

周期 $2 \pi$ のフーリエ級数展開の例

$ -\pi \leq x \leq \pi$ で定義された関数 $f(x) = x^2$

In [2]:
def f(x):
    return x**2
In [3]:
fig = plt.figure(figsize=[6.4, 4.8])
fig.tight_layout()

# f(x)
x = np.linspace(-np.pi, np.pi, 200)
plt.plot(x, f(x), lw = 2)

# 横軸縦軸の表示範囲
plt.xlim(-3*np.pi, 3*np.pi)
plt.ylim(-1, 12)

# グリッド（格子線）
plt.grid(ls=':')

# 軸ラベル，タイトル
plt.xlabel('$x$', fontsize=12)
plt.title('$f(x) = x^2 \ \ (-\pi \leq x \leq \pi)$', fontsize=14)

plt.xticks([i*np.pi for i in range(-3,4)], 
           ['$-3\pi$','$-2 \pi$','$-\pi$','$0$',
            '$\pi$',  '$2\pi$', '$3\pi$'], fontsize=12)
# x軸 y軸
plt.axhline(0, color='gray', ls = '--', lw=1)
plt.axvline(0, color='gray', ls = '--', lw=1);

区間外では周期 $2\pi$ の周期関数に

In [4]:
def xcyc(x, L=np.pi):
    return x - 2*L*np.floor((x+L)/(2*L))
In [5]:
fig = plt.figure(figsize=[6.4, 4.8])
fig.tight_layout()

# f(x)
x = np.linspace(-3*np.pi, 3*np.pi, 600)
# 指定しなければ L = np.pi
plt.plot(x, f(xcyc(x)), lw = 2)

# 横軸縦軸の表示範囲
plt.xlim(-3*np.pi, 3*np.pi)
plt.ylim(-1, 12)

# グリッド（格子線）
plt.grid(ls=':')

# 軸ラベル，タイトル
plt.xlabel('$x$', fontsize=12)
plt.title('周期 $2 \pi$ の周期関数 $f(x)$', fontsize=14)

plt.xticks([i*np.pi for i in range(-3,4)], 
           ['$-3\pi$','$-2 \pi$','$-\pi$','$0$',
            '$\pi$',  '$2\pi$', '$3\pi$'], fontsize=12)
# x軸 y軸
plt.axhline(0, color='gray', ls = '--', lw=1)
plt.axvline(0, color='gray', ls = '--', lw=1);

フーリエ係数，フーリエ級数

In [6]:
def a(n):
    return 4*(-1)**n/n**2

def Fourier(n, x):
    if n ==1:
        return np.pi**2/3 + a(1)*np.cos(1*x)
    else:
        return Fourier(n-1, x) + a(n)*np.cos(n*x)

$n = 10$ のグラフ

In [7]:
fig = plt.figure(figsize=[6.4, 4.8])
fig.tight_layout()

# f(x)
x = np.linspace(-3*np.pi, 3*np.pi, 600)
plt.plot(x, f(xcyc(x)), label="$f(x)$", lw = 2)

# Fourier(n, x)
n = 10
Label = '$n='+str(n)+'$'
Title = '$f(x)$ のフーリエ級数展開 $n='+str(n)+'$'
plt.plot(x, Fourier(n, x), label=Label, lw = 1)
plt.title(Title)

# 横軸縦軸の表示範囲
plt.xlim(-3*np.pi, 3*np.pi)
plt.ylim(-1, 12)

# グリッド（格子線）
plt.grid(ls=':')

# 軸ラベル，タイトル
plt.xlabel('$x$', fontsize=12)
plt.title('周期 $2 \pi$ の周期関数 $f(x)$', fontsize=14)

plt.xticks([i*np.pi for i in range(-3,4)], 
           ['$-3\pi$','$-2 \pi$','$-\pi$','$0$',
            '$\pi$',  '$2\pi$', '$3\pi$'], fontsize=12)
# x軸 y軸
plt.axhline(0, color='gray', ls = '--', lw=1)
plt.axvline(0, color='gray', ls = '--', lw=1);
plt.legend(fontsize=12);

filename = 'fourier-fig'+str(n).zfill(2)+'.svg'
plt.savefig(filename);

$n=1$ から $n=5$ までのグラフ

In [8]:
for n in range(1, 6):
    fig = plt.figure(figsize=[6.4, 4.8])
    fig.tight_layout()

    # 横軸縦軸の表示範囲
    plt.xlim(-3*np.pi, 3*np.pi)
    plt.ylim(-1, 12)

    # グリッド（格子線）
    plt.grid(ls=':')

    # 軸ラベル，タイトル
    plt.xlabel('$x$', fontsize=12)
    plt.title('周期 $2 \pi$ の周期関数 $f(x)$', fontsize=14)

    plt.xticks([i*np.pi for i in range(-3,4)], 
               ['$-3\pi$','$-2 \pi$', '$-\pi$', '$0$',
                '$\pi$',  '$2\pi$',   '$3\pi$'], fontsize=12)
    # x軸 y軸
    plt.axhline(0, color='gray', ls = '--', lw=1)
    plt.axvline(0, color='gray', ls = '--', lw=1);

    # f(x)
    x = np.linspace(-3*np.pi, 3*np.pi, 600)
    plt.plot(x, f(xcyc(x)), label="$f(x)$", lw = 2)

    # Fourier(n, x)
    Label = '$n='+str(n)+'$'
    Title = '$f(x)$ のフーリエ級数展開 $n='+str(n)+'$'
    plt.plot(x, Fourier(n, x), label=Label, lw = 1)
    plt.title(Title)
    plt.legend(fontsize=12);

    filename = 'fourier-fig'+str(n).zfill(2)+'.svg'
    plt.savefig(filename);

アニメーション

In [9]:
from matplotlib.animation import FuncAnimation
In [10]:
# 小綺麗な動画にするために解像度を設定
fig = plt.figure(dpi=288)
fig.tight_layout()

# n=i+1 のフーリエ級数展開を描くように
# func を定義
def func(i):
    # 前の frame を消す
    plt.cla()
    # 横軸縦軸の表示範囲
    plt.xlim(-3*np.pi, 3*np.pi)
    plt.ylim(-1, 12)

    # グリッド（格子線）
    plt.grid(ls=':')

    # 軸ラベル，タイトル
    plt.xlabel('$x$', fontsize=12)
    plt.title('周期 $2 \pi$ の周期関数 $f(x)$', fontsize=14)

    plt.xticks([i*np.pi for i in range(-3,4)], 
               ['$-3\pi$','$-2 \pi$', '$-\pi$', '$0$',
                '$\pi$',  '$2\pi$',   '$3\pi$'], fontsize=12)
    # x軸 y軸
    plt.axhline(0, color='gray', ls = '--', lw=1)
    plt.axvline(0, color='gray', ls = '--', lw=1);

    # f(x)
    x = np.linspace(-3*np.pi, 3*np.pi, 600)
    plt.plot(x, f(xcyc(x)), label="$f(x)$", lw = 2)

    # Fourier(n, x)
    Label = '$n='+str(i+1)+'$'
    Title = '$f(x)$ のフーリエ級数展開 $n='+str(i+1)+'$'
    plt.plot(x, Fourier(i+1, x), label=Label, lw = 1)
    plt.title(Title)
    plt.legend(fontsize=12)

# 変数名 frames は固定。
frames = 10
ani = FuncAnimation(fig, func,  
        # interval は frame 間の時間をミリ秒単位で。
        interval = 1500, 
        frames = range(frames))

# 動画を jupyterhub のホームに mp4 ファイルとして保存。
ani.save("Fourier-anim.mp4")

任意周期のフーリエ級数展開の例

$-1 \leq x \leq 1$ で定義された関数 $f(x) = x$

In [11]:
def f(x):
    return x
In [12]:
fig = plt.figure(figsize=[6.4, 4.8])
fig.tight_layout()

# f(x)
x = np.linspace(-1, 1, 200)
plt.plot(x, f(x), lw = 2)

# 横軸縦軸の表示範囲
plt.xlim(-5, 5)
plt.ylim(-1.2, 1.7)

# グリッド（格子線）
plt.grid(ls=':')

# 軸ラベル，タイトル
plt.xlabel('$x$', fontsize=12)
plt.title('$f(x) = x \ \ (-1 \leq x \leq 1)$', fontsize=14)

plt.xticks([i for i in range(-5,6)], fontsize=12)
# x軸 y軸
plt.axhline(0, color='gray', ls = '--', lw=1)
plt.axvline(0, color='gray', ls = '--', lw=1);

plt.savefig('fourier-Fx3.svg');

区間外では周期 $2$ の周期関数に

In [13]:
fig = plt.figure(figsize=[6.4, 4.8])
fig.tight_layout()

# f(x)
x = np.arange(-5, 5, 1/100)
# L = 1
plt.plot(x, f(xcyc(x, 1)), lw = 2)

# 横軸縦軸の表示範囲
plt.xlim(-5, 5)
plt.ylim(-1.2, 1.7)

# グリッド（格子線）
plt.grid(ls=':')

# 軸ラベル，タイトル
plt.xlabel('$x$', fontsize=12)
plt.title('周期 $2$ の周期関数 $f(x)$', fontsize=14)

plt.xticks([i for i in range(-5,6)], fontsize=12)
# x軸 y軸
plt.axhline(0, color='gray', ls = '--', lw=1)
plt.axvline(0, color='gray', ls = '--', lw=1);

フーリエ係数，フーリエ級数

In [14]:
def b(n):
    return 2*(-1)**(n+1)/(n*np.pi)

def Fourier(n, x):
    if n ==1:
        return b(1)*np.sin(np.pi*x)
    else:
        return Fourier(n-1, x) + b(n)*np.sin(n*np.pi*x)

$n=1$ から $n=3$ までのグラフ

In [15]:
fig = plt.figure(figsize=[6.4, 4.8])
fig.tight_layout()

# f(x)
x = np.arange(-5, 5, 1/100)
# L = 1
plt.plot(x, f(xcyc(x, 1)), lw = 2, label='$f(x)$')

for n in range(1, 4):
    Label = '$n='+str(n)+"$"
    plt.plot(x, Fourier(n, x), lw = 1, label=Label)

# 横軸縦軸の表示範囲
plt.xlim(-5, 5)
plt.ylim(-1.2, 1.7)

# グリッド（格子線）
plt.grid(ls=':')

# 軸ラベル，タイトル
plt.xlabel('$x$', fontsize=12)
plt.title('$f(x)$ のフーリエ級数展開', fontsize=14)

plt.xticks([i for i in range(-5,6)], fontsize=12)
# x軸 y軸
plt.axhline(0, color='gray', ls = '--', lw=1)
plt.axvline(0, color='gray', ls = '--', lw=1);

plt.legend(fontsize=10);

$n=10$ のグラフ

In [16]:
fig = plt.figure(figsize=[6.4, 4.8])
fig.tight_layout()

# f(x)
x = np.arange(-5, 5, 1/100)
# L = 1
plt.plot(x, f(xcyc(x, 1)), lw = 2, label='$f(x)$')

n = 10
Label = '$n='+str(n)+"$"
plt.plot(x, Fourier(n, x), lw = 1, label=Label)

# 横軸縦軸の表示範囲
plt.xlim(-5, 5)
plt.ylim(-1.2, 1.7)

# グリッド（格子線）
plt.grid(ls=':')

# 軸ラベル，タイトル
plt.xlabel('$x$', fontsize=12)
plt.title('$f(x)$ のフーリエ級数展開 $n=10$', fontsize=14)

plt.xticks([i for i in range(-5,6)], fontsize=12)
# x軸 y軸
plt.axhline(0, color='gray', ls = '--', lw=1)
plt.axvline(0, color='gray', ls = '--', lw=1);

plt.legend(fontsize=12);

アニメーション

In [17]:
# 小綺麗な動画にするために解像度を設定
fig = plt.figure(dpi=288)
fig.tight_layout()

cmap = plt.get_cmap("tab10")

# n=i+1 のフーリエ級数展開を描くように
# func を定義
def func(i):
    # 前の frame を消す
    plt.cla()
    # f(x)
    x = np.arange(-5, 5, 1/100)
    # L = 1
    plt.plot(x, f(xcyc(x, 1)), lw = 2, label='$f(x)$')

    Label1 = '$n='+str(i+1)+"$"
    plt.plot(x, Fourier(i+1, x), lw = 1, label=Label1, c=cmap(i))
    Label2 = '$n='+str(i+2)+"$"
    plt.plot(x, Fourier(i+2, x), lw = 1, label=Label2, c=cmap(i+1))

    # 横軸縦軸の表示範囲
    plt.xlim(-5, 5)
    plt.ylim(-1.2, 1.7)

    # グリッド（格子線）
    plt.grid(ls=':')

    # 軸ラベル，タイトル
    plt.xlabel('$x$', fontsize=12)
    plt.title('$f(x)$ のフーリエ級数展開', fontsize=14)

    plt.xticks([i for i in range(-5,6)], fontsize=12)
    # x軸 y軸
    plt.axhline(0, color='gray', ls = '--', lw=1)
    plt.axvline(0, color='gray', ls = '--', lw=1);

    plt.legend(fontsize=10)

# 変数名 frames は固定。
frames = 9
ani = FuncAnimation(fig, func,  
        # interval は frame 間の時間をミリ秒単位で。
        interval = 1500, 
        frames = range(frames))

# 動画を jupyterhub のホームに mp4 ファイルとして保存。
ani.save("fourier-Anim.mp4")