用Python和Matplotlib可视化理解向量场：从曲线积分到环量通量（附完整代码）

第一次接触向量场的概念时，那些抽象的数学公式总让我感到困惑——箭头在空间中的分布究竟代表什么？为什么曲线积分能描述"做功"？环量和通量又该如何直观理解？直到我开始用Python将这些概念可视化，一切才变得清晰起来。本文将带你用NumPy和Matplotlib，通过代码实现从二维到三维的向量场可视化，并动态演示曲线积分、环量和通量的计算过程。

1. 环境准备与基础向量场绘制

在开始前，确保已安装以下Python库：

pip install numpy matplotlib

让我们从一个简单的二维向量场开始。考虑风场模型，假设风速向量在平面上的分布为F(x,y) = (-y, x)。用Matplotlib绘制这个向量场：

import numpy as np
import matplotlib.pyplot as plt

# 创建网格
x = np.linspace(-5, 5, 10)
y = np.linspace(-5, 5, 10)
X, Y = np.meshgrid(x, y)

# 定义向量场
U = -Y  # x方向分量
V = X   # y方向分量

# 绘制向量场
plt.figure(figsize=(8, 6))
plt.quiver(X, Y, U, V, color='blue', scale=30)
plt.title('二维向量场示例: F(x,y) = (-y, x)')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

这段代码会生成一个旋转的向量场图案。 关键参数说明 ：

• np.linspace ：创建均匀分布的点
• np.meshgrid ：生成坐标网格
• plt.quiver ：绘制箭头， scale 控制箭头大小

提示：尝试修改U和V的表达式，观察不同向量场的形态差异，这是理解向量场特性的第一步。

2. 曲线积分的可视化实现

曲线积分计算的是向量场沿某条路径的"累积效应"。物理上可以理解为力场中移动物体所做的功。让我们实现一个二维曲线积分的可视化示例。

2.1 定义路径和向量场

# 定义参数曲线 r(t) = [t, t^2], t ∈ [0,1]
t = np.linspace(0, 1, 100)
x_path = t
y_path = t**2

# 定义向量场 F = [y, x]
def vector_field(x, y):
    return y, x

# 计算路径上各点的向量场值
Fx, Fy = vector_field(x_path, y_path)

2.2 绘制路径和向量场

plt.figure(figsize=(10, 6))

# 绘制向量场背景
x_grid = np.linspace(0, 1, 10)
y_grid = np.linspace(0, 1, 10)
X, Y = np.meshgrid(x_grid, y_grid)
U, V = vector_field(X, Y)
plt.quiver(X, Y, U, V, color='lightblue', scale=30)

# 绘制路径
plt.plot(x_path, y_path, 'r-', linewidth=2)

# 在路径上采样几个点显示向量
sample_indices = [10, 30, 50, 70, 90]
for i in sample_indices:
    plt.quiver(x_path[i], y_path[i], Fx[i], Fy[i], 
               color='red', scale=30, width=0.005)

plt.title('曲线积分可视化: 路径与向量场')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

2.3 计算曲线积分

曲线积分的计算公式为：∫F·dr = ∫F(r(t))·r'(t)dt

# 计算路径的导数
dx_dt = np.gradient(x_path, t)
dy_dt = np.gradient(y_path, t)

# 计算点积 F·dr/dt
integrand = Fx * dx_dt + Fy * dy_dt

# 数值积分
work = np.trapz(integrand, t)
print(f"曲线积分值(做功): {work:.4f}")

注意： np.gradient 用于数值计算导数， np.trapz 实现梯形法数值积分。对于简单函数，也可以解析求导。

3. 环量与旋度的三维可视化

环量是向量场沿闭合曲线的曲线积分，与旋度密切相关。让我们创建一个三维向量场并可视化其旋度。

3.1 三维向量场定义

from mpl_toolkits.mplot3d import Axes3D

# 创建三维网格
x = np.linspace(-2, 2, 8)
y = np.linspace(-2, 2, 8)
z = np.linspace(-2, 2, 8)
X, Y, Z = np.meshgrid(x, y, z)

# 定义三维向量场 F = [-y, x, z]
U = -Y
V = X
W = Z

# 绘制三维向量场
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.quiver(X, Y, Z, U, V, W, length=0.2, normalize=True)
ax.set_title('三维向量场 F = [-y, x, z]')
plt.show()

3.2 计算并可视化旋度

旋度计算公式：∇×F = (∂F₃/∂y - ∂F₂/∂z, ∂F₁/∂z - ∂F₃/∂x, ∂F₂/∂x - ∂F₁/∂y)

# 计算旋度各分量
dW_dy = np.gradient(W, axis=1)
dV_dz = np.gradient(V, axis=2)
curl_x = dW_dy - dV_dz

dU_dz = np.gradient(U, axis=2)
dW_dx = np.gradient(W, axis=0)
curl_y = dU_dz - dW_dx

dV_dx = np.gradient(V, axis=0)
dU_dy = np.gradient(U, axis=1)
curl_z = dV_dx - dU_dy

# 可视化旋度场
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# 绘制原始向量场(透明)
ax.quiver(X, Y, Z, U, V, W, length=0.2, 
          color='blue', alpha=0.3, normalize=True)

