常微分方程（十三）偏微分方程引论
Chen Kai BOSS
  2019-06-09 14:15:00 2019-06-09 14:15      6.6k 字  32 分钟  

当物理量不仅随时间变化，还随空间变化时，我们进入了偏微分方程（ PDE）的领域。热量在金属棒中传导、波在水面上扩散、电磁场在空间中传播——这些现象都需要 PDE 来描述。本章将介绍 PDE 的基本概念、分类，以及三大经典方程的求解方法。

从 ODE 到 PDE

为什么需要偏微分方程？

热传导方程的物理意义

在前面的章节中，我们研究的未知函数只依赖于一个自变量，如 或。但现实世界中，很多物理量同时依赖于多个变量：

  • 温度场：依赖于空间位置和时间
  • 波函数：依赖于位置和时间
  • 流体速度场：三维空间加时间

当未知函数依赖于多个自变量时，微分方程中出现偏导数，这就是偏微分方程（ PDE）。

基本概念

定义：偏微分方程是含有未知函数及其偏导数的方程。

阶：方程中出现的最高阶偏导数的阶数。

线性 vs 非线性： - 线性：未知函数及其导数以线性方式出现 - 非线性：存在未知函数或导数的乘积、幂次等

例子：

  1. 一维热方程（线性，二阶）：

  2. 一维波动方程（线性，二阶）：

  3. Laplace 方程（线性，二阶）：

  4. Burgers 方程（非线性，二阶）：

PDE 的分类

二阶线性 PDE 的标准形式

波动方程的解

二阶线性 PDE 的一般形式（两个自变量）：

分类标准

类比于二次曲线， PDE 根据判别式 分类：

判别式 类型 物理原型 特征
双曲型 波动方程 信息传播，有限速度
抛物型 热方程 扩散过程，不可逆
椭圆型 Laplace 方程 稳态，平衡

物理直觉

双曲型（波动）：信息以有限速度传播。如果你在湖中投一颗石子，涟漪会逐渐扩散，但远处暂时不受影响。

抛物型（扩散）：信息瞬间传播到无穷远（虽然影响很小）。热量会立即开始向周围扩散。

椭圆型（平衡）：没有时间演化，描述的是稳态。如静电场、稳态温度分布。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# 可视化三种类型 PDE 的典型解

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# 1. 双曲型：行波
x = np.linspace(0, 10, 200)
t_values = [0, 1, 2, 3]
for t in t_values:
    u = np.exp(-(x - 2*t - 2)**2)  # 右行波
    axes[0].plot(x, u, label=f't = {t}')
axes[0].set_xlabel('x')
axes[0].set_ylabel('u')
axes[0].set_title('双曲型：行波传播')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 2. 抛物型：热扩散
for t in [0.01, 0.1, 0.5, 1.0]:
    u = np.exp(-x**2 / (4*t)) / np.sqrt(4*np.pi*t)
    axes[1].plot(x - 5, u, label=f't = {t}')
axes[1].set_xlabel('x')
axes[1].set_ylabel('u')
axes[1].set_title('抛物型：热扩散')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim([-5, 5])

# 3. 椭圆型：稳态温度分布
x = np.linspace(0, 1, 50)
y = np.linspace(0, 1, 50)
X, Y = np.meshgrid(x, y)
# Laplace 方程的一个简单解
U = np.sin(np.pi * X) * np.sinh(np.pi * Y) / np.sinh(np.pi)
c = axes[2].contourf(X, Y, U, levels=20, cmap='hot')
plt.colorbar(c, ax=axes[2])
axes[2].set_xlabel('x')
axes[2].set_ylabel('y')
axes[2].set_title('椭圆型：稳态温度分布')

plt.tight_layout()
plt.savefig('pde_types.png', dpi=150, bbox_inches='tight')
plt.show()

初边值问题

为什么需要边界条件？

拉普拉斯方程与调和函数

ODE 需要初始条件来确定解。 PDE 不仅需要初始条件（对于演化型方程），还需要边界条件（因为有空间域）。

边界条件的类型

设区域边界为：

第一类（ Dirichlet）：给定函数值

第二类（ Neumann）：给定法向导数

第三类（ Robin/混合）：函数值和导数的线性组合

适定性

一个 PDE 问题是适定的（ well-posed），如果： 1. 存在性：解存在 2. 唯一性：解唯一 3. 稳定性：解连续依赖于数据

热方程

物理推导

分离变量法示意

考虑一维细杆的热传导。设 为 点在 时刻的温度。

Fourier 热传导定律：热流量与温度梯度成正比

能量守恒：

组合得到热方程：

其中 是热扩散率。

分离变量法

考虑问题：$

$

假设解的形式： 代入方程：

分离变量： （两边必须等于同一常数，因为左边只依赖，右边只依赖）

得到两个 ODE：$

$

求解空间部分（ Sturm-Liouville 问题）：

边界条件 导致：

求解时间部分：

一般解：

