💡 왼쪽 원문을 읽으면서 오른쪽에 따라 써보세요. Tab 키로 힌트를 받을 수 있습니다.

工程数学完全征服第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)$ 满足以下条件，其傅里叶级数就会收敛：

在一个周期内绝对可积： $\int_{-L}^{L}|f(x)|\,dx < \infty$
在一个周期内极值的数量有限
在一个周期内不连续点的数量有限

在不连续点处，级数收敛于左极限与右极限的平均值： $f(x_0) = \frac{f(x_0^-) + f(x_0^+)}{2}$

吉布斯现象：在不连续点附近，对有限项求和时会出现过冲（overshoot）。无论把项数增加到多少，在不连续点附近约 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{其中 } (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（快速傅里叶变换）：复杂度为 $O(N\log_2 N)$

当 $N = 2^{10} = 1024$ 时：

DFT：约 $1024^2 \approx 10^6$ 次运算
FFT：约 $1024 \times 10 \approx 10^4$ 次运算 → 快 100 倍

分治原理（时间抽取法，decimation-in-time）：

按奇偶索引分离：

$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 的分类

二阶线性 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 分离变量法（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$ 的二阶 ODE。结合边界条件求解，可得：

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

通过逆拉普拉斯变换即可求得时间响应。

5.2 Z 变换简介（数字系统）

模拟拉普拉斯变换在数字世界中的对应物，就是 Z 变换。

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

核心关系： $z = e^{sT_s}$ （ $T_s$ ：采样周期）

| 模拟（拉普拉斯） | 数字（Z 变换） | | ------------------------ | --------------- | --- | --- | | $s$ 平面 | $z$ 平面 | | $j\omega$ 轴（稳定边界） | 单位圆 $| z | =1$ | | 左半平面（稳定区域） | 单位圆内部 | | $1/s$ （积分器） | $z/(z-1)$ |

主要变换对：

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

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

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

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

Z 变换将在下一篇（第 3 篇：复分析与 Z 变换）中完整讲解。

6. 工程应用综合：信号处理流水线

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

# 完整的数字信号处理流水线示例

# 1. 生成信号
fs = 8000  # 音频采样频率
duration = 0.5
t = np.linspace(0, duration, int(fs * duration), endpoint=False)

# 440 Hz（A4 音）+ 880 Hz（A5 音）+ 高频噪声
audio = (np.sin(2 * np.pi * 440 * t) +
         0.5 * np.sin(2 * np.pi * 880 * t) +
         0.3 * np.random.randn(len(t)))

# 2. 用 FFT 做频谱分析
N = len(audio)
fft_audio = np.fft.rfft(audio)
freqs = np.fft.rfftfreq(N, 1/fs)
magnitude = np.abs(fft_audio)

# 3. 设计巴特沃斯滤波器（通过 1000 Hz 以下）
sos = sig.butter(8, 1000, fs=fs, btype='low', output='sos')

# 4. 应用滤波器
audio_filtered = sig.sosfilt(sos, audio)

# 5. 检查频率响应
w, h = sig.sosfreqz(sos, worN=2000, fs=fs)

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

axes[0, 0].plot(t[:400], audio[:400], 'b-', alpha=0.7)
axes[0, 0].set_title('原始信号（时域）')
axes[0, 0].set_xlabel('时间 (s)')
axes[0, 0].grid(True, alpha=0.3)

axes[0, 1].semilogy(freqs, magnitude)
axes[0, 1].set_title('原始信号 FFT 幅值谱')
axes[0, 1].set_xlabel('频率 (Hz)')
axes[0, 1].set_xlim([0, 2000])
axes[0, 1].grid(True, alpha=0.3)

axes[1, 0].plot(t[:400], audio_filtered[:400], 'r-', linewidth=1.5)
axes[1, 0].set_title('滤波后信号 (1000Hz LPF)')
axes[1, 0].set_xlabel('时间 (s)')
axes[1, 0].grid(True, alpha=0.3)

axes[1, 1].plot(w, 20 * np.log10(abs(h) + 1e-10), 'g-', linewidth=2)
axes[1, 1].axvline(x=1000, color='r', linestyle='--', label='截止频率')
axes[1, 1].set_title('巴特沃斯滤波器频率响应')
axes[1, 1].set_xlabel('频率 (Hz)')
axes[1, 1].set_ylabel('幅值 (dB)')
axes[1, 1].set_xlim([0, 3000])
axes[1, 1].set_ylim([-80, 5])
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

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

总结与下一步

本文涵盖的内容：

傅里叶级数：方波/三角波展开、收敛性、吉布斯现象、复数表示、帕塞瓦尔定理
傅里叶变换：定义、主要变换对（矩形脉冲、高斯函数、δ 函数）、9 条性质
DFT 与 FFT：定义、库利-图基算法、 $O(N\log N)$ 复杂度、Python 实现
PDE 分类：椭圆型/抛物型/双曲型
热方程：分离变量法、傅里叶级数解法、有限差分法
波动方程：驻波、达朗贝尔解、电磁波
拉普拉斯方程：矩形区域的解法、高斯-赛德尔数值解法
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 信号处理文档

测验

Q1. 吉布斯现象是什么？为什么会发生？

答案：吉布斯现象是指：当傅里叶级数被截断为有限项时，在不连续点附近会出现约 8.9% 的过冲。增加项数并不会消除这个过冲——它只会让过冲的位置更靠近不连续点。

解析：在跳跃不连续点处，傅里叶级数收敛于左右极限的平均值。由于级数在该点收敛较慢，有限项部分和会在不连续点附近出现振荡。这正是方波信号滤波器中出现振铃伪影的原因，因此实践中会使用窗函数（汉明窗、汉宁窗等）来削弱吉布斯效应。

Q2. 在傅里叶变换的语境下，卷积定理是什么？它为什么是信号处理的基础？

答案：傅里叶卷积定理指出：两个函数卷积的傅里叶变换，等于它们各自傅里叶变换的逐点乘积：F{(f*g)(t)} = F(omega) * G(omega)。

解析：在线性时不变（LTI）系统中，时域上的输出信号是输入与系统冲激响应的卷积。而在频域中，这变成了简单的乘法：H(omega) * X(omega) = Y(omega)。这就把计算代价高昂的卷积运算，转化成了简单的乘法运算——这正是实践中采用频域滤波的原因：对输入做 FFT，乘以滤波器响应，再做逆 FFT。

Q3. 分离变量法是如何应用于求解热方程的？

答案：假设 u(x,t) = X(x)T(t)。代入 PDE 并除以 XT，得到 T'/(c^2T) = X''/X = -lambda（一个常数）。由此分离为两个 ODE：X'' + lambdaX = 0（特征值问题）和 T' + c^2lambdaT = 0。代入边界条件，可得特征值 lambda_n = (npi/L)^2 与特征函数 X_n = sin(npix/L)。通解为叠加形式：u(x,t) = sum of b_n * sin(npix/L) _ exp(-c^2_(npi/L)^2 * t)。

解析：傅里叶系数 b_n 由初始条件 u(x,0) = f(x)，通过傅里叶正弦级数公式确定。这个解表明，高频模式衰减得快得多（指数因子中含有 n^2），因此温度分布会随时间趋于平滑。

Q4. 奈奎斯特采样定理是什么？违反它会发生什么？

答案：奈奎斯特定理指出，要如实重建一个信号，采样频率必须至少是信号中最高频率的两倍：f_s 大于等于 2*f_max。这个 f_max 被称为奈奎斯特频率。

解析：如果采样频率过低（f_s 小于 2*f_max），就会发生混叠（aliasing）。高频成分会折返到低频区域，表现为虚假的低频信号，导致原始信号无法恢复。实践中会在采样前使用抗混叠低通滤波器，去除高于 f_s/2 的频率成分。

Q5. 二阶线性 PDE 的三种分类是什么？各自代表哪类物理现象？

答案：二阶线性 PDE 按判别式 B^2 - 4AC 分类：

椭圆型（B^2 - 4AC 小于 0）：拉普拉斯方程与泊松方程——表示静电电位、稳态温度、不可压缩流体等稳态分布。
抛物型（B^2 - 4AC 等于 0）：热方程——表示随时间向稳态演化的扩散过程。
双曲型（B^2 - 4AC 大于 0）：波动方程——表示振动弦、声波、电磁波等传播型波动现象。

解析：这种分类决定了应采用的数值方法，以及边界/初始条件的要求。椭圆型方程需要在所有边界上给出边界条件；抛物型方程需要初始条件加边界条件；双曲型方程则需要两个初始条件（位移与速度）加边界条件。

工程数学完全征服 第2篇：傅里叶级数/变换与偏微分方程(PDE)