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

import numpy as np
import matplotlib.pyplot as plt
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

import numpy as np
import matplotlib.pyplot as plt

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

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:

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

import numpy as np
import matplotlib.pyplot as plt

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

import numpy as np
import matplotlib.pyplot as plt
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
• SciPy Signal Processing Documentation

Quiz