工業数学完全攻略2: フーリエ級数・変換と偏微分方程式（PDE）

信号処理、通信システム、電磁気学の数学的基盤は、まさに**フーリエ解析（Fourier Analysis）**です。どんな周期信号もサイン・コサインの和で表現できるという驚くべき事実が、現代工学の核心です。この記事ではフーリエ級数からFFT、偏微分方程式まで体系的にまとめます。

1. フーリエ級数

1.1 核心となるアイデア

1807年、ジョゼフ・フーリエ（Joseph Fourier）は熱方程式を解く中で驚くべき主張をしました: どんな周期関数も三角関数の無限級数で表現できるということです。

周期 $2L$ の関数 $f(x)$ のフーリエ級数：

$f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty}\left(a_n\cos\frac{n\pi x}{L} + b_n\sin\frac{n\pi x}{L}\right)$

フーリエ係数：

$a_0 = \frac{1}{L}\int_{-L}^{L}f(x)\,dx$

$a_n = \frac{1}{L}\int_{-L}^{L}f(x)\cos\frac{n\pi x}{L}\,dx \quad (n = 1, 2, 3, \ldots)$

$b_n = \frac{1}{L}\int_{-L}^{L}f(x)\sin\frac{n\pi x}{L}\,dx \quad (n = 1, 2, 3, \ldots)$

直交性（Orthogonality）: この公式が成立する理由は三角関数の直交性によります。

$\int_{-L}^{L}\cos\frac{m\pi x}{L}\cos\frac{n\pi x}{L}\,dx = \begin{cases} 2L & m = n = 0 \\ L & m = n \neq 0 \\ 0 & m \neq n \end{cases}$

$\int_{-L}^{L}\sin\frac{m\pi x}{L}\sin\frac{n\pi x}{L}\,dx = \begin{cases} L & m = n \\ 0 & m \neq n \end{cases}$

$\int_{-L}^{L}\cos\frac{m\pi x}{L}\sin\frac{n\pi x}{L}\,dx = 0 \quad (\text{常に})$

1.2 方形波の展開

問題: 周期 $2\pi$ の方形波

$f(x) = \begin{cases} 1 & 0 < x < \pi \\ -1 & -\pi < x < 0 \end{cases}$

計算:

$f(x)$ は奇関数なので $a_n = 0$ ($n \geq 0$)。

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

$= \frac{2}{\pi}\left[-\frac{\cos(nx)}{n}\right]_0^{\pi} = \frac{2}{n\pi}(1 - \cos n\pi) = \begin{cases} \frac{4}{n\pi} & n \text{ 奇数} \\ 0 & n \text{ 偶数} \end{cases}$

結果: $f(x) = \frac{4}{\pi}\left(\sin x + \frac{1}{3}\sin 3x + \frac{1}{5}\sin 5x + \cdots\right) = \frac{4}{\pi}\sum_{k=0}^{\infty}\frac{\sin(2k+1)x}{2k+1}$

最初の数項だけ足しても方形波に次第に近づきます。これがフーリエ級数の威力です。

ライプニッツの公式: $x = \pi/2$ を代入すると： $\frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \cdots$

1.3 三角波の展開

$f(x) = |x|, \quad -\pi \leq x \leq \pi$

偶関数なので $b_n = 0$。

$a_0 = \frac{2}{\pi}\int_0^{\pi} x\,dx = \pi$

$a_n = \frac{2}{\pi}\int_0^{\pi} x\cos(nx)\,dx = \frac{2}{\pi}\left[\frac{x\sin(nx)}{n} + \frac{\cos(nx)}{n^2}\right]_0^{\pi}$

$= \frac{2}{n^2\pi}(\cos n\pi - 1) = \begin{cases} -\frac{4}{n^2\pi} & n \text{ 奇数} \\ 0 & n \text{ 偶数} \end{cases}$

結果: $f(x) = \frac{\pi}{2} - \frac{4}{\pi}\left(\cos x + \frac{\cos 3x}{9} + \frac{\cos 5x}{25} + \cdots\right)$

方形波と違い、係数が $n^2$ で速く減少するため、より速く収束します。

1.4 収束とギブス現象