确定系数：由初始条件Fourier 正弦级数系数： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def heat_equation_analytical(f, L, alpha, x, t, n_terms=50):
    """
    热方程的解析解（分离变量法）
    u_t = alpha * u_xx
    u(0,t) = u(L,t) = 0
    u(x,0) = f(x)
    """
    # 计算 Fourier 系数
    from scipy.integrate import quad
    
    def integrand(xi, n):
        return f(xi) * np.sin(n * np.pi * xi / L)
    
    u = np.zeros_like(x, dtype=float)
    
    for n in range(1, n_terms + 1):
        Bn, _ = quad(integrand, 0, L, args=(n,))
        Bn *= 2 / L
        
        lambda_n = (n * np.pi / L)**2
        Xn = np.sin(n * np.pi * x / L)
        Tn = np.exp(-alpha * lambda_n * t)
        
        u += Bn * Xn * Tn
    
    return u

# 示例：初始温度分布为三角形
L = 1
alpha = 0.01

def f_initial(x):
    """三角形初始条件"""
    if x <= 0.5:
        return 2 * x
    else:
        return 2 * (1 - x)

f_vec = np.vectorize(f_initial)

# 绘制不同时刻的温度分布
x = np.linspace(0, L, 100)
times = [0, 0.1, 0.5, 1.0, 2.0, 5.0]

plt.figure(figsize=(10, 6))
for t in times:
    if t == 0:
        u = f_vec(x)
    else:
        u = heat_equation_analytical(f_initial, L, alpha, x, t)
    plt.plot(x, u, label=f't = {t}')

plt.xlabel('位置 x')
plt.ylabel('温度 u')
plt.title('热方程的解：三角形初始条件')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('heat_analytical.png', dpi=150, bbox_inches='tight')
plt.show()

有限差分法求解热方程

显式格式（ FTCS）：

用向前差分近似时间导数，中心差分近似空间导数：

整理得：

其中。