# 绘制旋度场
ax.quiver(X, Y, Z, curl_x, curl_y, curl_z, length=0.2,
          color='red', normalize=True)

ax.set_title('蓝色: 原始向量场 | 红色: 旋度场')
plt.show()

旋度的物理意义 ：旋度矢量方向表示旋转轴，大小表示旋转强度。在上面的例子中，可以看到z方向有明显的旋度分量。

4. 通量与散度的交互式演示

通量描述向量场通过曲面的"流量"，而散度则表示某点的"源强度"。让我们创建一个交互式示例来演示这些概念。

4.1 二维散度计算

# 定义新向量场 F = [x, y]
def div_vector_field(x, y):
    return x, y

# 计算散度 (∂F₁/∂x + ∂F₂/∂y)
def compute_divergence(Fx, Fy, x, y):
    dFx_dx = np.gradient(Fx, axis=1)
    dFy_dy = np.gradient(Fy, axis=0)
    return dFx_dx + dFy_dy

# 创建网格
x = np.linspace(-2, 2, 20)
y = np.linspace(-2, 2, 20)
X, Y = np.meshgrid(x, y)

# 计算向量场和散度
U, V = div_vector_field(X, Y)
div = compute_divergence(U, V, X, Y)

# 绘制
plt.figure(figsize=(10, 8))
plt.quiver(X, Y, U, V, color='blue', scale=30)
plt.contourf(X, Y, div, levels=20, cmap='RdBu')
plt.colorbar(label='散度值')
plt.title('向量场与散度热图 (蓝色箭头: 向量场)')
plt.show()

4.2 通量计算示例

考虑圆形闭合曲线，计算向量场通过它的通量：

# 定义圆形路径
theta = np.linspace(0, 2*np.pi, 100)
x_circle = 1.5 * np.cos(theta)
y_circle = 1.5 * np.sin(theta)

# 计算路径上的向量场
Fx_circle, Fy_circle = div_vector_field(x_circle, y_circle)

# 计算法向量 (单位圆的外法向量就是位置向量归一化)
normals = np.array([x_circle, y_circle]).T
normals = normals / np.linalg.norm(normals, axis=1)[:, np.newaxis]

# 计算通量 (F·n ds)
ds = 1.5 * (theta[1] - theta[0])  # 弧长微元
flux = np.sum(Fx_circle * normals[:, 0] + Fy_circle * normals[:, 1]) * ds

print(f"通过圆形边界的通量: {flux:.4f}")

散度定理验证 ：根据散度定理，通量应等于散度在圆内的积分。我们可以数值验证这一点：

# 创建圆形掩模
r = np.sqrt(X**2 + Y**2)
circle_mask = r <= 1.5

# 计算散度在圆内的积分
integral_div = np.sum(div[circle_mask]) * (x[1]-x[0])**2
print(f"散度在圆内的积分: {integral_div:.4f}")

5. 进阶应用：电磁场可视化案例

作为综合应用，让我们模拟一个简单的电磁场场景。考虑两个点电荷产生的电场：

5.1 点电荷电场定义

def electric_field(q, pos, X, Y):
    """计算点电荷产生的电场"""
    dx = X - pos[0]
    dy = Y - pos[1]
    r = np.sqrt(dx**2 + dy**2)
    Ex = q * dx / r**3
    Ey = q * dy / r**3
    return Ex, Ey

# 创建网格
x = np.linspace(-3, 3, 20)
y = np.linspace(-3, 3, 20)
X, Y = np.meshgrid(x, y)

# 定义两个点电荷 (+1在(-1,0), -1在(1,0))
Ex1, Ey1 = electric_field(1, [-1, 0], X, Y)
Ex2, Ey2 = electric_field(-1, [1, 0], X, Y)
Ex_total = Ex1 + Ex2
Ey_total = Ey1 + Ey2

# 绘制电场线
plt.figure(figsize=(10, 8))
plt.streamplot(X, Y, Ex_total, Ey_total, color='blue', density=2)
plt.scatter([-1, 1], [0, 0], c=['red', 'green'], s=200)
plt.title('两个点电荷的电场线 (红: 正电荷, 绿: 负电荷)')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

5.2 计算电通量

选择圆形高斯面，验证高斯定律：

# 定义高斯面 (圆心在原点，半径2)
theta = np.linspace(0, 2*np.pi, 100)
x_gauss = 2 * np.cos(theta)
y_gauss = 2 * np.sin(theta)

# 计算电场在高斯面上的值
Ex_gauss, Ey_gauss = electric_field(1, [-1, 0], x_gauss, y_gauss)
Ex_gauss += electric_field(-1, [1, 0], x_gauss, y_gauss)[0]
Ey_gauss += electric_field(-1, [1, 0], x_gauss, y_gauss)[1]

# 计算通量
normals = np.array([x_gauss, y_gauss]).T / 2  # 单位法向量
flux = np.sum(Ex_gauss * normals[:, 0] + Ey_gauss * normals[:, 1]) * (2 * np.pi / 100)

print(f"通过高斯面的电通量: {flux:.4f}")
print(f"高斯定律预测值 (Q_enclosed/ε0): {0.0}")  # 因为总电荷为+1-1=0

提示：尝试修改电荷位置和大小，观察电场线和通量的变化。对于更复杂的电荷分布，可以扩展此代码计算电势和电场能量密度。