ディリクレ条件（Dirichlet Conditions）:

関数 $f(x)$ が以下の条件を満たせばフーリエ級数は収束します：

1. 1周期内で絶対積分可能: $\int_{-L}^{L}|f(x)|\,dx < \infty$
2. 1周期内での極値の数が有限
3. 1周期内での不連続点の数が有限

不連続点では、級数は左極限と右極限の平均に収束します： $f(x_0) = \frac{f(x_0^-) + f(x_0^+)}{2}$

ギブス現象: 不連続点近くで有限個の項の和にオーバーシュートが発生します。項数をいくら増やしても不連続点近くで約**8.9%**のオーバーシュートが消えません。

$\text{オーバーシュート} \approx \frac{1}{\pi}\int_0^{\pi}\frac{\sin t}{t}\,dt - \frac{1}{2} \approx 0.0895$

これが方形波フィルタのリンギング（ringing）現象の数学的原因です。

1.5 複素フーリエ級数

オイラーの公式を利用すると：

$f(x) = \sum_{n=-\infty}^{\infty}c_n e^{in\pi x/L}$

複素係数： $c_n = \frac{1}{2L}\int_{-L}^{L}f(x)e^{-in\pi x/L}\,dx$

$a_n$、$b_n$ との関係： $c_n = \frac{a_n - ib_n}{2}, \quad c_{-n} = \frac{a_n + ib_n}{2} \quad (n > 0)$

パーセバルの定理（エネルギー保存）:

$\frac{1}{2L}\int_{-L}^{L}|f(x)|^2\,dx = \sum_{n=-\infty}^{\infty}|c_n|^2 = \frac{|a_0|^2}{4} + \frac{1}{2}\sum_{n=1}^{\infty}(|a_n|^2 + |b_n|^2)$

時間領域のエネルギーと周波数領域のエネルギーが等しいという意味です。

バーゼル問題への応用: 三角波のパーセバル定理を適用すると： $\sum_{n=1,3,5,\ldots}\frac{1}{n^4} = \frac{\pi^4}{96}$

2. フーリエ変換

2.1 周期関数から非周期関数へ

周期 $2L$ を $L \to \infty$ に持っていくと、フーリエ級数がフーリエ変換に発展します。

フーリエ変換（時間 -> 周波数）: $F(\omega) = \mathcal{F}\{f(t)\} = \int_{-\infty}^{\infty}f(t)e^{-i\omega t}\,dt$

逆フーリエ変換（周波数 -> 時間）: $f(t) = \mathcal{F}^{-1}\{F(\omega)\} = \frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)e^{i\omega t}\,d\omega$

周波数 $f$ を使う表現（工学でよく使用）: $F(f) = \int_{-\infty}^{\infty}f(t)e^{-i2\pi ft}\,dt, \quad f(t) = \int_{-\infty}^{\infty}F(f)e^{i2\pi ft}\,df$

2.2 主要変換対

矩形パルス:

$f(t) = \begin{cases} 1 & |t| < T/2 \\ 0 & |t| > T/2 \end{cases}$

$F(\omega) = \int_{-T/2}^{T/2}e^{-i\omega t}\,dt = \frac{2\sin(\omega T/2)}{\omega} = T\,\text{sinc}\!\left(\frac{\omega T}{2\pi}\right)$

時間領域の有限なパルスが周波数領域でsinc関数になります。

不確定性原理: 時間幅 $\Delta t$ と帯域幅 $\Delta\omega$ の積には下限があります。 $\Delta t \cdot \Delta\omega \geq \frac{1}{2}$

パルスが短いほど（時間分解能が高い）帯域幅が広くなります（周波数分解能が低い）。

ガウス関数: 自己双対変換対

$f(t) = e^{-at^2} \Leftrightarrow F(\omega) = \sqrt{\frac{\pi}{a}}e^{-\omega^2/(4a)}$

ガウス関数のフーリエ変換は再びガウス関数です。

ディラックのデルタ関数:

$\delta(t) = 0 \,(t \neq 0), \quad \int_{-\infty}^{\infty}\delta(t)\,dt = 1$

$\mathcal{F}\{\delta(t)\} = 1$