稳定性条件： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def heat_ftcs(alpha, L, T, Nx, Nt, f, u_left=0, u_right=0):
    """
    热方程的 FTCS 显式格式
    """
    dx = L / Nx
    dt = T / Nt
    r = alpha * dt / dx**2
    
    print(f"网格: dx = {dx:.4f}, dt = {dt:.4f}")
    print(f"稳定性参数 r = {r:.4f} (需要 r <= 0.5)")
    
    if r > 0.5:
        print("警告：可能不稳定！")
    
    # 初始化
    x = np.linspace(0, L, Nx + 1)
    u = f(x)
    u[0], u[-1] = u_left, u_right  # 边界条件
    
    # 存储用于动画
    u_history = [u.copy()]
    t_history = [0]
    
    for n in range(Nt):
        u_new = u.copy()
        for j in range(1, Nx):
            u_new[j] = u[j] + r * (u[j+1] - 2*u[j] + u[j-1])
        
        u_new[0], u_new[-1] = u_left, u_right
        u = u_new
        
        if (n + 1) % (Nt // 20) == 0:  # 每 5%保存一次
            u_history.append(u.copy())
            t_history.append((n + 1) * dt)
    
    return x, u_history, t_history

# 初始条件
f_init = lambda x: np.sin(np.pi * x)

# 稳定情况
L, T = 1, 0.5
alpha = 0.1
Nx, Nt = 50, 1000

x, u_history, t_history = heat_ftcs(alpha, L, T, Nx, Nt, f_init)

# 绘图
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左图：数值解随时间演化
for i in range(0, len(u_history), len(u_history)//5):
    axes[0].plot(x, u_history[i], label=f't = {t_history[i]:.3f}')

# 精确解比较
x_fine = np.linspace(0, L, 200)
for t in [0, T/4, T/2, 3*T/4, T]:
    u_exact = np.sin(np.pi * x_fine) * np.exp(-alpha * np.pi**2 * t)
    axes[0].plot(x_fine, u_exact, 'k--', alpha=0.3)

axes[0].set_xlabel('x')
axes[0].set_ylabel('u')
axes[0].set_title('FTCS 方法：热方程')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 右图：时空图
t_grid = np.linspace(0, T, len(u_history))
X, T_grid = np.meshgrid(x, t_grid)
U = np.array(u_history)

c = axes[1].contourf(X, T_grid, U, levels=20, cmap='hot')
plt.colorbar(c, ax=axes[1])
axes[1].set_xlabel('x')
axes[1].set_ylabel('t')
axes[1].set_title('温度的时空分布')

plt.tight_layout()
plt.savefig('heat_ftcs.png', dpi=150, bbox_inches='tight')
plt.show()

隐式格式（ Crank-Nicolson）

FTCS 的稳定性限制可能很严格。Crank-Nicolson 格式是无条件稳定的：

这是一个三对角线性系统。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def heat_crank_nicolson(alpha, L, T, Nx, Nt, f, u_left=0, u_right=0):
    """
    热方程的 Crank-Nicolson 隐式格式
    """
    dx = L / Nx
    dt = T / Nt
    r = alpha * dt / (2 * dx**2)
    
    # 初始化
    x = np.linspace(0, L, Nx + 1)
    u = f(x)
    u[0], u[-1] = u_left, u_right
    
    # 构造三对角矩阵
    n = Nx - 1  # 内部节点数
    
    # 左端矩阵 (I + r*A)
    A_left = np.diag([1 + 2*r] * n)
    A_left += np.diag([-r] * (n-1), 1)
    A_left += np.diag([-r] * (n-1), -1)
    
    # 右端矩阵 (I - r*A)
    A_right = np.diag([1 - 2*r] * n)
    A_right += np.diag([r] * (n-1), 1)
    A_right += np.diag([r] * (n-1), -1)
    
    u_history = [u.copy()]
    t_history = [0]
    
    for step in range(Nt):
        # 右端向量
        b = A_right @ u[1:-1]
        
        # 边界修正
        b[0] += r * (u_left + u_left)  # 两个时间层的边界
        b[-1] += r * (u_right + u_right)
        
        # 求解线性系统
        u_interior = np.linalg.solve(A_left, b)
        u[1:-1] = u_interior
        
        if (step + 1) % (Nt // 20) == 0:
            u_history.append(u.copy())
            t_history.append((step + 1) * dt)
    
    return x, u_history, t_history

# 比较 FTCS 和 Crank-Nicolson
# 使用大步长测试稳定性
Nx, Nt = 50, 50  # 大 dt， r > 0.5

x_ftcs, u_ftcs, t_ftcs = heat_ftcs(alpha, L, T, Nx, Nt, f_init)
x_cn, u_cn, t_cn = heat_crank_nicolson(alpha, L, T, Nx, Nt, f_init)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# FTCS 可能不稳定
for i in range(0, min(len(u_ftcs), 5)):
    axes[0].plot(x_ftcs, u_ftcs[i], label=f't = {t_ftcs[i]:.3f}')
axes[0].set_xlabel('x')
axes[0].set_ylabel('u')
axes[0].set_title('FTCS 方法（可能不稳定）')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim([-2, 2])

# Crank-Nicolson 总是稳定
for i in range(0, len(u_cn), len(u_cn)//5):
    axes[1].plot(x_cn, u_cn[i], label=f't = {t_cn[i]:.3f}')
axes[1].set_xlabel('x')
axes[1].set_ylabel('u')
axes[1].set_title('Crank-Nicolson 方法（无条件稳定）')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('heat_stability.png', dpi=150, bbox_inches='tight')
plt.show()

波动方程

物理推导

傅里叶级数展开

考虑绷紧的弦的小振动。设 为 点在 时刻的垂直位移。

由 Newton 第二定律：

得到波动方程：

其中 是波速。

d'Alembert 解

对于无界区域的初值问题：

解为：

这表明解由左行波和右行波的叠加组成。

分离变量法

对于有界区域，边界条件：

系数由初始条件确定。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
def wave_equation_analytical(f, g, L, c, x, t, n_terms=50):
    """
    波动方程的解析解（分离变量法）
    u_tt = c ² u_xx
    u(0,t) = u(L,t) = 0
    u(x,0) = f(x), u_t(x,0) = g(x)
    """
    from scipy.integrate import quad
    
    u = np.zeros_like(x, dtype=float)
    
    for n in range(1, n_terms + 1):
        # Fourier 系数
        def f_integrand(xi):
            return f(xi) * np.sin(n * np.pi * xi / L)
        def g_integrand(xi):
            return g(xi) * np.sin(n * np.pi * xi / L)
        
        An, _ = quad(f_integrand, 0, L)
        An *= 2 / L
        
        Bn, _ = quad(g_integrand, 0, L)
        omega_n = n * np.pi * c / L
        Bn *= 2 / (L * omega_n)
        
        Xn = np.sin(n * np.pi * x / L)
        Tn = An * np.cos(omega_n * t) + Bn * np.sin(omega_n * t)
        
        u += Xn * Tn
    
    return u

# 示例：拨弦
L = 1
c = 1

# 初始位移：三角形（拨弦）
def f_pluck(x):
    if x <= 0.3:
        return x / 0.3
    else:
        return (1 - x) / 0.7

f_pluck_vec = np.vectorize(f_pluck)
g_zero = lambda x: 0  # 初始速度为 0

x = np.linspace(0, L, 100)
times = [0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75]

fig, axes = plt.subplots(2, 4, figsize=(14, 6))
axes = axes.flatten()

for i, t in enumerate(times):
    if t == 0:
        u = f_pluck_vec(x)
    else:
        u = wave_equation_analytical(f_pluck, g_zero, L, c, x, t)
    
    axes[i].plot(x, u, 'b-', linewidth=2)
    axes[i].axhline(y=0, color='k', linestyle='-', alpha=0.3)
    axes[i].set_xlim([0, L])
    axes[i].set_ylim([-1.2, 1.2])
    axes[i].set_title(f't = {t:.2f}')
    axes[i].grid(True, alpha=0.3)

plt.suptitle('波动方程：拨弦的振动', fontsize=14)
plt.tight_layout()
plt.savefig('wave_analytical.png', dpi=150, bbox_inches='tight')
plt.show()

有限差分法求解波动方程

显式格式：

整理得：

其中。

CFL 条件：（ Courant-Friedrichs-Lewy）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
def wave_fd(c, L, T, Nx, Nt, f, g, u_left=0, u_right=0):
    """
    波动方程的有限差分法
    """
    dx = L / Nx
    dt = T / Nt
    r = c * dt / dx
    
    print(f"CFL 数 r = {r:.4f} (需要 r <= 1)")
    
    if r > 1:
        print("警告：不满足 CFL 条件，可能不稳定！")
    
    x = np.linspace(0, L, Nx + 1)
    
    # 初始化
    u_prev = f(x)
    u_curr = np.zeros(Nx + 1)
    
    # 第一步用初始速度
    for j in range(1, Nx):
        u_curr[j] = u_prev[j] + dt * g(x[j]) + \
                    0.5 * r**2 * (u_prev[j+1] - 2*u_prev[j] + u_prev[j-1])
    u_curr[0], u_curr[-1] = u_left, u_right
    
    u_history = [u_prev.copy(), u_curr.copy()]
    t_history = [0, dt]
    
    # 时间推进
    for n in range(2, Nt + 1):
        u_next = np.zeros(Nx + 1)
        for j in range(1, Nx):
            u_next[j] = 2*u_curr[j] - u_prev[j] + \
                        r**2 * (u_curr[j+1] - 2*u_curr[j] + u_curr[j-1])
        u_next[0], u_next[-1] = u_left, u_right
        
        u_prev = u_curr
        u_curr = u_next
        
        if n % (Nt // 20) == 0:
            u_history.append(u_curr.copy())
            t_history.append(n * dt)
    
    return x, u_history, t_history

# 初始条件：高斯脉冲
f_gauss = lambda x: np.exp(-100*(x - 0.5)**2)
g_zero = lambda x: 0

L, T = 1, 2
c = 1
Nx, Nt = 100, 200

x, u_history, t_history = wave_fd(c, L, T, Nx, Nt, f_gauss, g_zero)

# 创建时空图
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左图：几个时刻的快照
for i in range(0, len(u_history), max(1, len(u_history)//6)):
    axes[0].plot(x, u_history[i], label=f't = {t_history[i]:.2f}')
axes[0].set_xlabel('x')
axes[0].set_ylabel('u')
axes[0].set_title('波动方程：有限差分解')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 右图：时空图
t_grid = np.array(t_history)
X, T_grid = np.meshgrid(x, t_grid)
U = np.array(u_history)

c = axes[1].contourf(X, T_grid, U, levels=20, cmap='seismic')
plt.colorbar(c, ax=axes[1])
axes[1].set_xlabel('x')
axes[1].set_ylabel('t')
axes[1].set_title('位移的时空分布')

plt.tight_layout()
plt.savefig('wave_fd.png', dpi=150, bbox_inches='tight')
plt.show()

Laplace 方程与 Poisson 方程

物理背景

Laplace 方程：

描述稳态、无源的场： - 稳态温度分布（无内热源） - 静电势（无电荷区域） - 不可压缩流的速度势

Poisson 方程：

当有源时使用（如有电荷、有热源）。

性质

极值原理：调和函数（满足 Laplace 方程的函数）在区域内部不取极值，极值只在边界上达到。

均值性质：圆心处的值等于圆周上的平均值。

有限差分法

用五点差分格式：

即：

这是一个大型稀疏线性系统。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
def laplace_2d(Nx, Ny, bc_left, bc_right, bc_bottom, bc_top):
    """
    二维 Laplace 方程 ∇² u = 0
    在单位正方形上，给定 Dirichlet 边界条件
    """
    hx = 1.0 / Nx
    hy = 1.0 / Ny
    
    # 内部节点数
    n = (Nx - 1) * (Ny - 1)
    
    # 坐标
    x = np.linspace(0, 1, Nx + 1)
    y = np.linspace(0, 1, Ny + 1)
    
    # 构造系数矩阵和右端向量
    A = np.zeros((n, n))
    b = np.zeros(n)
    
    def index(i, j):
        """将 2D 索引转换为 1D"""
        return (j - 1) * (Nx - 1) + (i - 1)
    
    for j in range(1, Ny):
        for i in range(1, Nx):
            k = index(i, j)
            
            # 主对角线
            A[k, k] = -4
            
            # 左邻居
            if i > 1:
                A[k, index(i-1, j)] = 1
            else:
                b[k] -= bc_left(y[j])
            
            # 右邻居
            if i < Nx - 1:
                A[k, index(i+1, j)] = 1
            else:
                b[k] -= bc_right(y[j])
            
            # 下邻居
            if j > 1:
                A[k, index(i, j-1)] = 1
            else:
                b[k] -= bc_bottom(x[i])
            
            # 上邻居
            if j < Ny - 1:
                A[k, index(i, j+1)] = 1
            else:
                b[k] -= bc_top(x[i])
    
    # 求解
    u_interior = np.linalg.solve(A, b)
    
    # 重构完整解
    u = np.zeros((Ny + 1, Nx + 1))
    
    # 填充内部
    for j in range(1, Ny):
        for i in range(1, Nx):
            u[j, i] = u_interior[index(i, j)]
    
    # 填充边界
    u[:, 0] = [bc_left(y[j]) for j in range(Ny + 1)]
    u[:, -1] = [bc_right(y[j]) for j in range(Ny + 1)]
    u[0, :] = [bc_bottom(x[i]) for i in range(Nx + 1)]
    u[-1, :] = [bc_top(x[i]) for i in range(Nx + 1)]
    
    return x, y, u

# 边界条件
bc_left = lambda y: 0
bc_right = lambda y: 0
bc_bottom = lambda x: 0
bc_top = lambda x: np.sin(np.pi * x)  # 顶部加热

Nx, Ny = 50, 50
x, y, u = laplace_2d(Nx, Ny, bc_left, bc_right, bc_bottom, bc_top)

# 精确解（分离变量）
X, Y = np.meshgrid(x, y)
u_exact = np.sin(np.pi * X) * np.sinh(np.pi * Y) / np.sinh(np.pi)

fig = plt.figure(figsize=(14, 5))

# 数值解
ax1 = fig.add_subplot(131)
c1 = ax1.contourf(X, Y, u, levels=20, cmap='hot')
plt.colorbar(c1, ax=ax1)
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_title('数值解')
ax1.set_aspect('equal')

# 精确解
ax2 = fig.add_subplot(132)
c2 = ax2.contourf(X, Y, u_exact, levels=20, cmap='hot')
plt.colorbar(c2, ax=ax2)
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_title('精确解')
ax2.set_aspect('equal')

# 误差
ax3 = fig.add_subplot(133)
c3 = ax3.contourf(X, Y, abs(u - u_exact), levels=20, cmap='viridis')
plt.colorbar(c3, ax=ax3)
ax3.set_xlabel('x')
ax3.set_ylabel('y')
ax3.set_title('绝对误差')
ax3.set_aspect('equal')

plt.tight_layout()
plt.savefig('laplace_2d.png', dpi=150, bbox_inches='tight')
plt.show()

print(f"最大误差: {np.max(abs(u - u_exact)):.6e}")

迭代方法

对于大规模问题，直接求解太慢。常用迭代方法：

Jacobi 迭代：

Gauss-Seidel 迭代：

SOR（逐次超松弛）：

最优松弛因子约为。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
def laplace_sor(Nx, Ny, bc, f=None, omega=None, tol=1e-6, max_iter=10000):
    """
    SOR 方法求解 Laplace/Poisson 方程
    bc: 字典 {'left': func, 'right': func, 'bottom': func, 'top': func}
    f: 源项函数（ Poisson 方程），默认为 0（ Laplace 方程）
    """
    h = 1.0 / Nx
    
    if omega is None:
        omega = 2 / (1 + np.sin(np.pi * h))  # 最优 SOR 因子
    
    x = np.linspace(0, 1, Nx + 1)
    y = np.linspace(0, 1, Ny + 1)
    
    # 初始化
    u = np.zeros((Ny + 1, Nx + 1))
    
    # 边界条件
    u[:, 0] = [bc['left'](y[j]) for j in range(Ny + 1)]
    u[:, -1] = [bc['right'](y[j]) for j in range(Ny + 1)]
    u[0, :] = [bc['bottom'](x[i]) for i in range(Nx + 1)]
    u[-1, :] = [bc['top'](x[i]) for i in range(Nx + 1)]
    
    residuals = []
    
    for iteration in range(max_iter):
        max_change = 0
        
        for j in range(1, Ny):
            for i in range(1, Nx):
                old_u = u[j, i]
                
                # Gauss-Seidel 更新值
                f_ij = f(x[i], y[j]) if f else 0
                new_u = 0.25 * (u[j, i+1] + u[j, i-1] + 
                               u[j+1, i] + u[j-1, i] - h**2 * f_ij)
                
                # SOR 更新
                u[j, i] = (1 - omega) * old_u + omega * new_u
                
                max_change = max(max_change, abs(u[j, i] - old_u))
        
        residuals.append(max_change)
        
        if max_change < tol:
            print(f"SOR 收敛于第 {iteration + 1} 次迭代")
            break
    
    return x, y, u, residuals

# 测试 SOR
bc = {
    'left': lambda y: 0,
    'right': lambda y: 0,
    'bottom': lambda x: 0,
    'top': lambda x: np.sin(np.pi * x)
}

Nx, Ny = 50, 50
x, y, u_sor, residuals = laplace_sor(Nx, Ny, bc)

# 收敛历史
plt.figure(figsize=(10, 5))
plt.semilogy(residuals)
plt.xlabel('迭代次数')
plt.ylabel('最大变化量')
plt.title('SOR 方法的收敛历史')
plt.grid(True, alpha=0.3)
plt.savefig('sor_convergence.png', dpi=150, bbox_inches='tight')
plt.show()

数值方法的稳定性分析

Von Neumann 稳定性分析

基本思想：将误差展开为 Fourier 模式，分析每个模式的增长因子。

设误差，代入差分格式，求解 的条件。

例：热方程 FTCS 格式

代入：

稳定性要求：

当 时，。

要使，需要。

CFL 条件

对于波动方程， CFL 条件说：数值信息传播速度必须大于等于物理信息传播速度。

物理波速： 数值信息传播速度：CFL 条件：，即。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
# 演示 CFL 条件的重要性
c = 1
L, T = 1, 0.5
Nx = 50

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# 情况 1： r < 1（稳定）
Nt_stable = 100
r_stable = c * (T/Nt_stable) / (L/Nx)
x1, u1, t1 = wave_fd(c, L, T, Nx, Nt_stable, f_gauss, g_zero)
axes[0].plot(x1, u1[-1], 'b-')
axes[0].set_title(f'稳定： r = {r_stable:.2f}')
axes[0].set_ylim([-1.5, 1.5])
axes[0].grid(True, alpha=0.3)

# 情况 2： r = 1（边界稳定，实际上是精确的！）
Nt_boundary = 50
r_boundary = c * (T/Nt_boundary) / (L/Nx)
x2, u2, t2 = wave_fd(c, L, T, Nx, Nt_boundary, f_gauss, g_zero)
axes[1].plot(x2, u2[-1], 'g-')
axes[1].set_title(f'边界稳定： r = {r_boundary:.2f}')
axes[1].set_ylim([-1.5, 1.5])
axes[1].grid(True, alpha=0.3)

# 情况 3： r > 1（不稳定）
Nt_unstable = 30
r_unstable = c * (T/Nt_unstable) / (L/Nx)
try:
    x3, u3, t3 = wave_fd(c, L, T, Nx, Nt_unstable, f_gauss, g_zero)
    # 限制显示范围
    u3_clipped = np.clip(u3[-1], -10, 10)
    axes[2].plot(x3, u3_clipped, 'r-')
except:
    pass
axes[2].set_title(f'不稳定： r = {r_unstable:.2f}')
axes[2].set_ylim([-1.5, 1.5])
axes[2].grid(True, alpha=0.3)

plt.suptitle('波动方程： CFL 条件的影响')
plt.tight_layout()
plt.savefig('cfl_condition.png', dpi=150, bbox_inches='tight')
plt.show()

高级话题

谱方法

思想：用全局基函数（如 Fourier 级数、 Chebyshev 多项式）展开解。

优点：对于光滑解，可以达到指数收敛（ spectral accuracy）。

缺点：处理复杂边界困难；不适合不光滑解。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
def heat_spectral(alpha, L, T, N, Nt, f):
    """
    热方程的谱方法（ Fourier）
    周期边界条件
    """
    x = np.linspace(0, L, N, endpoint=False)
    dx = L / N
    dt = T / Nt
    
    # 初始条件的 FFT
    u_hat = np.fft.fft(f(x))
    
    # 波数
    k = np.fft.fftfreq(N, dx) * 2 * np.pi
    
    u_history = [f(x)]
    t_history = [0]
    
    # 时间演化
    for n in range(1, Nt + 1):
        # 精确时间演化（对于线性方程）
        u_hat = u_hat * np.exp(-alpha * k**2 * dt)
        
        if n % (Nt // 20) == 0:
            u = np.real(np.fft.ifft(u_hat))
            u_history.append(u)
            t_history.append(n * dt)
    
    return x, u_history, t_history

# 比较谱方法和有限差分
L, T = 2*np.pi, 1
alpha = 0.1
f_smooth = lambda x: np.sin(x) + 0.5*np.sin(3*x)

# 谱方法
x_spec, u_spec, t_spec = heat_spectral(alpha, L, T, 32, 100, f_smooth)

# 精确解
t_final = T
u_exact = np.sin(x_spec)*np.exp(-alpha*t_final) + 0.5*np.sin(3*x_spec)*np.exp(-9*alpha*t_final)

plt.figure(figsize=(10, 5))
plt.plot(x_spec, u_spec[-1], 'b.-', label='谱方法 (N=32)')
plt.plot(x_spec, u_exact, 'r--', label='精确解')
plt.xlabel('x')
plt.ylabel('u')
plt.title('谱方法 vs 精确解')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('spectral_method.png', dpi=150, bbox_inches='tight')
plt.show()

print(f"谱方法误差: {np.max(abs(u_spec[-1] - u_exact)):.2e}")

有限元方法简介

有限元方法（ FEM）是工程中最常用的 PDE 数值方法。

基本步骤：

  1. 弱形式：将 PDE 乘以测试函数并积分
  2. 离散化：用有限维空间近似
  3. 组装：构造刚度矩阵和载荷向量
  4. 求解：解线性系统

优点： - 处理复杂几何灵活 - 理论基础完善 - 可以处理非均匀网格

多重网格方法

问题：迭代方法（如 Jacobi 、 Gauss-Seidel）对低频误差收敛慢。

多重网格思想： 1. 在细网格上迭代几步（消除高频误差） 2. 将残差限制到粗网格 3. 在粗网格上求解（低频误差在粗网格上变成高频） 4. 将粗网格修正延拓回细网格 5. 重复

复杂度：！（最优）

实际应用

热扩散：冷却鳍片

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def cooling_fin(k, h_conv, T_base, T_ambient, L, A, P, Nx):
    """
    冷却鳍片的稳态温度分布
    -kA d ² T/dx ² + hP(T - T_ambient) = 0
    T(0) = T_base, dT/dx(L) = 0 (绝热端)
    """
    dx = L / Nx
    m2 = h_conv * P / (k * A)  # 特征参数
    
    # 有限差分
    T = np.zeros(Nx + 1)
    
    # 系数矩阵
    A_mat = np.zeros((Nx + 1, Nx + 1))
    b = np.zeros(Nx + 1)
    
    # 边界条件 T(0) = T_base
    A_mat[0, 0] = 1
    b[0] = T_base
    
    # 内部节点
    for i in range(1, Nx):
        A_mat[i, i-1] = 1
        A_mat[i, i] = -2 - m2 * dx**2
        A_mat[i, i+1] = 1
        b[i] = -m2 * dx**2 * T_ambient
    
    # 边界条件 dT/dx(L) = 0 (用向后差分)
    A_mat[Nx, Nx] = 1
    A_mat[Nx, Nx-1] = -1
    b[Nx] = 0
    
    T = np.linalg.solve(A_mat, b)
    x = np.linspace(0, L, Nx + 1)
    
    # 精确解
    m = np.sqrt(m2)
    T_exact = T_ambient + (T_base - T_ambient) * np.cosh(m*(L-x)) / np.cosh(m*L)
    
    return x, T, T_exact

# 参数
k = 200  # 热导率 W/(m · K)（铝）
h = 25   # 对流换热系数 W/(m ²· K)
T_base = 100  # 基底温度 ° C
T_ambient = 25  # 环境温度 ° C
L = 0.1  # 长度 m
width, thickness = 0.05, 0.005  # 宽度和厚度 m
A = width * thickness  # 截面积
P = 2 * (width + thickness)  # 周长

x, T_num, T_exact = cooling_fin(k, h, T_base, T_ambient, L, A, P, Nx=50)

plt.figure(figsize=(10, 5))
plt.plot(x*100, T_num, 'b.-', label='数值解')
plt.plot(x*100, T_exact, 'r--', label='精确解')
plt.xlabel('位置 (cm)')
plt.ylabel('温度 (° C)')
plt.title('冷却鳍片的温度分布')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('cooling_fin.png', dpi=150, bbox_inches='tight')
plt.show()

二维稳态传热

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def heat_2d_source(Nx, Ny, k, Q, T_boundary):
    """
    二维稳态传热，带内热源
    -k ∇² T = Q
    """
    hx = 1.0 / Nx
    hy = 1.0 / Ny
    
    x = np.linspace(0, 1, Nx + 1)
    y = np.linspace(0, 1, Ny + 1)
    X, Y = np.meshgrid(x, y)
    
    # 热源分布（圆形区域）
    Q_field = np.zeros((Ny + 1, Nx + 1))
    center = (0.5, 0.5)
    radius = 0.2
    mask = (X - center[0])**2 + (Y - center[1])**2 < radius**2
    Q_field[mask] = Q
    
    # 用 SOR 求解
    T = np.ones((Ny + 1, Nx + 1)) * T_boundary
    
    omega = 1.8
    
    for _ in range(5000):
        T_old = T.copy()
        
        for j in range(1, Ny):
            for i in range(1, Nx):
                T_new = 0.25 * (T[j, i+1] + T[j, i-1] + 
                               T[j+1, i] + T[j-1, i] + 
                               hx**2 * Q_field[j, i] / k)
                T[j, i] = (1 - omega) * T[j, i] + omega * T_new
        
        if np.max(abs(T - T_old)) < 1e-6:
            break
    
    return X, Y, T, Q_field

# 参数
k = 1  # 热导率
Q = 100  # 热源强度
T_boundary = 0  # 边界温度

X, Y, T, Q_field = heat_2d_source(50, 50, k, Q, T_boundary)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 热源分布
c1 = axes[0].contourf(X, Y, Q_field, levels=20, cmap='Reds')
plt.colorbar(c1, ax=axes[0])
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_title('热源分布')
axes[0].set_aspect('equal')

# 温度分布
c2 = axes[1].contourf(X, Y, T, levels=20, cmap='hot')
plt.colorbar(c2, ax=axes[1])
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].set_title('稳态温度分布')
axes[1].set_aspect('equal')

plt.tight_layout()
plt.savefig('heat_2d_source.png', dpi=150, bbox_inches='tight')
plt.show()

波动方程：鼓膜振动

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def drum_vibration(Nx, Ny, c, T, Nt, f_initial):
    """
    圆形鼓膜的振动（在正方形网格上近似）
    使用掩模处理圆形边界
    """
    dx = 2.0 / Nx  # 域[-1,1]×[-1,1]
    dy = 2.0 / Ny
    dt = T / Nt
    
    r = c * dt / dx
    
    x = np.linspace(-1, 1, Nx + 1)
    y = np.linspace(-1, 1, Ny + 1)
    X, Y = np.meshgrid(x, y)
    
    # 圆形掩模
    mask = X**2 + Y**2 <= 1
    
    # 初始化
    u_prev = f_initial(X, Y) * mask
    u_curr = u_prev.copy()  # 假设初始速度为 0
    
    u_history = [u_prev.copy()]
    
    for n in range(Nt):
        u_next = np.zeros_like(u_curr)
        
        for j in range(1, Ny):
            for i in range(1, Nx):
                if mask[j, i]:
                    u_next[j, i] = (2*u_curr[j, i] - u_prev[j, i] +
                                   r**2 * (u_curr[j, i+1] + u_curr[j, i-1] +
                                          u_curr[j+1, i] + u_curr[j-1, i] - 
                                          4*u_curr[j, i]))
        
        # 边界条件（圆周上 u=0）
        u_next[~mask] = 0
        
        u_prev = u_curr
        u_curr = u_next
        
        if n % (Nt // 8) == 0:
            u_history.append(u_curr.copy())
    
    return X, Y, u_history

# 初始条件：中心凸起
f_init = lambda X, Y: np.maximum(0, 0.5 - 2*np.sqrt(X**2 + Y**2)) * (X**2 + Y**2 <= 1)

c = 1
T = 3
X, Y, u_hist = drum_vibration(50, 50, c, T, 300, f_init)

# 绘制振动过程
fig, axes = plt.subplots(2, 4, figsize=(16, 8))
axes = axes.flatten()

for i, u in enumerate(u_hist):
    ax = axes[i]
    c = ax.contourf(X, Y, u, levels=20, cmap='seismic', vmin=-0.3, vmax=0.3)
    circle = plt.Circle((0, 0), 1, fill=False, color='k', linewidth=2)
    ax.add_patch(circle)
    ax.set_aspect('equal')
    ax.set_title(f'帧 {i}')
    ax.set_xlim([-1.2, 1.2])
    ax.set_ylim([-1.2, 1.2])

plt.suptitle('圆形鼓膜的振动', fontsize=14)
plt.tight_layout()
plt.savefig('drum_vibration.png', dpi=150, bbox_inches='tight')
plt.show()

练习题

基础题

练习 1：对方程 进行分类（双曲/抛物/椭圆）。

练习 2：用分离变量法求解：$

$

练习 3：实现热方程的 BTCS（向后时间中心空间）格式，验证它是无条件稳定的。

练习 4：用有限差分法求解 Poisson 方程：

在单位正方形上，边界为 0 。精确解是什么？

练习 5：验证 d'Alembert 公式对于，，。

中级题

练习 6：实现波动方程的 Newmark-β方法，比较不同, 参数的稳定性。

练习 7：用 ADI（交替方向隐式）方法求解二维热方程。

练习 8：对于不规则边界（如 L 形区域），实现有限差分法求解 Laplace 方程。

练习 9：实现热方程的谱方法（用 FFT），比较与有限差分的精度。

练习 10：研究波动方程在不同边界条件下的行为： - 固定端： - 自由端：

高级题

练习 11：实现简单的多重网格方法（ V-cycle）求解 Poisson 方程，验证 复杂度。

练习 12：求解非线性热方程：

研究解的行为。

练习 13：用有限差分法求解 Schr ö dinger 方程：

观察波包的演化。

练习 14：实现 Burgers 方程的数值解：

观察激波的形成。

练习 15：用边界元方法（ BEM）求解外部 Laplace 问题，与有限差分比较。

总结

本章我们进入了偏微分方程的世界：

方程类型 物理模型 典型数值方法 稳定性条件
抛物型（热方程） 扩散、热传导 FTCS, CN (FTCS)
双曲型（波动方程） 波传播、振动 中心差分 CFL:
椭圆型（ Laplace） 稳态、平衡 迭代法、直接法 总是稳定

关键概念：

  1. 分类：判别式 决定方程类型
  2. 适定性：存在性、唯一性、稳定性
  3. 稳定性分析： Von Neumann 分析、 CFL 条件
  4. 数值方法：有限差分、有限元、谱方法

实践建议：

  1. 先从简单的显式格式开始
  2. 检查 CFL 条件
  3. 验证收敛阶数
  4. 对于刚性问题使用隐式方法
  5. 大规模问题考虑迭代方法或多重网格

这一章标志着我们从 ODE 到 PDE 的跨越。 PDE 是一个广阔的领域，我们只是触及了冰山一角。希望这个引论能为你进一步学习计算数学、流体力学、电磁学等领域打下基础！

  • 本文标题：常微分方程（十三）偏微分方程引论
  • 本文作者：Chen Kai
  • 创建时间：2019-06-09 14:15:00
  • 本文链接：https://www.chenk.top/%E5%B8%B8%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%EF%BC%88%E5%8D%81%E4%B8%89%EF%BC%89%E5%81%8F%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%E5%BC%95%E8%AE%BA/
  • 版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
 评论