工業数学完全攻略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を実装して信号を解析する

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を使った信号フィルタリング

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$ に応じて：

| 条件 | 分類 | 代表的な方程式 |

| --------------- | ---------------------- | -------------- |

| $B^2 - 4AC < 0$ | 楕円型（Elliptic） | ラプラス方程式 |

| $B^2 - 4AC = 0$ | 放物線型（Parabolic） | 熱方程式 |

| $B^2 - 4AC > 0$ | 双曲線型（Hyperbolic） | 波動方程式 |

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を数値的に解く

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

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}")

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

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. 工学応用総合：信号処理パイプライン

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 公式ドキュメント](https://numpy.org/doc/stable/reference/routines.fft.html)

- [SciPy 信号処理ドキュメント](https://docs.scipy.org/doc/scipy/reference/signal.html)

練習問題

**答え**: ギブス現象は、フーリエ級数を有限個の項で打ち切ったとき、不連続点近くで約8.9%のオーバーシュートが発生する現象です。項数をいくら増やしてもこのオーバーシュートは消えません（不連続点に近づくだけです）。

**解説**: 不連続点ではフーリエ級数は左極限と右極限の平均に収束します。有限個の部分和は不連続点近くで振動が起きます。これが方形波フィルタでリンギング（ringing）が発生する数学的原因です。実際にはハミング窓やハニング窓などの窓関数を使ってギブス効果を軽減します。

**答え**: 畳み込み定理は、2つの関数の畳み込みのフーリエ変換が、それぞれの関数のフーリエ変換の点ごとの積に等しいことを述べます: F\{(f*g)(t)\} = F(omega) * G(omega)

**解説**: 線形時不変（LTI）システムでは、時間領域の出力信号は入力とシステムのインパルス応答の畳み込みです。周波数領域ではこれが単純な乗算 H(omega) \* X(omega) = Y(omega) になります。これにより周波数領域フィルタリングが実用的になりました：入力のFFTを計算し、フィルタ応答を乗算し、逆FFTを計算するだけです。

**答え**: u(x,t) = X(x)*T(t) と仮定します。PDEに代入してX*Tで割ると T'/(c^2*T) = X''/X = -lambda （定数）となります。これにより2つのODEに分離されます：X'' + lambda*X = 0（固有値問題）とT' + c^2*lambda*T = 0。境界条件を適用すると固有値 lambda_n = (n*pi/L)^2 と固有関数 X_n = sin(n*pi*x/L) が得られます。一般解は重ね合わせ: u(x,t) = sum of b_n * sin(n*pi*x/L) _ exp(-c^2_(n*pi/L)^2 * t)

**解説**: フーリエ係数 b_n は初期条件 u(x,0) = f(x) からフーリエ正弦級数の公式で決まります。高周波モードが速く減衰するため（指数項にn^2がある）、時間が経つと温度分布は滑らかになります。

**答え**: ナイキスト定理は、信号を忠実に再現するには、サンプリング周波数が信号の最大周波数の少なくとも2倍以上必要であることを述べます: f_s は 2\*f_max 以上。この f_max をナイキスト周波数といいます。

**解説**: サンプリング周波数が低すぎると（f_s が 2\*f_max より低い）エイリアシングが発生します。高周波成分が低周波帯域に折り返し、偽の低周波信号として現れるため、元の信号を復元することが不可能になります。実際にはサンプリング前にアンチエイリアシング低域通過フィルタを適用してf_s/2以上の周波数成分を除去します。

**答え**: 2階線形PDEは判別式 B^2 - 4AC で分類されます：

- 楕円型（B^2 - 4AC が0より小さい）：ラプラス方程式・ポアソン方程式 — 定常状態の分布（静電位、定常温度、非圧縮性流体など）を表します。

- 放物線型（B^2 - 4AC が0に等しい）：熱方程式 — 定常状態に向かって時間的に発展する拡散現象を表します。

- 双曲線型（B^2 - 4AC が0より大きい）：波動方程式 — 振動する弦、音波、電磁波などの伝播波動現象を表します。

**解説**: 分類によって適切な数値解法と境界条件・初期条件の要件が決まります。楕円型は全境界に境界条件が必要、放物線型は初期条件と境界条件が必要、双曲線型は2つの初期条件（変位と速度）と境界条件が必要です。