공업수학 완전 정복 2편: 푸리에 급수/변환과 편미분방정식(PDE)

공업수학 완전 정복 2편: 푸리에 급수/변환과 편미분방정식(PDE)

신호 처리, 통신 시스템, 전자기학의 수학적 기반은 단연 **푸리에 해석(Fourier Analysis)**입니다. 어떤 주기 신호도 사인·코사인의 합으로 표현할 수 있다는 놀라운 사실이 현대 공학의 핵심입니다. 이번 글에서는 푸리에 급수부터 FFT, 그리고 편미분방정식까지 체계적으로 정리합니다.

1. 푸리에 급수 (Fourier Series)

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 사각파(Square Wave) 전개

문제: 주기 $2\pi$인 사각파

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

계산:

$f(x)$는 기함수(odd function)이므로 $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$를 대입하면: $1 = \frac{4}{\pi}\left(1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \cdots\right)$ $\frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \cdots$

1.3 삼각파(Triangle Wave) 전개

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

우함수(even function)이므로 $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 수렴과 깁스 현상 (Gibbs Phenomenon)

디리클레 조건(Dirichlet Conditions):

함수 $f(x)$가 다음 조건을 만족하면 푸리에 급수는 수렴합니다:

1. 한 주기 내에서 절대 적분 가능: $\int_{-L}^{L}|f(x)|\,dx < \infty$
2. 한 주기 내에서 극값의 수가 유한
3. 한 주기 내에서 불연속점의 수가 유한

불연속점에서 급수는 좌극한과 우극한의 평균으로 수렴합니다: $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)$

파시발(Parseval) 정리 (에너지 보존):

$\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. 푸리에 변환 (Fourier Transform)

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}$

펄스가 짧을수록(시간 분해능 높음) 대역폭이 넓어집니다(주파수 분해능 낮음).

가우시안 함수: 자기 쌍대(self-dual) 변환쌍

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

가우시안의 푸리에 변환은 다시 가우시안입니다.

디락 델타 함수(Dirac Delta):

$\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{where } (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$

주파수 응답(Frequency Response):

LTI 시스템에 $e^{i\omega t}$를 입력하면: $y(t) = H(\omega)e^{i\omega t}$

$H(\omega)$는 복소 크기(gain)와 위상 이동을 동시에 나타냅니다: $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 (Fast Fourier Transform): $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 변수분리법 (Separation of Variables)

기본 아이디어: $u(x,t) = X(x)T(t)$로 가정하여 PDE를 두 개의 ODE로 분리합니다.

4.3 열방정식 (Heat Equation)

물리적 배경: 막대의 열전도

$\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 파동방정식 (Wave Equation)

물리적 배경: 팽팽한 줄의 횡진동

$\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 라플라스 방정식 (Laplace Equation)

$\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$

정전기 응용: 두 도체 평판 사이의 전위 분포, PCB 도선 사이의 전기장 해석

4.6 Python으로 PDE 수치 풀기

유한 차분법(Finite Difference Method)으로 열방정식 풀기:

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):
    """
    가우스-자이델 반복법으로 라플라스 방정식 풀이
    경계 조건: boundary_func(x, y, side) 반환
    """
    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}$

다음 편(3편: 복소해석학과 Z-변환)에서 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 Signal Processing 문서