インパルスの周波数スペクトルはすべての周波数で一様です（白色雑音と類似）。

$\mathcal{F}\{1\} = 2\pi\delta(\omega)$

純粋なDC信号は $\omega = 0$ でのみエネルギーを持ちます。

2.3 主要性質

線形性: $\mathcal{F}\{af(t) + bg(t)\} = aF(\omega) + bG(\omega)$

時間移動: $\mathcal{F}\{f(t-t_0)\} = e^{-i\omega t_0}F(\omega)$

遅延は振幅を変えず位相のみを変化させます。

周波数移動（変調）: $\mathcal{F}\{e^{i\omega_0 t}f(t)\} = F(\omega - \omega_0)$

時間スケーリング: $\mathcal{F}\{f(at)\} = \frac{1}{|a|}F\!\left(\frac{\omega}{a}\right)$

信号を時間圧縮すると帯域幅が広がります。

微分: $\mathcal{F}\{f^{(n)}(t)\} = (i\omega)^n F(\omega)$

微分が周波数領域で乗算に変換されます。ODEが代数方程式に！

積分: $\mathcal{F}\left\{\int_{-\infty}^{t}f(\tau)\,d\tau\right\} = \frac{F(\omega)}{i\omega} + \pi F(0)\delta(\omega)$

畳み込み定理（信号処理の核心）: $\mathcal{F}\{(f*g)(t)\} = F(\omega) \cdot G(\omega)$

$\text{ここで } (f*g)(t) = \int_{-\infty}^{\infty}f(\tau)g(t-\tau)\,d\tau$

線形時不変（LTI）システムでは：出力 = 入力とインパルス応答の畳み込み → 周波数領域では乗算！

パーセバルの定理（エネルギー保存）: $\int_{-\infty}^{\infty}|f(t)|^2\,dt = \frac{1}{2\pi}\int_{-\infty}^{\infty}|F(\omega)|^2\,d\omega$

2.4 応用：フィルタ設計

理想的な低域通過フィルタ（LPF）:

$H(\omega) = \begin{cases} 1 & |\omega| \leq \omega_c \\ 0 & |\omega| > \omega_c \end{cases}$

インパルス応答： $h(t) = \mathcal{F}^{-1}\{H(\omega)\} = \frac{\omega_c}{\pi}\text{sinc}\!\left(\frac{\omega_c t}{\pi}\right)$

因果性がない（未来のサンプルが必要）ため実際には実装不可 → 実用フィルタ（バターワース、チェビシェフなど）を使用。

帯域通過フィルタ:

$H(\omega) = H_{LPF}(\omega - \omega_0) + H_{LPF}(\omega + \omega_0)$

インパルス応答: $h(t) = 2h_{LPF}(t)\cos\omega_0 t$

周波数応答:

LTIシステムに $e^{i\omega t}$ を入力すると： $y(t) = H(\omega)e^{i\omega t}$

$H(\omega)$ は複素ゲインと位相ずれを同時に表します： $H(\omega) = |H(\omega)|e^{i\angle H(\omega)}$

3. 離散フーリエ変換（DFT）とFFT

3.1 DFTの定義

$N$ 個の離散サンプル $x[0], x[1], \ldots, x[N-1]$ に対して：

$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/N}, \quad k = 0, 1, \ldots, N-1$

$x[n] = \frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi kn/N}, \quad n = 0, 1, \ldots, N-1$

回転因子（twiddle factor）: $W_N = e^{-j2\pi/N}$

$X[k] = \sum_{n=0}^{N-1}x[n]W_N^{kn}$

周波数分解能: $\Delta f = f_s / N$ (サンプリング周波数をNで割る)

ナイキスト限界: $f_s \geq 2 f_{max}$ (エイリアシング防止)

3.2 FFTアルゴリズム（クーリー-テューキーアルゴリズム）

直接DFT計算: $O(N^2)$ 計算量

FFT（高速フーリエ変換）: $O(N\log_2 N)$ 計算量

$N = 2^{10} = 1024$ の場合：

• DFT: $1024^2 \approx 10^6$ 演算
• FFT: $1024 \times 10 \approx 10^4$ 演算 → 100倍速い

分割統治原理（時間間引き）:

