不知道你有没有遇到过这样的场景：上课讲电偶极子，书上画的等势线挺漂亮，可自己在草稿纸上推了半天，脑子里还是那块球形电场绕来绕去的样子转不过弯来。我当时学电磁学就是这样，公式能背，题会算，但要让我说说电偶极子在空间里到底长什么样，电场线从哪出发到哪结束，脑子里一片空白。后来我干脆写了个MATLAB仿真，把电位和电场全部画出来，才算是把电偶极子彻底“看”明白了。这篇文章我就把这个思路和完整源码分享出来，不光是贴代码，连每一步的数学原理、画图参数、常见报错都一起聊清楚，保证你能直接运行出图，也能看懂每行代码背后的道理。
这套内容适合三类人：正在学大学物理或电磁场理论的学生，想用MATLAB做电磁场可视化但不知道怎么下手的初学者，以及要做科研演示、课程报告的朋友。你不需要有多深的MATLAB功底，只要会基本语法，跟着思路走就能跑通。
1. 电偶极子模拟的整体设计与数学基础
1.1 电偶极子的物理模型回顾
电偶极子说白了就是两个电荷量相等、符号相反的点电荷，距离靠得很近。一个带正电+q，一个带负电-q，它们之间的距离是d。这个系统最重要的物理量叫电偶极矩，记作p，大小是 p = q·d，方向从负电荷指向正电荷。很多教材里也会用p = q·a来表示，其中a是从负电荷到正电荷的位移矢量。
之所以要把“间距很小”这个条件明确出来，是因为电偶极子的很多性质都是在小间距近似下才成立的。比如说远处某点的电位，就近似为：
[ V(\mathbf{r}) = \frac{1}{4\pi\varepsilon_0} \cdot \frac{\mathbf{p} \cdot \hat{\mathbf{r}}}{r^2} ]
这个式子的含义很直观：某点的电位与偶极矩在该点方向上的投影成正比，与距离的平方成反比。对比单个点电荷的电位 V = q/(4πε₀r)，偶极子的电位随距离衰减得更快，这是很多教材强调的重点。
但注意，这个公式是远场近似。如果场点离偶极子很近，比如距离和电荷间距d一个量级，那就不能再用这个近似公式了，必须用两个点电荷电位的直接叠加：
[ V(x,y) = \frac{1}{4\pi\varepsilon_0} \left( \frac{q}{r_+} - \frac{q}{r_-} \right) ]
其中r₊是场点到正电荷的距离，r₋是场点到负电荷的距离。我们仿真时优先采用这种精确表达式，因为它能画出近场区域的真实分布，远场自然也会退化为近似公式。
具体到二维平面（比如x-y平面），假设正电荷在 (d/2, 0)，负电荷在 (-d/2, 0)，那么空间任意一点 (x, y) 的电位就是：
[ V(x,y) = \frac{1}{4\pi\varepsilon_0} \left( \frac{q}{\sqrt{(x-d/2)^2 + y^2}} - \frac{q}{\sqrt{(x+d/2)^2 + y^2}} \right) ]
1.2 为什么要用MATLAB做场模拟
市面上能画电磁场的软件不少，COMSOL、Ansys这些专业有限元工具功能强大，但它们的定位是工程计算，建模流程繁琐，而且初学者很难在半小时内把结果跑出来。MATLAB做这件事的优势在于“脚本即实验”：你写几行代码，马上能看到结果，改参数重新运行，思路连续不断。
另一个关键点是，MATLAB的矩阵运算天然适配网格化场计算。电场模拟的本质就是在空间离散点上计算物理量，这正好是矩阵的高维运算。你不需要写for循环一个个点去算，用向量化操作一次就能算出整个平面的数据，速度快，代码也简洁。
还有一点是MATLAB的可视化工具相当成熟。contour、quiver、surf这几个函数分别对应等势线、矢量场、三维曲面，配合起来可以一图看全电偶极子的所有信息。这是做物理演示时特别实用的功能。
1.3 模拟方案的选型思路
我见过不少人一上来就搞三维体积渲染，把空间划分成三维网格，然后每个点算电位。这个方案视觉效果好，但计算量大，而且三维图里等势面不好观察。我更推荐的方案是先做二维平面模拟，把x-y平面上的电位分布画出来，再叠加电场矢量箭头。这个方案的好处有三个：
第一，二维图信息完整。电偶极子具有轴对称性，在包含偶极矩方向的平面里已经包含了所有关键信息；另一个垂直平面的情况几乎相同，只是旋转了一下。
第二，contour和quiver两种图可以叠加在同一坐标轴上，等势线和电场线的正交关系一眼就能看出来，这是理解电磁场的重要直观抓手。
第三，计算量小，普通笔记本完全没压力。三维图可以作为补充展示，但没必要一开始就上。
2. 核心代码实现：从参数设置到电位计算
2.1 参数设置与网格生成
先看完整的初始化和网格生成代码：
% 电偶极子参数 q = 1e-9; % 电荷量，单位C，这里取1nC d = 0.1; % 电荷间距，单位m epsilon0 = 8.854187817e-12; % 真空介电常数 % 正负电荷位置 pos_charge = [d/2, 0]; neg_charge = [-d/2, 0]; % 网格范围与密度 x_range = [-0.5, 0.5]; y_range = [-0.5, 0.5]; N = 200; % 网格点数，点数越多图越细腻 % 生成二维网格 x = linspace(x_range(1), x_range(2), N); y = linspace(y_range(1), y_range(2), N); [X, Y] = meshgrid(x, y);这里有两个参数很容易调错，我提醒一下。
第一个是网格点数N。点数太少，等势线会又粗又毛糙；点数太多，计算时间会明显上升，而且生成的图片文件也大。200左右是经验值，够细腻且运行快。如果你要出印刷级别的图，可以调到400或500，但普通演示200就行。
第二个是电荷间距d和网格范围的比例。d如果取得过大，比如0.2，而网格范围只有 -0.5 到 0.5，那正负电荷之间会出现剧烈的电位变化，画图时颜色映射会比较难看，等势线间距也不均匀。我的经验是让电荷间距占整个网格范围的五分之一到十分之一，这样既能看清近场细节，又能看到远场的衰减趋势。在这个例子里，d=0.1，网格范围是 [-0.5, 0.5]，比例是1:10，效果正好。
2.2 电位计算的精确公式实现
有了网格和电荷位置，直接算每个格点到两个电荷的距离，再叠加电位：
% 计算每个网格点到正负电荷的距离 r_plus = sqrt((X - pos_charge(1)).^2 + (Y - pos_charge(2)).^2); r_minus = sqrt((X - neg_charge(1)).^2 + (Y - neg_charge(2)).^2); % 防止距离为零导致除零错误 r_plus(r_plus < 1e-6) = 1e-6; r_minus(r_minus < 1e-6) = 1e-6; % 电位叠加 V = (1 / (4 * pi * epsilon0)) * (q ./ r_plus - q ./ r_minus);这一小段是整个模拟的核心，原理就是静电场的叠加定理：空间中任意一点的电位等于各电荷单独存在时产生电位的代数和。电位是标量，所以直接加减就行。
q 取 1e-9（1纳库仑），d取0.1米，算出来V的量级在几伏到几百伏之间，画图时需要处理一下数值范围，不然大部分区域的颜色会非常接近，细节看不清。常用的办法是只画有限范围的等势线，或者对电位做符号处理后再画。后面可视化部分我会专门说。
关于除零错误的处理：虽然理论上电荷只在两个孤立的点上有值，但网格点有可能恰好落在电荷所在坐标上。一旦出现距离为0，除以0会得到Inf，后面整个图都会出现异常。用r_plus(r_plus < 1e-6) = 1e-6这句话把过于小的距离替换掉，是常规做法。注意阈值不能设太大，否则会扭曲电荷附近的真实电位。
2.3 电场计算的两种方式
画电场线需要知道每个点的电场矢量E。有两种方式可以获得，我建议都用一下，各有用处。
第一种方式是从电位梯度计算，( \mathbf{E} = -\nabla V )。在数值上可以用gradient函数实现：
% 用电位梯度计算电场 [Ex, Ey] = gradient(-V, x(2) - x(1), y(2) - y(1));注意这里gradient的第一个参数是 -V，因为电场是电位的负梯度。x(2)-x(1) 是x方向的步长，y(2)-y(1) 是y方向的步长，这两个参数不能省，不然梯度计算的尺度不对，画出来的箭头方向虽然正确但长度会失真。
第二种方式是从叠加定理直接计算：
% 直接叠加两个点电荷的电场 Ex = zeros(size(X)); Ey = zeros(size(Y)); % 正电荷贡献 Ex = Ex + (1/(4*pi*epsilon0)) * q * (X - pos_charge(1)) ./ r_plus.^3; Ey = Ey + (1/(4*pi*epsilon0)) * q * (Y - pos_charge(2)) ./ r_plus.^3; % 负电荷贡献（方向指向负电荷，所以是减法） Ex = Ex - (1/(4*pi*epsilon0)) * q * (X - neg_charge(1)) ./ r_minus.^3; Ey = Ey - (1/(4*pi*epsilon0)) * q * (Y - neg_charge(2)) ./ r_minus.^3;这里用到的公式是点电荷电场公式，正电荷产生的电场方向从电荷指向场点，负电荷产生的电场方向从场点指向电荷，所以代码里有一个减号。
实测下来，两种方式的远场结果几乎一致，近场会有一点点差别，因为梯度法在电位变化特别剧烈的地方会有数值误差。如果只是画矢量场看方向，直接用梯度法快且省事；如果要做定量分析，建议用叠加法。
2.4 电场强度幅值与归一化
画quiver图的时候有个视觉陷阱：如果直接用真实的电场幅值来画箭头，靠近电荷的地方箭头会特别长，远离电荷的地方箭头几乎看不见，整张图的信息会严重失衡。解决办法是归一化处理，只保留方向信息，统一箭头长度：
% 计算电场强度幅值 E_mag = sqrt(Ex.^2 + Ey.^2); % 归一化，用于绘制箭头方向 Ex_norm = Ex ./ E_mag; Ey_norm = Ey ./ E_mag;归一化之后，每个箭头的长度都一样，视觉上看的是方向分布。要显示幅值大小，我会用背景色或者箭头颜色来表示，这个后面会讲到。
补充一句：E_mag这个变量不要扔掉，后面画三维曲面、展示电场强度分布的时候还要用。
3. 可视化方案：一图看懂电偶极子的全部信息
3.1 等势线（contour）绘制技巧
画等势线最直接的方式是contour(X, Y, V, levels)，其中levels是一个向量，指定要画的等势线数值。但这里有个坑：V的范围跨度很大，正电荷附近能到几千伏，远处只剩下几伏，直接用线性刻度画，远处的等势线会全部挤在一起。
我的解决方法是把电位用sigmoid函数压缩一下：
% 对电位做非线性映射，压缩动态范围 V_scaled = 2 ./ (1 + exp(-V / 50)) - 1;这个公式把V映射到 (-1, 1) 区间，V=0 处映射为0。这样做的好处是：近处的大电位不会把颜色标尺撑爆，远处的细微变化也能显示出来。这里的50是缩放因子，具体值取决于你算出的V的量级。我试下来，V_max大约几百伏的情况下，50效果不错；如果V_max是几万伏，就要把缩放因子调大到500左右。你可以先跑一下看V的范围，再调这个参数。
然后画出等势线：
figure('Color', 'w'); % 等势线 [C, h] = contour(X, Y, V_scaled, 30); % 30条等势线 clabel(C, h, 'FontSize', 8, 'LabelSpacing', 150); xlabel('x (m)'); ylabel('y (m)'); title('电偶极子等势线分布'); axis equal; grid on;clabel函数会在等势线上标注数值，LabelSpacing控制标注之间的距离，数值越大标注越稀疏，避免标签重叠。这里有个很容易忽略的点：axis equal 必须加，不然坐标轴的刻度比例不一致，圆形的等势线会被拉成椭圆形，物理图像就错了。
3.2 电场矢量场（quiver）绘制
画电场方向用quiver。为了图面清晰，我不会在每一个网格点上都画箭头，那样200×200个箭头会把整个图涂黑。通常每隔几个点采样一次：
% 每隔4个点取一个箭头，避免箭头过密 skip = 4; figure('Color', 'w'); quiver(X(1:skip:end, 1:skip:end), ... Y(1:skip:end, 1:skip:end), ... Ex_norm(1:skip:end, 1:skip:end), ... Ey_norm(1:skip:end, 1:skip:end), ... 0.5); % 0.5是箭头缩放因子 hold on;这个0.5的缩放因子很关键。quiver默认会按数据自动缩放箭头长度，但归一化之后箭头的“长度”信息已经没意义了，你再让它自动缩放就会得到一堆莫名其妙的短箭头或者长箭头。手动给定缩放因子之后，箭头长度基本一致，只表达方向。你可以根据实际效果调0.3到0.8之间。
另外注意quiver画出来的箭头只是方向，不是电场线的连续轨迹。要得到真正的电场线（从正电荷出发到负电荷结束的连续曲线），得用streamline函数或者自己写流线追踪算法。这个我放在扩展内容里讲，基础版先用quiver就够了。
3.3 三维电位曲面：直观感受“势阱”与“势垒”
除了二维图，三维曲面图也是理解电位分布的好帮手。surf可以把电位值映射成曲面高度和颜色：
figure('Color', 'w'); surf(X, Y, V, 'EdgeColor', 'none'); colormap(jet); colorbar; xlabel('x (m)'); ylabel('y (m)'); zlabel('电位 V (V)'); title('电偶极子电位三维曲面'); view(35, 45);从这张图能明显看到两个“尖峰”：正电荷处电位冲向正无穷（图中显示为红色高峰），负电荷处电位跌向负无穷（蓝色深谷），远处电位逐渐趋于零。中间电中性的地方电位为零，形成一个“鞍点”，这个鞍点的位置在原点，是很多教材里讲的重点，三维图上看得特别清晰。
surf默认会画出网格线，数据点一多会显得杂乱，所以加上'EdgeColor', 'none'关掉网格线，只保留颜色映射，图面会干净很多。
3.4 合并画图：等势线+电场矢量一图流
把等势线和电场矢量画在同一张图上，是最推荐的展示方式：
figure('Color', 'w'); [C, h] = contour(X, Y, V_scaled, 30); clabel(C, h, 'FontSize', 8, 'LabelSpacing', 150); hold on; quiver(X(1:skip:end, 1:skip:end), ... Y(1:skip:end, 1:skip:end), ... Ex_norm(1:skip:end, 1:skip:end), ... Ey_norm(1:skip:end, 1:skip:end), ... 0.5, 'Color', [0.6, 0.2, 0.2]); hold off; xlabel('x (m)'); ylabel('y (m)'); title('电偶极子等势线与电场矢量分布'); axis equal; grid on;这张图能非常直观地验证两个基本定理：
第一，电场矢量和等势线处处正交。你会发现箭头的方向总是垂直于等势线的切线方向，这正是电场与等势面的几何关系。
第二，电场从正电荷出发、终止于负电荷。箭头的方向在正电荷附近是向外辐射的，在负电荷附近是指向电荷的，整张图呈现出从正到负的“流线”感。
我把quiver箭头的颜色设成暗红色，和黑色的等势线区分开，图面层次更分明。如果你喜欢其他配色，也可以根据实际需要调整RGB值。
3.5 电荷位置标记
图上最好把正负电荷的位置标出来，直观对应：
plot(pos_charge(1), pos_charge(2), 'r+', 'MarkerSize', 12, 'LineWidth', 2); plot(neg_charge(1), neg_charge(2), 'bo', 'MarkerSize', 12, 'LineWidth', 2);这里用红色加号表示正电荷，蓝色圆圈表示负电荷。如果你的图上已经有quiver的红色箭头了，可以给正电荷换一个颜色，比如用绿色五角星，避免同色混淆。根据自己的审美调整就行。
4. 常见问题与调试技巧实录
4.1 等势线图“糊成一团”怎么办
这是出现频率最高的问题。原因几乎都是V的动态范围太大，电荷附近电位极高，远处的电位在颜色映射中被压到同一个色阶。解决办法就是我在3.1节说的非线性映射V_scaled = 2 ./ (1 + exp(-V / 50)) - 1。这个公式处理之后的等势线在近处和远处都能清晰区分。
如果你觉得 sigmoid 映射还是不够好，也可以试另一个思路：只在某个有限的电位范围内画等势线。比如：
V_min = -100; V_max = 100; levels = linspace(V_min, V_max, 20); contour(X, Y, V, levels);这种方式会把超出范围的区域留白，等于只截取中间电位区域来观察。对于想看近场细节的场景，这个方案更直接。
4.2 quiver箭头方向“乱跳”是什么原因
箭头方向看起来乱，多数情况下是因为箭头采样过密或者缩放因子设置不当。skip参数取4、5、6都可以尝试，要保证箭头之间有足够间距。如果箭头还是看起来乱，可能是纯数值假象：在E_mag极小的地方（比如电位鞍点附近），归一化之后的Ex_norm、Ey_norm方向会发生突变，因为真实的电场幅值趋近于零，方向本身就不稳定，再除以很小的E_mag就会放大误差。
处理方法是给E_mag设一个下限，低于下限的点不画箭头：
valid = E_mag > 1e-3; % 阈值根据实际数据调整 Ex_norm(~valid) = NaN; Ey_norm(~valid) = NaN;把无效点设为NaN，quiver就不会在这些位置画箭头了。这个技巧在等势线密集的区域特别有用，能让图面干净很多。
4.3 电位数值太大或太小，colorbar显示不正常
如果电位量级是几千伏或者微伏，前面说的V/50这个缩放因子就需要调整。你可以先跑一行代码查看范围：
fprintf('V min = %.3e, V max = %.3e\n', min(V(:)), max(V(:)));然后根据量级设置缩放因子。我一般是让缩放因子等于V_max的十分之一左右，然后再微调。这个值的大小直接影响等势线的疏密分布，值得花一点时间调到视觉效果最佳。
4.4 跑得慢怎么办
如果N取得很大（比如500以上），计算量会明显增加。其实没必要优化代码，最简单的方法就是先把N设小一点（150左右）快速调试，确认代码没问题、图形符合预期后，再调大N出最终图。开发流程上这叫“先跑通再跑好”，不要一开始就用高精度参数。
如果需要经常调整参数，你可以把代码封装成函数：
function plot_dipole(q, d, N, x_range, y_range) % 画电偶极子的等势线、电场矢量、三维曲面 ... end后续只要改参数调用函数就行，不用每次都对着脚本改来改去。
4.5 单位制混乱导致结果不对
这个问题一般不报错，但对错完全看不出来。我见过有人把电荷间距写成厘米，电荷量写成微库仑，导致最后画的电场分布比例失衡。建议一开始就用国际单位制：电荷用库仑（C），长度用米（m），介电常数用真空介电常数8.854187817e-12 F/m。全都统一，不要混用。
如果你做的是理论示意而非真实物理量计算，也可以把常数简化。比如令1/(4πε₀) = 1，只看场分布的形状。这样数值上更清爽，画出来的图形形状完全一样。但简化会失去物理量级的参考，用于课程展示时最好保留真实值，方便算具体电压。
5. 深度扩展：从静态图到真正可用的研究工具
5.1 用streamline画真正的电场线
quiver展示的是离散的矢量场信息，而电场线是连续的曲线。MATLAB的streamline函数可以根据矢量场数据追踪流线，效果类似于手绘的电场线：
% 设置起始点（围绕正电荷一圈） theta = linspace(0, 2*pi, 16); start_x = pos_charge(1) + 0.02 * cos(theta); start_y = 0.02 * sin(theta); figure('Color', 'w'); hold on; streamline(X, Y, Ex, Ey, start_x, start_y); contour(X, Y, V_scaled, 20, 'LineWidth', 0.5); hold off;这里的start_x和start_y是流线的起点坐标，我取的是以正电荷为中心、半径0.02米的圆周上的16个点。streamline会从这些点出发，沿着矢量场的方向追踪电场线。由于电场线从正电荷出发终止于负电荷，所以这些流线最终都会收束到负电荷附近。
这个效果比quiver更有“物理味道”，推荐你要做演示时用。
5.2 电偶极子的动态旋转演示
静态图看多了，可以把偶极子旋转起来，观察空间场分布随方向的变化。只需要在函数外加一层角度循环，把电荷坐标用旋转矩阵更新一下：
for angle = 0:5:180 theta_rad = deg2rad(angle); rot = [cos(theta_rad), -sin(theta_rad); sin(theta_rad), cos(theta_rad)]; pos_charge = rot * [d/2; 0]; neg_charge = rot * [-d/2; 0]; % 重算电位和电场，重新绘图 % ... drawnow; end这样就能生成一个动态旋转的偶极子场分布。如果保存成gif，还可以直接放到PPT里当课件动画。
5.3 偶极子阵列的应用
实际工程问题里，单个偶极子很少，更多是偶极子阵列。比如天线阵列中，多个偶极子叠加形成特定的辐射方向图。MATLAB里实现阵列叠加热别简单，只要把每个偶极子的电位相加就行：
function V_total = dipole_array_potential(X, Y, charges) % charges 是 Nx3 矩阵，每行表示一个偶极子的 [x, y, theta] % 其中 x, y 是偶极子中心位置，theta 是偶极子方向角 V_total = zeros(size(X)); for i = 1:size(charges, 1) cx = charges(i, 1); cy = charges(i, 2); theta = charges(i, 3); % 计算局部坐标 X_local = (X - cx) * cos(theta) + (Y - cy) * sin(theta); Y_local = -(X - cx) * sin(theta) + (Y - cy) * cos(theta); % 叠加该偶极子的电位 V_total = V_total + single_dipole_potential(X_local, Y_local); end end这种叠加逻辑对理解天线阵列的工作原理很有帮助。
5.4 加上介质材料后的变化
如果在偶极子附近放置一块介质材料，电位分布会改变，原因是介质在电场作用下发生极化，在表面产生束缚电荷。严格计算需要求解拉普拉斯方程，并考虑不同区域的介电常数差异。MATLAB可以用PDE工具箱处理这类问题，但代码复杂度明显上升，不是一两个脚本能解决的。
如果你只是想在课上展示介质的影响，可以用近似模拟方式：在介质区域人为添加一些等效的极化电荷，观察等势线因极化电荷扭曲的效果。这个方案物理上不是特别精确，但视觉上很有说服力。做科研时还是老老实实上专业工具吧。
6. 好习惯与收尾经验
最后分享几个我在写这类仿真代码时的习惯，能省下不少调试时间。
第一个习惯是尽量变量化所有参数。q、d、N、范围、缩放因子全部定义在文件最前面，不要算到一半发现想改电荷间距，还得去代码中间找。这看起来像废话，但实际操作中很多人偷懒，写着写着又回到硬编码的老路上。
第二个习惯是保留V和E的原始数据，不要只用画图处理之后的数据。有些时候你需要回到原始数据做定量分析，比如计算某个点的精确电场值，或者沿着某条线提取电位分布。如果原始数据被覆盖了，就得重跑一遍。
第三个习惯是注释要写“为什么”，不要写“是什么”。比如r_plus(r_plus < 1e-6) = 1e-6旁边，注释应该写“防止网格点落在电荷上导致除零”，而不是写“赋值1e-6”。做代码重构或者回头看旧代码的时候，这类注释能救你命。
第四个习惯是把成品代码保存成函数，加上默认参数。我的做法是写一个demo_dipole()函数，内部调用plot_dipole(q, d, N, range)，这样既能一键运行看效果，又能灵活改参数做实验。
这些习惯一开始坚持起来会觉得麻烦，但坚持一个学期之后，你会发现自己写代码的速度和调试效率都上来了。回到电偶极子本身，真正把这个模型跑通、画清楚之后，我对电磁场的理解确实上了一个台阶。很多东西光靠公式推导确实抽象，但一旦把等势线和电场矢量画在屏幕上，脑子里就建立起图像了，后面处理更复杂的电磁问题也会更有底气。希望这套代码和思路对你也有帮助。
如果你跑的过程中遇到报错或者图形不符合预期，欢迎评论交流，我看见了会回复。
