Engineering Mathematics 2: Fourier Series/Transform and Partial Differential Equations (PDE)
The mathematical foundation of signal processing, communication systems, and electromagnetics is undoubtedly **Fourier Analysis**. The remarkable fact that any periodic signal can be expressed as a sum of sines and cosines is the core of modern engineering. This article systematically covers Fourier series through FFT and partial differential equations.
1. Fourier Series
1.1 The Core Idea
In 1807, Joseph Fourier made a remarkable claim while solving the heat equation: **any periodic function can be expressed as an infinite series of trigonometric functions**.
**Fourier series** of a function $f(x)$ with period $2L$:
$$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)$$
**Fourier coefficients**:
$$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**: This formula works because of the orthogonality of trigonometric functions.
$$\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{always})$$
1.2 Square Wave Expansion
**Problem**: A square wave with period $2\pi$
$$f(x) = \begin{cases} 1 & 0 < x < \pi \\ -1 & -\pi < x < 0 \end{cases}$$
**Calculation**:
Since $f(x)$ is an odd function, $a_n = 0$ for $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{ odd} \\ 0 & n \text{ even} \end{cases}$$
**Result**:
$$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}$$
Even with just the first few terms, the sum approaches the square wave increasingly closely. This is the power of the Fourier series.
**Leibniz formula**: Substituting $x = \pi/2$:
$$\frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \cdots$$
1.3 Triangle Wave Expansion
$$f(x) = |x|, \quad -\pi \leq x \leq \pi$$
Since it is an 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{ odd} \\ 0 & n \text{ even} \end{cases}$$
**Result**:
$$f(x) = \frac{\pi}{2} - \frac{4}{\pi}\left(\cos x + \frac{\cos 3x}{9} + \frac{\cos 5x}{25} + \cdots\right)$$
Unlike the square wave, the coefficients decrease rapidly as $n^2$, giving faster convergence.
1.4 Convergence and the Gibbs Phenomenon
**Dirichlet Conditions**:
A function $f(x)$ converges in its Fourier series if:
1. Absolutely integrable over one period: $\int_{-L}^{L}|f(x)|\,dx < \infty$
2. Finite number of extrema within one period
3. Finite number of discontinuities within one period
At a point of discontinuity, the series converges to the average of the left and right limits:
$$f(x_0) = \frac{f(x_0^-) + f(x_0^+)}{2}$$
**Gibbs Phenomenon**: Near discontinuities, overshoot occurs when summing a finite number of terms. No matter how many terms are added, approximately **8.9%** overshoot near discontinuities persists.
$$\text{Overshoot} \approx \frac{1}{\pi}\int_0^{\pi}\frac{\sin t}{t}\,dt - \frac{1}{2} \approx 0.0895$$
This is the mathematical cause of the ringing phenomenon in square-wave filters.
1.5 Complex Fourier Series
Using Euler's formula:
$$f(x) = \sum_{n=-\infty}^{\infty}c_n e^{in\pi x/L}$$
Complex coefficients:
$$c_n = \frac{1}{2L}\int_{-L}^{L}f(x)e^{-in\pi x/L}\,dx$$
Relationship with $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's Theorem** (energy conservation):
$$\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)$$
The energy in the time domain equals the energy in the frequency domain.
**Basel problem application**: Applying Parseval's theorem to the triangle wave:
$$\sum_{n=1,3,5,\ldots}\frac{1}{n^4} = \frac{\pi^4}{96}$$
2. Fourier Transform
2.1 From Periodic to Aperiodic Functions
Taking the period $2L$ to $L \to \infty$, the Fourier series develops into the Fourier transform.
**Fourier Transform (time to frequency)**:
$$F(\omega) = \mathcal{F}\{f(t)\} = \int_{-\infty}^{\infty}f(t)e^{-i\omega t}\,dt$$
**Inverse Fourier Transform (frequency to time)**:
$$f(t) = \mathcal{F}^{-1}\{F(\omega)\} = \frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)e^{i\omega t}\,d\omega$$
**Expression using frequency $f$** (commonly used in engineering):
$$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 Key Transform Pairs
**Rectangular pulse**:
$$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)$$
A finite pulse in the time domain becomes a sinc function in the frequency domain.
**Uncertainty principle**: The product of the time width and bandwidth has a lower bound.
$$\Delta t \cdot \Delta\omega \geq \frac{1}{2}$$
A shorter pulse (higher time resolution) has wider bandwidth (lower frequency resolution).
**Gaussian function**: Self-dual transform pair
$$f(t) = e^{-at^2} \Leftrightarrow F(\omega) = \sqrt{\frac{\pi}{a}}e^{-\omega^2/(4a)}$$
The Fourier transform of a Gaussian is again a Gaussian.
**Dirac Delta Function**:
$$\delta(t) = 0 \,(t \neq 0), \quad \int_{-\infty}^{\infty}\delta(t)\,dt = 1$$
$$\mathcal{F}\{\delta(t)\} = 1$$
The frequency spectrum of an impulse is uniform at all frequencies (similar to white noise).
$$\mathcal{F}\{1\} = 2\pi\delta(\omega)$$
A pure DC signal has energy only at $\omega = 0$.
2.3 Key Properties
**Linearity**:
$$\mathcal{F}\{af(t) + bg(t)\} = aF(\omega) + bG(\omega)$$
**Time shift**:
$$\mathcal{F}\{f(t-t_0)\} = e^{-i\omega t_0}F(\omega)$$
Delay changes only the phase, not the magnitude.
**Frequency shift (modulation)**:
$$\mathcal{F}\{e^{i\omega_0 t}f(t)\} = F(\omega - \omega_0)$$
**Time scaling**:
$$\mathcal{F}\{f(at)\} = \frac{1}{|a|}F\!\left(\frac{\omega}{a}\right)$$
Compressing a signal in time widens its bandwidth.
**Differentiation**:
$$\mathcal{F}\{f^{(n)}(t)\} = (i\omega)^n F(\omega)$$
Differentiation becomes multiplication in the frequency domain. ODEs become algebraic equations!
**Integration**:
$$\mathcal{F}\left\{\int_{-\infty}^{t}f(\tau)\,d\tau\right\} = \frac{F(\omega)}{i\omega} + \pi F(0)\delta(\omega)$$
**Convolution theorem** (core of signal processing):
$$\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$$
In linear time-invariant (LTI) systems: output = convolution of input with impulse response. In the frequency domain, this is just multiplication!
**Parseval's theorem (energy conservation)**:
$$\int_{-\infty}^{\infty}|f(t)|^2\,dt = \frac{1}{2\pi}\int_{-\infty}^{\infty}|F(\omega)|^2\,d\omega$$
2.4 Application: Filter Design
**Ideal Low-Pass Filter (LPF)**:
$$H(\omega) = \begin{cases} 1 & |\omega| \leq \omega_c \\ 0 & |\omega| > \omega_c \end{cases}$$
Impulse response:
$$h(t) = \mathcal{F}^{-1}\{H(\omega)\} = \frac{\omega_c}{\pi}\text{sinc}\!\left(\frac{\omega_c t}{\pi}\right)$$
Non-causal (requires future samples), so cannot be implemented in practice. Practical filters (Butterworth, Chebyshev, etc.) are used.
**Bandpass filter**:
$$H(\omega) = H_{LPF}(\omega - \omega_0) + H_{LPF}(\omega + \omega_0)$$
Impulse response: $h(t) = 2h_{LPF}(t)\cos\omega_0 t$
**Frequency Response**:
When $e^{i\omega t}$ is input to an LTI system:
$$y(t) = H(\omega)e^{i\omega t}$$
$H(\omega)$ simultaneously represents complex magnitude (gain) and phase shift:
$$H(\omega) = |H(\omega)|e^{i\angle H(\omega)}$$
3. Discrete Fourier Transform (DFT) and FFT
3.1 DFT Definition
For $N$ discrete samples $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}$$
**Frequency resolution**: $\Delta f = f_s / N$ (sampling frequency divided by N)
**Nyquist limit**: $f_s \geq 2 f_{max}$ (to prevent aliasing)
3.2 FFT Algorithm (Cooley-Tukey Algorithm)
Direct DFT computation: $O(N^2)$ complexity
FFT (Fast Fourier Transform): $O(N\log_2 N)$ complexity
For $N = 2^{10} = 1024$:
- DFT: $1024^2 \approx 10^6$ operations
- FFT: $1024 \times 10 \approx 10^4$ operations — **100x faster**
**Divide and conquer principle** (decimation-in-time):
Separate into even/odd indices:
$$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]$$
Since $W_N^2 = W_{N/2}$, $E[k]$ and $O[k]$ are each $N/2$-point DFTs.
Applying this recursively gives $O(N\log N)$.
3.3 FFT Implementation and Signal Analysis with Python
from scipy import signal as sig
Generate composite signal (50 Hz + 120 Hz mixed)
fs = 1000 # Sampling frequency (Hz)
T = 1.0 # Signal length (seconds)
N = int(fs * T)
t = np.linspace(0, T, N, endpoint=False)
Signal: 50 Hz (amplitude 1.0) + 120 Hz (amplitude 0.5) + noise
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
Perform FFT
fft_result = np.fft.fft(signal_noisy)
freqs = np.fft.fftfreq(N, 1/fs)
Two-sided spectrum -> one-sided spectrum
magnitude = (2.0 / N) * np.abs(fft_result[:N//2])
freqs_pos = freqs[:N//2]
Visualization
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('Time (s)')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Time Domain Signal (blue: with noise, red: clean signal)')
axes[0].legend(['Noisy', 'Clean'])
axes[0].grid(True, alpha=0.3)
axes[1].plot(freqs_pos, magnitude, 'g-')
axes[1].set_xlabel('Frequency (Hz)')
axes[1].set_ylabel('Amplitude Spectrum')
axes[1].set_title('Frequency Domain (FFT Magnitude)')
axes[1].set_xlim([0, 200])
axes[1].grid(True, alpha=0.3)
Spectrogram
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('Frequency (Hz)')
axes[2].set_xlabel('Time (s)')
axes[2].set_title('Spectrogram')
axes[2].set_ylim([0, 200])
plt.tight_layout()
plt.savefig('fft_analysis.png', dpi=150)
plt.show()
Print dominant frequency components
peak_indices = np.where(magnitude > 0.1)[0]
print("Dominant frequency components:")
for idx in peak_indices:
print(f" {freqs_pos[idx]:.1f} Hz: amplitude {magnitude[idx]:.3f}")
3.4 Signal Filtering with FFT
def fft_lowpass_filter(signal_in, fs, cutoff_freq):
"""FFT-based ideal low-pass filter"""
N = len(signal_in)
fft_result = np.fft.fft(signal_in)
freqs = np.fft.fftfreq(N, 1/fs)
Zero out components above cutoff frequency
fft_filtered = fft_result.copy()
fft_filtered[np.abs(freqs) > cutoff_freq] = 0
Inverse transform
return np.real(np.fft.ifft(fft_filtered))
Filtering example
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='Original signal')
plt.plot(t[:200], filtered[:200], 'r-', linewidth=2, label='Filtered signal (80Hz LPF)')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.title('FFT-Based Low-Pass Filtering')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
4. Partial Differential Equations (PDE)
4.1 PDE Classification
Second-order linear PDE:
$$Au_{xx} + Bu_{xy} + Cu_{yy} + Du_x + Eu_y + Fu = G$$
Based on discriminant $B^2 - 4AC$:
| Condition | Classification | Representative Equation |
| --------------- | -------------- | ----------------------- |
| $B^2 - 4AC < 0$ | Elliptic | Laplace equation |
| $B^2 - 4AC = 0$ | Parabolic | Heat equation |
| $B^2 - 4AC > 0$ | Hyperbolic | Wave equation |
4.2 Separation of Variables
**Basic idea**: Assume $u(x,t) = X(x)T(t)$ to separate the PDE into two ODEs.
4.3 Heat Equation
**Physical background**: Heat conduction in a rod
$$\frac{\partial u}{\partial t} = c^2\frac{\partial^2 u}{\partial x^2}, \quad 0 < x < L, \quad t > 0$$
**Boundary conditions**: $u(0,t) = 0$, $u(L,t) = 0$ (both ends fixed at temperature 0)
**Initial condition**: $u(x,0) = f(x)$
**Solution procedure**:
Substituting $u = X(x)T(t)$:
$$XT' = c^2 X''T$$
$$\frac{T'}{c^2 T} = \frac{X''}{X} = -\lambda$$
The separation constant $-\lambda$ is determined by the boundary conditions.
**X equation**: $X'' + \lambda X = 0$, $X(0) = X(L) = 0$
Eigenvalue problem: $\lambda_n = (n\pi/L)^2$, eigenfunctions: $X_n(x) = \sin(n\pi x/L)$
**T equation**: $T' + c^2\lambda_n T = 0$
$$T_n(t) = e^{-c^2(n\pi/L)^2 t}$$
**General solution** (superposition principle):
$$u(x,t) = \sum_{n=1}^{\infty}b_n \sin\frac{n\pi x}{L} e^{-c^2(n\pi/L)^2 t}$$
**Coefficient determination** (applying $t=0$):
$$u(x,0) = f(x) = \sum_{n=1}^{\infty}b_n\sin\frac{n\pi x}{L}$$
This is a Fourier sine series, so:
$$b_n = \frac{2}{L}\int_0^L f(x)\sin\frac{n\pi x}{L}\,dx$$
**Application to electronic component thermal analysis**:
The heat equation is used in CPU heatsink design, power transistor heat distribution, and PCB thermal management.
**Physical interpretation**:
- High-frequency modes (large $n$): decay rapidly (the $n^2$ dependence in the exponential term)
- Low-frequency mode ($n = 1$): decays most slowly
- Over time, heat distributes uniformly regardless of the initial distribution
4.4 Wave Equation
**Physical background**: Transverse vibration of a taut string
$$\frac{\partial^2 u}{\partial t^2} = c^2\frac{\partial^2 u}{\partial x^2}, \quad 0 < x < L$$
Wave speed $c = \sqrt{T/\rho}$ (tension/linear mass density)
**Boundary conditions**: $u(0,t) = u(L,t) = 0$ (both ends fixed)
**Initial conditions**: $u(x,0) = f(x)$, $u_t(x,0) = g(x)$
**Solution** (separation of variables):
$$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}$$
Coefficient determination:
$$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**:
For the case of initial displacement only and zero initial velocity:
$$u(x,t) = \frac{1}{2}[f(x+ct) + f(x-ct)]$$
This is the superposition of **traveling waves** — the sum of a wave moving to the right by $ct$ and a wave moving to the left.
**Standing Wave**:
$$u_n(x,t) = A_n\cos\frac{n\pi ct}{L}\sin\frac{n\pi x}{L}$$
Nodes: points where $\sin(n\pi x/L) = 0$ — at $x = 0, L/n, 2L/n, \ldots, L$
Antinodes: points where $|\sin(n\pi x/L)| = 1$
**Electromagnetic wave propagation**: Wave equation on a transmission line:
$$\frac{\partial^2 V}{\partial x^2} = LC\frac{\partial^2 V}{\partial t^2}$$
Where $c_{em} = 1/\sqrt{LC}$ (per-unit-length inductance/capacitance).
4.5 Laplace Equation
$$\nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0$$
**Physical meaning**: Steady-state (time-invariant) heat distribution, electrostatic potential, velocity potential in incompressible flow.
**Rectangular domain** ($0 \leq x \leq a$, $0 \leq y \leq b$):
Boundary conditions:
$$u(0,y) = 0, \quad u(a,y) = 0, \quad u(x,0) = 0, \quad u(x,b) = f(x)$$
**Solution**:
Separating $u = X(x)Y(y)$:
$$\frac{X''}{X} = -\frac{Y''}{Y} = -\lambda$$
X equation: $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 equation: $Y'' - \lambda_n Y = 0$, $Y(0) = 0$
$$Y_n = \sinh\frac{n\pi y}{a}$$
General solution:
$$u(x,y) = \sum_{n=1}^{\infty}c_n\sinh\frac{n\pi y}{a}\sin\frac{n\pi x}{a}$$
At $y = b$:
$$c_n = \frac{2}{a\sinh(n\pi b/a)}\int_0^a f(x)\sin\frac{n\pi x}{a}\,dx$$
**Electrostatic application**: Potential distribution between two conducting plates, electric field analysis between PCB traces.
4.6 Numerical PDE Solutions with Python
**Solving the heat equation with the Finite Difference Method**:
from matplotlib.animation import FuncAnimation
def solve_heat_equation_fd(L, T_total, Nx, Nt, c, initial_func):
"""
Numerical solution of heat equation using finite difference method
u_t = c^2 * u_xx
Boundary conditions: u(0,t) = u(L,t) = 0
"""
dx = L / (Nx - 1)
dt = T_total / (Nt - 1)
Check stability condition (CFL condition)
r = c**2 * dt / dx**2
if r > 0.5:
print(f"Warning: r = {r:.3f} > 0.5, may be unstable!")
else:
print(f"r = {r:.3f} (stability condition satisfied)")
x = np.linspace(0, L, Nx)
u = initial_func(x)
u[0] = u[-1] = 0 # Boundary conditions
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)
Parameter setup
L = 1.0 # Rod length (m)
T_total = 0.1 # Simulation time (s)
c = 1.0 # Thermal diffusivity
def initial_temp(x):
"""Initial temperature distribution (half sine wave)"""
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)
Visualization
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('Position x (m)')
plt.ylabel('Temperature u(x,t)')
plt.title('Heat Equation Numerical Solution - Finite Difference Method')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('heat_equation.png', dpi=150)
plt.show()
Compare with analytical solution
print("\nAnalytical vs numerical comparison (t=0.05, n=1 term):")
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"Maximum error: {max_error:.6f}")
**Solving the Laplace equation (static potential) with finite difference**:
def solve_laplace_fd(Nx, Ny, boundary_func, tol=1e-6, max_iter=10000):
"""
Solving Laplace equation using Gauss-Seidel iteration
"""
u = np.zeros((Ny, Nx))
x = np.linspace(0, 1, Nx)
y = np.linspace(0, 1, Ny)
Set boundary conditions
for j in range(Nx):
u[0, j] = 0 # Bottom boundary
u[-1, j] = np.sin(np.pi * x[j]) # Top boundary
for i in range(Ny):
u[i, 0] = 0 # Left boundary
u[i, -1] = 0 # Right boundary
Iteration
for iteration in range(max_iter):
u_old = u.copy()
Update interior points
u[1:-1, 1:-1] = 0.25 * (u[2:, 1:-1] + u[:-2, 1:-1] +
u[1:-1, 2:] + u[1:-1, :-2])
Re-apply boundary conditions
u[0, :] = 0
u[-1, :] = np.sin(np.pi * x)
u[:, 0] = 0
u[:, -1] = 0
Check convergence
residual = np.max(np.abs(u - u_old))
if residual < tol:
print(f"Converged: residual {residual:.2e} after {iteration+1} iterations")
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='Potential (V)')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Laplace Equation - Equipotential Lines')
plt.subplot(1, 2, 2)
plt.contour(X, Y, u, 20)
Electric field vector (negative gradient of potential)
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('Electric Field Vectors')
plt.tight_layout()
plt.savefig('laplace_solution.png', dpi=150)
plt.show()
5. Laplace Transform and PDE, Introduction to Z-Transform
5.1 Solving PDE with the Laplace Transform
Applying the Laplace transform to the heat equation with respect to the time variable:
$$\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}$$
This becomes a second-order ODE in $x$. Solving with boundary conditions:
$$U(x,s) = \frac{f_0}{s}\frac{\sinh(\sqrt{s/c^2}\cdot x)}{\sinh(\sqrt{s/c^2}\cdot L)}$$
The time response can be obtained by inverse Laplace transform.
5.2 Introduction to Z-Transform (Digital Systems)
The digital counterpart of the analog Laplace transform is the **Z-transform**.
**Definition**:
$$X(z) = \mathcal{Z}\{x[n]\} = \sum_{n=-\infty}^{\infty}x[n]z^{-n}$$
**Key relationship**: $z = e^{sT_s}$ ($T_s$: sampling period)
| Analog (Laplace) | Digital (Z-transform) |
| ----------------------------------- | --------------------- | --- | --- |
| $s$ plane | $z$ plane |
| $j\omega$ axis (stability boundary) | Unit circle $ | z | =1$ |
| Left half-plane (stable region) | Inside unit circle |
| $1/s$ (integrator) | $z/(z-1)$ |
**Key transform pairs**:
$$\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}$$
The Z-transform will be fully covered in the next article (Part 3: Complex Analysis and Z-Transform).
6. Engineering Application: Signal Processing Pipeline
from scipy import signal as sig
Complete digital signal processing pipeline example
1. Signal generation
fs = 8000 # Audio sampling frequency
duration = 0.5
t = np.linspace(0, duration, int(fs * duration), endpoint=False)
440 Hz (A4 note) + 880 Hz (A5 note) + high-frequency noise
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. Spectrum analysis with FFT
N = len(audio)
fft_audio = np.fft.rfft(audio)
freqs = np.fft.rfftfreq(N, 1/fs)
magnitude = np.abs(fft_audio)
3. Butterworth filter design (pass below 1000 Hz)
sos = sig.butter(8, 1000, fs=fs, btype='low', output='sos')
4. Apply filter
audio_filtered = sig.sosfilt(sos, audio)
5. Check frequency response
w, h = sig.sosfreqz(sos, worN=2000, fs=fs)
Visualization
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('Original Signal (Time Domain)')
axes[0, 0].set_xlabel('Time (s)')
axes[0, 0].grid(True, alpha=0.3)
axes[0, 1].semilogy(freqs, magnitude)
axes[0, 1].set_title('Original Signal FFT Magnitude Spectrum')
axes[0, 1].set_xlabel('Frequency (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('Filtered Signal (1000Hz LPF)')
axes[1, 0].set_xlabel('Time (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='Cutoff frequency')
axes[1, 1].set_title('Butterworth Filter Frequency Response')
axes[1, 1].set_xlabel('Frequency (Hz)')
axes[1, 1].set_ylabel('Magnitude (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()
Summary and Next Steps
Topics covered in this article:
1. **Fourier series**: Square/triangle wave expansion, convergence, Gibbs phenomenon, complex representation, Parseval's theorem
2. **Fourier transform**: Definition, key transform pairs (rectangular pulse, Gaussian, delta function), 9 properties
3. **DFT and FFT**: Definition, Cooley-Tukey algorithm, $O(N\log N)$ complexity, Python implementation
4. **PDE classification**: Elliptic/Parabolic/Hyperbolic
5. **Heat equation**: Separation of variables, Fourier series solution, finite difference method
6. **Wave equation**: Standing waves, d'Alembert solution, electromagnetic waves
7. **Laplace equation**: Rectangular domain solution, Gauss-Seidel numerical solution
8. **Introduction to Z-transform**: Analog-digital correspondence
The next article will fully cover **Complex Analysis and Z-Transform** — using the residue theorem to solve difficult real integrals and designing digital filters with the Z-transform.
References
- 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 Documentation](https://numpy.org/doc/stable/reference/routines.fft.html)
- [SciPy Signal Processing Documentation](https://docs.scipy.org/doc/scipy/reference/signal.html)
Quiz
**Answer**: The Gibbs phenomenon is an overshoot of approximately 8.9% that occurs near discontinuities when a Fourier series is truncated to a finite number of terms. Increasing the number of terms does not eliminate this overshoot — it merely moves closer to the discontinuity.
**Explanation**: At a jump discontinuity, the Fourier series converges to the average of the left and right limits. The finite partial sum exhibits oscillations near the discontinuity due to the slow convergence of the series at that point. This causes ringing artifacts in square-wave signal filters, which is why windowing functions (Hamming, Hanning, etc.) are used in practice to reduce the Gibbs effect.
**Answer**: The Fourier convolution theorem states that the Fourier transform of the convolution of two functions equals the pointwise product of their individual Fourier transforms: F\{(f*g)(t)\} = F(omega) * G(omega).
**Explanation**: In linear time-invariant (LTI) systems, the output signal is the convolution of the input with the system's impulse response in the time domain. In the frequency domain, this becomes simple multiplication H(omega) \* X(omega) = Y(omega). This transforms a computationally expensive convolution operation into a simple multiplication, which is why frequency-domain filtering is used in practice: compute FFT of input, multiply by filter response, compute inverse FFT.
**Answer**: Assume u(x,t) = X(x)*T(t). Substituting into the PDE and dividing by X*T gives T'/(c^2*T) = X''/X = -lambda (a constant). This separates into two ODEs: X'' + lambda*X = 0 (eigenvalue problem) and T' + c^2*lambda*T = 0. Applying boundary conditions gives eigenvalues lambda_n = (n*pi/L)^2 and eigenfunctions X_n = sin(n*pi*x/L). The general solution is the superposition u(x,t) = sum of b_n * sin(n*pi*x/L) _ exp(-c^2_(n*pi/L)^2 * t).
**Explanation**: The Fourier coefficients b_n are determined from the initial condition u(x,0) = f(x) using the Fourier sine series formula. The solution shows that higher-frequency modes decay much faster (exponential factor has n^2), so the temperature distribution smooths out over time.
**Answer**: The Nyquist theorem states that to faithfully reconstruct a signal, the sampling frequency must be at least twice the maximum frequency in the signal: f_s greater than or equal to 2\*f_max. The frequency f_max is called the Nyquist frequency.
**Explanation**: If the sampling frequency is too low (f_s less than 2\*f_max), aliasing occurs. High-frequency components fold back into the lower frequency range and appear as false low-frequency signals, making it impossible to recover the original signal. In practice, anti-aliasing low-pass filters are applied before sampling to remove frequency components above f_s/2.
**Answer**: Second-order linear PDEs are classified by the discriminant B^2 - 4AC:
- Elliptic (B^2 - 4AC less than 0): Laplace and Poisson equations — represent steady-state distributions such as static electric potential, steady-state temperature, and incompressible fluid flow.
- Parabolic (B^2 - 4AC equals 0): Heat equation — represents diffusion processes that evolve over time toward a steady state.
- Hyperbolic (B^2 - 4AC greater than 0): Wave equation — represents propagating wave phenomena such as vibrating strings, sound waves, and electromagnetic waves.
**Explanation**: The classification determines the appropriate numerical methods and boundary/initial condition requirements. Elliptic equations need boundary conditions on all sides; parabolic equations need initial conditions plus boundary conditions; hyperbolic equations need two initial conditions (displacement and velocity) plus boundary conditions.
현재 단락 (1/456)
The mathematical foundation of signal processing, communication systems, and electromagnetics is und...