偶数・奇数インデックスに分離：

$X[k] = \sum_{n=0}^{N/2-1}x[2n]W_N^{2kn} + W_N^k\sum_{n=0}^{N/2-1}x[2n+1]W_N^{2kn}$

$= E[k] + W_N^k O[k]$

$W_N^2 = W_{N/2}$ なので $E[k]$ と $O[k]$ はそれぞれ $N/2$ 点DFTです。

これを再帰的に適用すると $O(N\log N)$ になります。

3.3 PythonでFFTを実装して信号を解析する

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as sig

# 複合信号の生成（50 Hz + 120 Hz 混合）
fs = 1000  # サンプリング周波数 (Hz)
T = 1.0    # 信号長 (秒)
N = int(fs * T)
t = np.linspace(0, T, N, endpoint=False)

# 信号: 50 Hz (振幅 1.0) + 120 Hz (振幅 0.5) + ノイズ
np.random.seed(42)
signal_clean = (np.sin(2 * np.pi * 50 * t) +
                0.5 * np.sin(2 * np.pi * 120 * t))
noise = 0.2 * np.random.randn(N)
signal_noisy = signal_clean + noise

# FFT実行
fft_result = np.fft.fft(signal_noisy)
freqs = np.fft.fftfreq(N, 1/fs)

