从零开始，用Python和Matplotlib可视化库仑定律与电场线（附完整代码）
从零开始用Python和Matplotlib可视化库仑定律与电场线
当物理公式遇上Python代码，抽象的电场概念突然变得触手可及。本文将带您用不到100行代码，构建一个完整的静电场可视化系统——从单个点电荷的辐射状电场线，到复杂电荷分布的动态力场模拟，让库仑定律的数学之美通过编程语言跃然屏上。
1. 环境准备与基础原理
1.1 安装必要的Python库
在开始前，请确保已安装以下库（可通过pip安装）：
pip install numpy matplotlib scipy库仑定律的数学表达为：
F = k * q1 * q2 / r²其中：
k为静电力常数（8.9875×10⁹ N·m²/C²）q1,q2为点电荷量（单位：库仑）r为电荷间距离（单位：米）
1.2 电场强度计算原理
单个点电荷在空间某点产生的电场强度：
def electric_field(q, r0, x, y): r = np.sqrt((x - r0[0])**2 + (y - r0[1])**2) return q * (x - r0[0]) / r**3, q * (y - r0[1]) / r**3注意：当r接近零时需要特殊处理以避免除零错误
2. 单点电荷电场可视化
2.1 创建二维网格
import numpy as np import matplotlib.pyplot as plt x = np.linspace(-5, 5, 100) y = np.linspace(-5, 5, 100) X, Y = np.meshgrid(x, y)2.2 计算并绘制电场线
# 正电荷位于原点 Ex, Ey = electric_field(1, (0, 0), X, Y) plt.figure(figsize=(8, 6)) plt.streamplot(X, Y, Ex, Ey, density=1.5, color='b', linewidth=1) plt.scatter([0], [0], c='red', s=200) # 绘制正电荷 plt.title('单个正点电荷的电场线分布') plt.xlabel('x (m)') plt.ylabel('y (m)') plt.grid() plt.show()执行结果将显示典型的辐射状电场线分布，直观展示：
- 电场线从正电荷向外发散
- 场强随距离增大而减弱（通过线密度表现）
3. 多点电荷系统模拟
3.1 多电荷场强叠加
def multi_charges_field(charges, X, Y): Ex_total = np.zeros_like(X) Ey_total = np.zeros_like(Y) for q, pos in charges: Ex, Ey = electric_field(q, pos, X, Y) Ex_total += Ex Ey_total += Ey return Ex_total, Ey_total3.2 电偶极子案例
charges = [(1, (-1, 0)), (-1, (1, 0))] # 电偶极子配置 Ex, Ey = multi_charges_field(charges, X, Y) plt.figure(figsize=(10, 8)) plt.streamplot(X, Y, Ex, Ey, density=2, color=np.sqrt(Ex**2 + Ey**2), cmap='viridis') plt.scatter([-1, 1], [0, 0], c=['red', 'blue'], s=200) plt.colorbar(label='电场强度') plt.title('电偶极子电场线分布') plt.show()关键观察点：
- 正负电荷间的电场线连接
- 远离系统的场线分布特征
- 中心区域的场强变化梯度
4. 高级可视化技巧
4.1 电势与等势面绘制
电势计算公式：
def electric_potential(q, r0, x, y): r = np.sqrt((x - r0[0])**2 + (y - r0[1])**2) return q / r绘制等势线：
from matplotlib import cm phi = sum(electric_potential(q, pos, X, Y) for q, pos in charges) plt.contour(X, Y, phi, 20, cmap=cm.plasma) plt.contourf(X, Y, phi, 20, alpha=0.3, cmap=cm.plasma) plt.colorbar(label='电势 (V)')4.2 交互式3D电势图
from mpl_toolkits.mplot3d import Axes3D fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111, projection='3d') ax.plot_surface(X, Y, phi, cmap=cm.viridis, alpha=0.8) ax.set_zlabel('电势 (V)')5. 实战应用：复杂电荷分布
5.1 均匀带电圆环
def ring_charge_field(Q, a, X, Y, n=100): theta = np.linspace(0, 2*np.pi, n) dq = Q/n Ex = np.zeros_like(X) Ey = np.zeros_like(Y) for th in theta: pos = (a*np.cos(th), a*np.sin(th)) dEx, dEy = electric_field(dq, pos, X, Y) Ex += dEx Ey += dEy return Ex, Ey5.2 可视化对比
| 电荷分布类型 | 代码复杂度 | 典型特征 |
|---|---|---|
| 单点电荷 | 简单 | 径向对称 |
| 电偶极子 | 中等 | 连接场线 |
| 带电圆环 | 复杂 | 轴向对称 |
# 比较不同半径圆环的场分布 for a in [1, 2, 3]: Ex, Ey = ring_charge_field(1, a, X, Y) plt.streamplot(X, Y, Ex, Ey, density=1.5)6. 性能优化技巧
当处理大量电荷时，可采用矢量化计算：
def vectorized_field(charges, X, Y): Ex = np.zeros_like(X) Ey = np.zeros_like(Y) for q, (x0, y0) in charges: r = np.sqrt((X - x0)**2 + (Y - y0)**2) with np.errstate(divide='ignore', invalid='ignore'): Ex += q * (X - x0) / r**3 Ey += q * (Y - y0) / r**3 Ex = np.nan_to_num(Ex) Ey = np.nan_to_num(Ey) return Ex, Ey提示：对于超大规模计算，可考虑使用Numba加速或GPU计算
7. 扩展应用：电场中的粒子运动
模拟带电粒子在电场中的轨迹：
def particle_trace(q, m, initial_pos, initial_v, E_func, t_max=10, dt=0.01): t = np.arange(0, t_max, dt) r = np.zeros((len(t), 2)) v = np.zeros_like(r) r[0] = initial_pos v[0] = initial_v for i in range(1, len(t)): Ex, Ey = E_func(r[i-1][0], r[i-1][1]) a = np.array([Ex, Ey]) * q / m v[i] = v[i-1] + a * dt r[i] = r[i-1] + v[i] * dt return r实际项目中，这种可视化方法曾帮助我快速验证粒子加速器的电场设计——通过调整代码中的电荷分布参数，能在几分钟内观察到不同配置下的场线变化，比传统手工绘图效率提升数十倍。