# 両側スペクトル -> 片側スペクトル
magnitude = (2.0 / N) * np.abs(fft_result[:N//2])
freqs_pos = freqs[:N//2]

# 可視化
fig, axes = plt.subplots(3, 1, figsize=(12, 10))

axes[0].plot(t[:200], signal_noisy[:200], 'b-', alpha=0.7)
axes[0].plot(t[:200], signal_clean[:200], 'r-', linewidth=2)
axes[0].set_xlabel('時間 (s)')
axes[0].set_ylabel('振幅')
axes[0].set_title('時間領域信号 (青: ノイズあり、赤: 元信号)')
axes[0].legend(['ノイズあり', '元信号'])
axes[0].grid(True, alpha=0.3)

axes[1].plot(freqs_pos, magnitude, 'g-')
axes[1].set_xlabel('周波数 (Hz)')
axes[1].set_ylabel('振幅スペクトル')
axes[1].set_title('周波数領域 (FFT 振幅)')
axes[1].set_xlim([0, 200])
axes[1].grid(True, alpha=0.3)

# スペクトログラム
f_spec, t_spec, Sxx = sig.spectrogram(signal_noisy, fs, nperseg=64)
axes[2].pcolormesh(t_spec, f_spec, 10*np.log10(Sxx + 1e-10), shading='gouraud')
axes[2].set_ylabel('周波数 (Hz)')
axes[2].set_xlabel('時間 (s)')
axes[2].set_title('スペクトログラム')
axes[2].set_ylim([0, 200])

plt.tight_layout()
plt.savefig('fft_analysis.png', dpi=150)
plt.show()

# 主要周波数成分の出力
peak_indices = np.where(magnitude > 0.1)[0]
print("主要周波数成分:")
for idx in peak_indices:
    print(f"  {freqs_pos[idx]:.1f} Hz: 振幅 {magnitude[idx]:.3f}")

3.4 FFTを使った信号フィルタリング

import numpy as np
import matplotlib.pyplot as plt

def fft_lowpass_filter(signal_in, fs, cutoff_freq):
    """FFTベースの理想低域通過フィルタ"""
    N = len(signal_in)
    fft_result = np.fft.fft(signal_in)
    freqs = np.fft.fftfreq(N, 1/fs)

    # 遮断周波数以上を除去
    fft_filtered = fft_result.copy()
    fft_filtered[np.abs(freqs) > cutoff_freq] = 0

    # 逆変換
    return np.real(np.fft.ifft(fft_filtered))

# フィルタリング例
fs = 1000
t = np.linspace(0, 1, fs, endpoint=False)
signal_in = (np.sin(2 * np.pi * 50 * t) +
             0.3 * np.sin(2 * np.pi * 200 * t) +
             0.1 * np.random.randn(fs))

filtered = fft_lowpass_filter(signal_in, fs, cutoff_freq=80)

plt.figure(figsize=(10, 5))
plt.plot(t[:200], signal_in[:200], 'b-', alpha=0.5, label='元信号')
plt.plot(t[:200], filtered[:200], 'r-', linewidth=2, label='フィルタ後の信号 (80Hz LPF)')
plt.xlabel('時間 (s)')
plt.ylabel('振幅')
plt.title('FFTベースの低域通過フィルタリング')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

4. 偏微分方程式（PDE）

4.1 PDEの分類

2階線形PDE： $Au_{xx} + Bu_{xy} + Cu_{yy} + Du_x + Eu_y + Fu = G$

判別式 $B^2 - 4AC$ に応じて：

4.2 変数分離法

基本アイデア: $u(x,t) = X(x)T(t)$ と仮定してPDEを2つのODEに分離します。

4.3 熱方程式

物理的背景: 棒の熱伝導

$\frac{\partial u}{\partial t} = c^2\frac{\partial^2 u}{\partial x^2}, \quad 0 < x < L, \quad t > 0$

境界条件: $u(0,t) = 0$、$u(L,t) = 0$ (両端が温度0に固定)

初期条件: $u(x,0) = f(x)$

解法:

$u = X(x)T(t)$ を代入：

$XT' = c^2 X''T$

$\frac{T'}{c^2 T} = \frac{X''}{X} = -\lambda$

分離定数 $-\lambda$ は境界条件から決まります。

X方程式: $X'' + \lambda X = 0$、$X(0) = X(L) = 0$

固有値問題: $\lambda_n = (n\pi/L)^2$、固有関数: $X_n(x) = \sin(n\pi x/L)$

T方程式: $T' + c^2\lambda_n T = 0$

$T_n(t) = e^{-c^2(n\pi/L)^2 t}$

一般解（重ね合わせの原理）:

$u(x,t) = \sum_{n=1}^{\infty}b_n \sin\frac{n\pi x}{L} e^{-c^2(n\pi/L)^2 t}$

係数決定（$t=0$ を適用）:

$u(x,0) = f(x) = \sum_{n=1}^{\infty}b_n\sin\frac{n\pi x}{L}$

これはフーリエ正弦級数なので：

$b_n = \frac{2}{L}\int_0^L f(x)\sin\frac{n\pi x}{L}\,dx$

電子デバイスの熱解析への応用:

CPUのヒートシンク、パワートランジスタの熱分布、PCB熱管理などに熱方程式が使われます。

物理的解釈:

• 高周波モード（$n$ が大きい項）: 速く減衰（指数項の $n^2$ 依存性）
• 低周波モード（$n = 1$）: 最も遅く減衰
• 時間が経つと初期分布によらず熱が一様に分布

4.4 波動方程式

物理的背景: 張ったひもの横振動

$\frac{\partial^2 u}{\partial t^2} = c^2\frac{\partial^2 u}{\partial x^2}, \quad 0 < x < L$

波速 $c = \sqrt{T/\rho}$ (張力/線密度)

境界条件: $u(0,t) = u(L,t) = 0$ (両端固定)

初期条件: $u(x,0) = f(x)$、$u_t(x,0) = g(x)$

解法（変数分離）:

$u(x,t) = \sum_{n=1}^{\infty}\left(A_n\cos\frac{n\pi ct}{L} + B_n\sin\frac{n\pi ct}{L}\right)\sin\frac{n\pi x}{L}$

係数決定： $A_n = \frac{2}{L}\int_0^L f(x)\sin\frac{n\pi x}{L}\,dx$

$B_n = \frac{2}{n\pi c}\int_0^L g(x)\sin\frac{n\pi x}{L}\,dx$

ダランベールの解（d'Alembert Solution）:

初期変位のみで初期速度がない場合： $u(x,t) = \frac{1}{2}[f(x+ct) + f(x-ct)]$

これは進行波の重ね合わせです。右に $ct$ だけ進んだ波と左に進んだ波の和。

定在波（Standing Wave）:

$u_n(x,t) = A_n\cos\frac{n\pi ct}{L}\sin\frac{n\pi x}{L}$

節（node）: $\sin(n\pi x/L) = 0$ の点 → $x = 0, L/n, 2L/n, \ldots, L$

腹（antinode）: $|\sin(n\pi x/L)| = 1$ の点

電磁波の伝播: 伝送線路（transmission line）の波動方程式：

$\frac{\partial^2 V}{\partial x^2} = LC\frac{\partial^2 V}{\partial t^2}$

ここで $c_{em} = 1/\sqrt{LC}$ (単位長さあたりのインダクタンス/キャパシタンス)

4.5 ラプラス方程式

$\nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0$

物理的意味: 定常状態（時間不変）の熱分布、静電位、非圧縮性流体の速度ポテンシャル

矩形領域（$0 \leq x \leq a$、$0 \leq y \leq b$）:

境界条件： $u(0,y) = 0, \quad u(a,y) = 0, \quad u(x,0) = 0, \quad u(x,b) = f(x)$

解法:

$u = X(x)Y(y)$ に分離：

$\frac{X''}{X} = -\frac{Y''}{Y} = -\lambda$

X方程式: $X'' + \lambda X = 0$、$X(0) = X(a) = 0$

$\lambda_n = \left(\frac{n\pi}{a}\right)^2, \quad X_n = \sin\frac{n\pi x}{a}$

Y方程式: $Y'' - \lambda_n Y = 0$、$Y(0) = 0$

$Y_n = \sinh\frac{n\pi y}{a}$

一般解： $u(x,y) = \sum_{n=1}^{\infty}c_n\sinh\frac{n\pi y}{a}\sin\frac{n\pi x}{a}$

$y = b$ で： $c_n = \frac{2}{a\sinh(n\pi b/a)}\int_0^a f(x)\sin\frac{n\pi x}{a}\,dx$

静電気への応用: 2枚の導体平板間の電位分布、PCB配線間の電界解析

4.6 PythonでPDEを数値的に解く

有限差分法で熱方程式を解く:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

def solve_heat_equation_fd(L, T_total, Nx, Nt, c, initial_func):
    """
    有限差分法で熱方程式を数値的に解く
    u_t = c^2 * u_xx
    境界条件: u(0,t) = u(L,t) = 0
    """
    dx = L / (Nx - 1)
    dt = T_total / (Nt - 1)

    # 安定性条件の確認（CFL条件）
    r = c**2 * dt / dx**2
    if r > 0.5:
        print(f"警告: r = {r:.3f} > 0.5、不安定になる可能性があります！")
    else:
        print(f"r = {r:.3f} (安定条件を満足)")

    x = np.linspace(0, L, Nx)
    u = initial_func(x)
    u[0] = u[-1] = 0  # 境界条件

    u_history = [u.copy()]

    for _ in range(Nt - 1):
        u_new = u.copy()
        u_new[1:-1] = u[1:-1] + r * (u[2:] - 2*u[1:-1] + u[:-2])
        u_new[0] = u_new[-1] = 0
        u = u_new
        u_history.append(u.copy())

    return x, np.array(u_history)

# パラメータ設定
L = 1.0       # 棒の長さ (m)
T_total = 0.1 # シミュレーション時間 (s)
c = 1.0       # 熱拡散係数

def initial_temp(x):
    """初期温度分布（正弦半波）"""
    return np.sin(np.pi * x / L)

Nx, Nt = 50, 500
x, u_hist = solve_heat_equation_fd(L, T_total, Nx, Nt, c, initial_temp)

# 可視化
t_points = np.linspace(0, T_total, Nt)
t_indices = [0, Nt//10, Nt//4, Nt//2, Nt-1]

plt.figure(figsize=(10, 6))
for idx in t_indices:
    plt.plot(x, u_hist[idx],
             label=f't = {t_points[idx]:.3f} s')

plt.xlabel('位置 x (m)')
plt.ylabel('温度 u(x,t)')
plt.title('熱方程式の数値解 - 有限差分法')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('heat_equation.png', dpi=150)
plt.show()

# 解析解との比較
print("\n解析解 vs 数値解の比較 (t=0.05, n=1項):")
t_check = 0.05
idx_check = int(t_check / T_total * (Nt - 1))
u_analytical = np.sin(np.pi * x / L) * np.exp(-c**2 * (np.pi/L)**2 * t_check)
max_error = np.max(np.abs(u_hist[idx_check] - u_analytical))
print(f"最大誤差: {max_error:.6f}")

有限差分法でラプラス方程式（静的電位）を解く:

import numpy as np
import matplotlib.pyplot as plt

def solve_laplace_fd(Nx, Ny, boundary_func, tol=1e-6, max_iter=10000):
    """
    ガウス-ザイデル反復法でラプラス方程式を解く
    """
    u = np.zeros((Ny, Nx))
    x = np.linspace(0, 1, Nx)
    y = np.linspace(0, 1, Ny)

    # 境界条件の設定
    for j in range(Nx):
        u[0, j] = 0              # 下側境界
        u[-1, j] = np.sin(np.pi * x[j])  # 上側境界
    for i in range(Ny):
        u[i, 0] = 0              # 左側境界
        u[i, -1] = 0             # 右側境界

    # 反復法
    for iteration in range(max_iter):
        u_old = u.copy()
        # 内部点の更新
        u[1:-1, 1:-1] = 0.25 * (u[2:, 1:-1] + u[:-2, 1:-1] +
                                  u[1:-1, 2:] + u[1:-1, :-2])
        # 境界条件の再設定
        u[0, :] = 0
        u[-1, :] = np.sin(np.pi * x)
        u[:, 0] = 0
        u[:, -1] = 0

        # 収束確認
        residual = np.max(np.abs(u - u_old))
        if residual < tol:
            print(f"収束: {iteration+1}回の反復後、残差 {residual:.2e}")
            break

    return x, y, u

x, y, u = solve_laplace_fd(50, 50, None)
X, Y = np.meshgrid(x, y)

plt.figure(figsize=(10, 8))
plt.subplot(1, 2, 1)
plt.contourf(X, Y, u, 20, cmap='hot')
plt.colorbar(label='電位 (V)')
plt.xlabel('x')
plt.ylabel('y')
plt.title('ラプラス方程式 - 等電位線')

plt.subplot(1, 2, 2)
plt.contour(X, Y, u, 20)
# 電界ベクトル（電位の負の勾配）
Ey, Ex = np.gradient(-u)
plt.quiver(X[::4, ::4], Y[::4, ::4],
           Ex[::4, ::4], Ey[::4, ::4],
           alpha=0.7)
plt.xlabel('x')
plt.ylabel('y')
plt.title('電界ベクトル')

plt.tight_layout()
plt.savefig('laplace_solution.png', dpi=150)
plt.show()

5. ラプラス変換とPDE、Z変換の紹介

5.1 ラプラス変換でPDEを解く

熱方程式の時間変数についてラプラス変換を適用します：

$\frac{\partial u}{\partial t} = c^2\frac{\partial^2 u}{\partial x^2}$

$sU(x,s) - u(x,0) = c^2\frac{d^2U}{dx^2}$

これは $x$ に関する2階ODEになります。境界条件と一緒に解くと：

$U(x,s) = \frac{f_0}{s}\frac{\sinh(\sqrt{s/c^2}\cdot x)}{\sinh(\sqrt{s/c^2}\cdot L)}$

逆ラプラス変換で時間応答を求めることができます。

5.2 Z変換の紹介（デジタルシステム）

アナログのラプラス変換のデジタル版がZ変換です。

定義: $X(z) = \mathcal{Z}\{x[n]\} = \sum_{n=-\infty}^{\infty}x[n]z^{-n}$

核心的関係: $z = e^{sT_s}$ ($T_s$: サンプリング周期)

| アナログ（ラプラス変換） | デジタル（Z変換） | | ------------------------ | ----------------- | --- | --- | | $s$ 平面 | $z$ 平面 | | $j\omega$ 軸（安定境界） | 単位円 $| z | =1$ | | 左半平面（安定領域） | 単位円の内部 | | $1/s$ (積分器) | $z/(z-1)$ |

主要変換対:

$\mathcal{Z}\{\delta[n]\} = 1$

$\mathcal{Z}\{u[n]\} = \frac{z}{z-1} \quad (|z| > 1)$

$\mathcal{Z}\{a^n u[n]\} = \frac{z}{z-a} \quad (|z| > |a|)$

$\mathcal{Z}\{\sin\Omega_0 n \cdot u[n]\} = \frac{z\sin\Omega_0}{z^2 - 2z\cos\Omega_0 + 1}$

Z変換は次回（3編：複素解析学とZ変換）で完全に扱います。

6. 工学応用総合：信号処理パイプライン

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as sig

# 完全なデジタル信号処理パイプラインの例

# 1. 信号の生成
fs = 8000  # オーディオサンプリング周波数
duration = 0.5
t = np.linspace(0, duration, int(fs * duration), endpoint=False)

# 440 Hz (A4音) + 880 Hz (A5音) + 高周波ノイズ
audio = (np.sin(2 * np.pi * 440 * t) +
         0.5 * np.sin(2 * np.pi * 880 * t) +
         0.3 * np.random.randn(len(t)))

# 2. FFTでスペクトル解析
N = len(audio)
fft_audio = np.fft.rfft(audio)
freqs = np.fft.rfftfreq(N, 1/fs)
magnitude = np.abs(fft_audio)

# 3. バターワースフィルタ設計（1000 Hz以下を通過）
sos = sig.butter(8, 1000, fs=fs, btype='low', output='sos')

# 4. フィルタ適用
audio_filtered = sig.sosfilt(sos, audio)

# 5. 周波数応答確認
w, h = sig.sosfreqz(sos, worN=2000, fs=fs)

# 可視化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

axes[0, 0].plot(t[:400], audio[:400], 'b-', alpha=0.7)
axes[0, 0].set_title('元信号（時間領域）')
axes[0, 0].set_xlabel('時間 (s)')
axes[0, 0].grid(True, alpha=0.3)

axes[0, 1].semilogy(freqs, magnitude)
axes[0, 1].set_title('元信号FFT振幅スペクトル')
axes[0, 1].set_xlabel('周波数 (Hz)')
axes[0, 1].set_xlim([0, 2000])
axes[0, 1].grid(True, alpha=0.3)

axes[1, 0].plot(t[:400], audio_filtered[:400], 'r-', linewidth=1.5)
axes[1, 0].set_title('フィルタ後の信号 (1000Hz LPF)')
axes[1, 0].set_xlabel('時間 (s)')
axes[1, 0].grid(True, alpha=0.3)

axes[1, 1].plot(w, 20 * np.log10(abs(h) + 1e-10), 'g-', linewidth=2)
axes[1, 1].axvline(x=1000, color='r', linestyle='--', label='遮断周波数')
axes[1, 1].set_title('バターワースフィルタの周波数応答')
axes[1, 1].set_xlabel('周波数 (Hz)')
axes[1, 1].set_ylabel('振幅 (dB)')
axes[1, 1].set_xlim([0, 3000])
axes[1, 1].set_ylim([-80, 5])
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('signal_processing_pipeline.png', dpi=150)
plt.show()

まとめと次のステップ

この記事で扱った内容：

1. フーリエ級数: 方形波・三角波の展開、収束、ギブス現象、複素表現、パーセバルの定理
2. フーリエ変換: 定義、主要変換対（矩形波、ガウス関数、デルタ関数）、9つの性質
3. DFTとFFT: 定義、クーリー-テューキーアルゴリズム、$O(N\log N)$ 計算量、Python実装
4. PDEの分類: 楕円型・放物線型・双曲線型
5. 熱方程式: 変数分離法、フーリエ級数による解法、有限差分法
6. 波動方程式: 定在波、ダランベールの解、電磁波
7. ラプラス方程式: 矩形領域の解法、ガウス-ザイデル数値解法
8. Z変換の紹介: アナログ-デジタル対応関係

次回は複素解析学とZ変換を完全に扱い、留数定理で難しい実積分を解き、Z変換でデジタルフィルタを設計する方法を学びます。

参考文献

• Oppenheim, A. & Willsky, A. "Signals and Systems", 2nd Edition, Prentice Hall
• Kreyszig, E. "Advanced Engineering Mathematics", 10th Edition, Wiley
• Proakis, J. & Manolakis, D. "Digital Signal Processing", 4th Edition
• NumPy FFT 公式ドキュメント
• SciPy 信号処理ドキュメント

練習問題