自己画过电偶极子的电场线吗我指的不是课本上那种标注好的示意图而是真的打开MATLAB让屏幕上几千个箭头、几十条流线自己排布出来的那种。我第一次跑通这个模拟的时候最直观的感受是书里那页图终于变成了我自己的东西。电偶极子看起来是电磁场理论里最基础的教学模型但把它做成MATLAB场模拟一次性能打通“物理公式—数值计算—可视化表达”整条链路特别适合正在学电磁场的大学生、要给课件配动态图的老师以及想在MATLAB里练手矢量场可视化的人。这篇文章我会直接给出一套可运行的源码并拆解每一步为什么这么写。你拿过去哪怕之前没怎么碰过MATLAB跟着敲一遍也能跑出等势面、电场矢量、电场线和三维电位曲面。后面我还会聊几个真踩过的坑以及怎么在这个模型上做动画和扩展帮你把这一个小demo玩出花来。1. 为什么我建议第一个场模拟就做电偶极子1.1 公式简单扩展性却不简单电偶极子说白了就是两个相距很近的等量异号点电荷一个带正电一个带负电。在电磁场里它是暴露度极高的一个简化模型天线的近场可以看成振荡偶极子介质极化后可以用偶极子密度来描述分子物理里很多电学性质也拿它做起点。关键的是它的解析表达式足够简单。空间任意一点的电位就是两个点电荷电位直接叠加V(x,y) k * q / R1 - k * q / R2其中 R1、R2 分别是从正负电荷到场点的距离。电场也不用绕弯子求解什么方程直接用库仑定律的矢量形式做叠加就可以了。这种“公式足够简单、物理图像足够丰富”的模型最适合拿来当第一个数值模拟项目。我见过不少同学一上来就抱着有限元工具箱去算场结果被网格剖分、边界条件搞得头大最后连场长什么样都没看清。不如先从电偶极子这种闭式解入手把核心逻辑跑明白再往复杂模型走。1.2 一次打通三层能力做这个项目实际上是同时练了三件事。第一把连续公式变成离散网格。书本上的场是连续的任意一点的电位都能算但计算机没法处理无穷多个点只能在有限的网格上采样。这一步要绕明白 meshgrid、索引、矩阵运算是很多数值模拟的公共基础课。第二把标量场和矢量场分开理解和表达。电位是标量只有大小没有方向电场是矢量既有大小又有方向。两者通过梯度关系关联但可视化手段完全不同。标量场适合用色块、等高线、三维曲面表达矢量场适合用箭头、流线表达。这个区分想清楚了到流体模拟、热场分析里也通用。第三学会迭代调优图形效果。算对和画好看是两回事。颜色映射的范围怎么定、箭头取多密、三维曲面怎么看才不骗人这些没有哪个教科书会系统教你只能靠试。电偶极子模型小、计算快参数怎么改都不会崩简直是练可视化的天然沙盒。2. 从解析公式到网格先把场算对2.1 电位和电场的叠加式怎么落到代码我习惯把整个模拟分成三段参数定义、物理量计算、可视化。第一段定义电荷量、间距、网格范围第二段算 V、Ex、Ey第三段才去画图。核心计算部分代码可以写成这样clear; clc; % 物理参数 k 8.99e9; % 库仑常数N·m^2/C^2 q 1e-9; % 电荷量C d 0.1; % 偶极子间距m % 正负电荷的位置 r1 [ d/2, 0]; % q r2 [-d/2, 0]; % -q % 网格定义 x -1:0.02:1; % x轴采样点 y -1:0.02:1; % y轴采样点 [X, Y] meshgrid(x, y); % 场量计算 R1 sqrt((X - r1(1)).^2 (Y - r1(2)).^2 1e-10); R2 sqrt((X - r2(1)).^2 (Y - r2(2)).^2 1e-10); % 电位叠加 V k * q ./ R1 - k * q ./ R2; % 电场叠加库仑定律矢量形式 Ex k * q .* (X - r1(1)) ./ R1.^3 - k * q .* (X - r2(1)) ./ R2.^3; Ey k * q .* (Y - r1(2)) ./ R1.^3 - k * q .* (Y - r2(2)) ./ R2.^3; % 电场模长 E sqrt(Ex.^2 Ey.^2);这里很多人会问为什么算距离的时候在根号里面加一个1e-10其实就是为了防止场点和电荷位置完全重合时出现除零。网格步长取 0.02采样范围是 -1 到 1得到的矩阵是 101×101总共一万多个点。这个规模对普通电脑来说毫无压力计算部分是瞬间完成的。你可以试试把步长改成 0.005分辨率更高但计算量会变成原来的16倍画图时也会明显变慢。2.2 meshgrid 到底帮你做了什么meshgrid这个函数新手第一次用往往会懵。它的作用是把一维的 x 向量和 y 向量扩展成二维的坐标矩阵[X, Y] meshgrid(-1:1:1, -2:1:2)返回的 X 和 Y 是两个同尺寸的矩阵X 里的每一行都等于原始的 x 向量Y 里的每一列都等于原始的 y 向量。于是 X(i,j)、Y(i,j) 这两个矩阵元素合在一起就唯一确定了网格上第 i 行第 j 列那个点的坐标。大多数人踩的第一个坑是不理解为什么 X 和 Y 的尺寸是 length(y) 行、length(x) 列。记住一句话行对应的是 y 方向列对应的是 x 方向。所以在画图时contourf(X, Y, V)里三个矩阵尺寸必须完全一致而xlabel、ylabel的方向和自然坐标系一致不会转置。如果你不太习惯这种行和列的概念也可以直接把 x、y 两个向量传给网格函数。MATLAB 的很多绘图函数会自动调用 meshgrid但我还是建议你显式写出来因为后面算 R1、R2、Ex、Ey 都要用矩阵逐点操作。2.3 点乘和矩阵乘法别搞混这段代码里密集使用.*、./、.^它们的意思是逐元素运算。也就是说两个矩阵相同位置的元素逐一计算结果矩阵的形状不变。如果没有那个点MATLAB 会尝试做矩阵乘法或矩阵除法那完全就是另一回事了轻则报错重则算出离谱的结果。我见过最经典的低级错误是Ex k * q * (X - r1(1)) / R1.^3;写成了矩阵除法。表面看似乎没问题但 MATLAB 会把它按线性方程组求解处理输出一个标量或者一个意料之外的矩阵。等你画图的时候要么维度对不上要么图形彻底乱掉。所以只要涉及场点坐标和距离矩阵的运算一律无脑用点运算符。3. 可视化三步走等势面、矢量场和三维曲面3.1 用 contourf 把电位铺成云图算完 V 之后最简单的展示方式就是contourf填充等势面。它会按电位大小把整个平面染成不同颜色直接用色块告诉你哪里电位高、哪里电位低。figure(Color, w); contourf(X, Y, V, 40, LineStyle, none); colorbar; axis equal; xlabel(x / m); ylabel(y / m); title(电偶极子电位等值面云图);40表示画 40 层等势间隔。数字越大颜色过渡越细腻。不过这里有个细节因为点电荷附近电位趋近于正负无穷直接按线性间隔画等势线会非常密集地压在电荷附近而远处区域几乎看不出梯度差异。我的处理方法是把颜色轴范围限制住只显示我们有兴趣的那一段。MATLAB 新版本用clim老版本用caxisvmax 20; % 只显示 -20V 到 20V 之间的电位 clim([-vmax vmax]);这样图面会干净很多正负电位的分布对比也清楚。不信你跑一下对比限制前后的颜色条马上就明白这个经验的价值。3.2 quiver 画电场方向归一化是基本操作电位图只能看到大小方向还得靠矢量箭头。quiver是 MATLAB 里画二维箭头场的标准函数figure(Color, w); contourf(X, Y, V, 40, LineStyle, none); colorbar; hold on; % 每隔 5 个点取一个箭头太密会糊成一团 step 5; quiver(X(1:step:end, 1:step:end), ... Y(1:step:end, 1:step:end), ... Ex(1:step:end, 1:step:end), ... Ey(1:step:end, 1:step:end), ... 1.2, k, LineWidth, 1); plot(r1(1), r1(2), ro, MarkerSize, 12, MarkerFaceColor, r); plot(r2(1), r2(2), bo, MarkerSize, 12, MarkerFaceColor, b); axis equal; title(电位云图叠加电场矢量);这里我要重点说一个坑如果你直接用quiver(X, Y, Ex, Ey)而不做处理得到的结果会非常糟糕。因为电场强度随距离衰减是平方反比的电荷附近的箭头大得吓人边界处的箭头又小到看不见整张图的动态范围极大中间大部分区域全是密密麻麻的小针完全没法看。我一般把矢量拆成“方向”和“大小”两层方向用归一化箭头展示大小用云图底色展示。方法很简单每个方向分量除以模长En max(E, 1e-12); quiver(X(1:step:end, 1:step:end), Y(1:step:end, 1:step:end), ... Ex(1:step:end, 1:step:end) ./ En(1:step:end, 1:step:end), ... Ey(1:step:end, 1:step:end) ./ En(1:step:end, 1:step:end), ... 0.6, k);这样所有箭头长度保持一致图上表达的是该点的场方向而场强大小由背景色表达。这是论文级可视化里非常常见的做法自己写论文配图时也用得上。3.3 streamslice 画电场线矢量箭头表达的是离散方向采样而电场线这种连续的流线用streamslice更合适。它是专门为二维矢量场设计的一条线绘制函数自动从网格边界或者指定位置发线figure(Color, w); contourf(X, Y, V, 30, LineStyle, -, LineWidth, 0.5); hold on; h streamslice(X, Y, Ex, Ey, 2); set(h, Color, [0.2 0.2 0.2], LineWidth, 1.2); plot(r1(1), r1(2), ro, MarkerSize, 12, MarkerFaceColor, r); plot(r2(1), r2(2), bo, MarkerSize, 12, MarkerFaceColor, b); axis equal tight; title(电偶极子电场线streamslice);streamslice的第二个参数是一个密度系数我一般取 1 到 3 之间。取值越大线条越密。它会自动绕开电位变化不连续的奇异区域线条非常平滑。有一点必须提醒如果你的网格范围取得太小或者电场分量里有 NaN、Infstreamslice可能什么都不画出来也不报错。遇到这种情况先用sum(isnan(Ex(:)))查一下数据多半是距离矩阵在有电荷位置出现了除零。3.4 三维电位曲面二维图看平面分布足够了但想要更直观的“高度感”可以用surf画三维电位曲面。严格来说这个面不是物理上的某个面只是把电位数值当高度来展示figure(Color, w); surf(X, Y, V, EdgeColor, none); colormap(jet); colorbar; lighting gouraud; light(Position, [0 0 1]); xlabel(x / m); ylabel(y / m); zlabel(电位 / V); title(电偶极子三维电位曲面);正电荷的位置拔起一座尖峰负电荷的位置陷下一个深坑中间鞍部就是电位为零的分界线。这个三维图拿去给不懂电磁场的朋友看也能一眼理解偶极子“一头高一头低”的电位分布。不过三维图有个隐含陷阱观察角度会骗人。MATLAB 默认视角下曲面高度很容易被理解为“电场强度”但实际它表示的是电位。所以如果你要放在报告里一定在图注里写清楚纵轴是电位 V而不是电场 E。这一点我在第 5 章还会展开说。4. 让场动起来参数扫描与动画输出4.1 循环体里改变偶极子间距静态图看多了你会发现最有冲击力的其实是动态演示。让偶极子间距周期性地变大变小整个场结构跟着呼吸一样地变化课上展示效果相当好。基本思路就是在一个循环里重新计算场量然后更新图形x -1:0.02:1; y -1:0.02:1; [X, Y] meshgrid(x, y); k 8.99e9; q 1e-9; figure(Color, w); for t 1:150 % 间距在 0.02 到 0.28 之间变化 d 0.02 0.26 * abs(sin(t / 20)); r1 [ d/2, 0]; r2 [-d/2, 0]; R1 sqrt((X - r1(1)).^2 (Y - r1(2)).^2 1e-10); R2 sqrt((X - r2(1)).^2 (Y - r2(2)).^2 1e-10); V k * q ./ R1 - k * q ./ R2; Ex k * q .* (X - r1(1)) ./ R1.^3 - k * q .* (X - r2(1)) ./ R2.^3; Ey k * q .* (Y - r1(2)) ./ R1.^3 - k * q .* (Y - r2(2)) ./ R2.^3; clf; contourf(X, Y, V, 40, LineStyle, none); hold on; plot(r1(1), r1(2), ro, MarkerSize, 12, MarkerFaceColor, r); plot(r2(1), r2(2), bo, MarkerSize, 12, MarkerFaceColor, b); axis equal tight; colorbar; clim([-20 20]); title(sprintf(电偶极子动态模拟 d %.3f m, d)); drawnow; end这段代码运行时正负电荷的位置不断分开、靠拢等势面的形状也实时跟着变化效果非常直观。注意循环里每次都用clf清空当前图形然后重新画虽然简单但性能一般。性能优化我在 4.2 里单独讲。4.2 帧率控制和性能优化drawnow的作用是强制刷新图形窗口没有它MATLAB 会攒着所有绘图操作等循环结束了一次性画你就看不到动画过程了。我习惯在drawnow后面加一句pause(0.05)一方面是让动画不会闪得太快另一方面是给图形系统一点缓冲时间避免 UI 卡死。帧率不用太激进15 到 20 帧每秒就足够流畅了。如果追求更高性能尽量不要用clf清空整个图窗再重绘。更好的办法是创建一个图形对象然后在循环里更新它的XData、YData、CData等属性。这种做法在写 GUI、做实时数据采集显示时会频繁用到。figure(Color, w); contourf(X, Y, V, 40, LineStyle, none); hold on; hPlot plot(r1(1), r1(2), ro, MarkerSize, 12, MarkerFaceColor, r); plot(r2(1), r2(2), bo, MarkerSize, 12, MarkerFaceColor, b); colorbar; for t 1:150 d 0.02 0.26 * abs(sin(t / 20)); % 重新计算 R1, R2, V, Ex, Ey ... % 更新等势图颜色数据 contourf(X, Y, V, 40, LineStyle, none); % 更新电荷位置 set(hPlot, XData, d/2, YData, 0); title(sprintf(d %.3f m, d)); drawnow; pause(0.05); end不过说实话对电偶极子这个量级的模拟性能差别感受不明显。我通常只在图形对象数量非常多、需要实时响应时才认真做对象化重构。4.3 用 VideoWriter 把动画存成文件动态效果只在屏幕上播放还不够课上展示、实验报告、小组汇报都需要视频文件。MATLAB 的VideoWriter可以非常方便地逐帧保存v VideoWriter(dipole_demo.avi); v.FrameRate 20; open(v); % 在动画循环内部每帧绘制完成后加两行 frame getframe(gcf); writeVideo(v, frame); % 循环结束后关闭 close(v);FrameRate设置 20 帧每秒比较合适。如果想让视频节奏更慢可以改成 10 到 15。getframe会把当前图形内容截成一张图像逐帧写入后就是一个可直接播放的 AVI 文件。想导出 GIF 的话MATLAB 没有完全原生的一键功能。可以借助exportgraphics配合循环逐帧写或者干脆生成 PNG 序列再用第三方工具合成。命令行下很多小工具都能完成这里就不展开了。5. 实测踩坑五个反复出现的隐性错误5.1 距离矩阵里的“幽灵”零值这个坑我在第一次自己写的时候也踩过。计算 R1、R2 时如果网格点恰好落在电荷所在的坐标上距离就是 0分母除零导致 V、Ex、Ey 出现 Inf 或 NaN。很多人在公式上反复检查也没发现问题最后把网格一放大发现某几个点数值乱七八糟。解决办法就是我在代码里写的 1e-10。这个微小偏移只影响电荷附近极其局部的位置对整个场的分布几乎不产生可感知的误差但可以彻底避免除零问题。如果你在算更复杂的电场比如多个电荷叠加也一定要在每一个距离项后面都加上这个小量。5.2 quiver 箭头全部一样大这个坑相当隐蔽。默认情况下quiver会按矢量长度自动缩放箭头本来是个好事。但当你对Ex、Ey做了归一化之后所有箭头长度都会变成 1此时如果还保留默认缩放参数很可能所有箭头画出来一样大。如果你想让某些区域强调方向、某些区域强调大小那就别做归一化而是用原始的Ex、Ey加限制长度的方式quiver(X, Y, Ex, Ey, 0.8);0.8是自动缩放系数小于 1 会缩小箭头。多试几个值找到一个观感最好的缩放比例。如果图上同时需要等势面和箭头我建议等势面表达大小归一化箭头表达方向这样各司其职。5.3 三维纵轴误导电位面高度不等于电场强度三维电位曲面很容易被外行误读。有一次我用这个图去给课题组做组会展示讲完以后有同事问我“为什么电场强度在这个位置不是最大”其实是因为曲面 y 轴高度是电位 V电场强度是 V 的梯度也就是曲面最陡的方向和陡峭程度而不是高度本身。如果你想让三维图直接展示电场强度可以画surf(X, Y, E)而不是surf(X, Y, V)。但要注意 E 在电荷附近也会趋近无穷同样需要clim截断。搞清楚“你画的是哪个物理量”比任何调参技巧都重要。5.4 动态循环里 clf 和 cla 的区别clf清空整个 figurecla清空当前坐标轴。动态动画里如果用clf坐标轴、colorbar、标题等所有对象都会消失需要重新创建用cla则只清空曲线保留坐标轴属性。我自己测试下来contourf这类底层绘图对象更新时cla并不总是能完全释放旧对象的句柄偶尔会内存持续增长。如果你观察到循环越跑越慢多半就是这里的问题。最简单的兜底方案就是每帧clf在这个模型的计算量下完全没问题。5.5 脚本命名别跟内置函数撞车最后说一个很基础但特别容易忽略的如果你把脚本命名为dipole.m或者field.m有可能和 MATLAB 自带的工具箱函数冲突。我建议用dipole_sim.m、dipole_visualize.m这种带下划线的具体命名避免在命令行敲dipole时调出奇怪的东西。遇到“函数名与 MATLAB 内置函数重名”导致的诡异报错时可以先在命令行输入which 脚本名看看 MATLAB 到底解析到哪个路径的文件。这个排查思路在大型项目里同样有效。6. 从偶极子到更多扩展思路与函数化6.1 改成四极子只需加两个电荷偶极子的代码框架是通用的。想改成电四极子只需要把电荷数目从 2 个变成 4 个在对应的位置再叠加两次就行。比如线性四极子中间是 -2q两端是 q 和 q或者更直观的四角布局charges [ q, d/2, 0; q, -d/2, 0; -q, 0, d/2; -q, 0, -d/2 ];然后用一个循环叠加所有电荷的贡献。你会发现四极子的场形状比偶极子复杂得多等势面会出现四个象限的对称结构。这一步扩展做下来你对“叠加原理”的理解会彻底从公式层面落到图像层面。6.2 旋转偶极子的模拟边界想让偶极子绕原点旋转只需要在循环里让电荷坐标随时间变化theta t * 0.05; r1 [ d/2 * cos(theta), d/2 * sin(theta)]; r2 [-d/2 * cos(theta), -d/2 * sin(theta)];这样能模拟一个旋转偶极子的近场变化。但必须提醒一点这里用到的只是静态库仑场的叠加不是完整的辐射电磁波模拟。没有考虑推迟势和电磁波传播不能直接当成天线辐射方向图来看。如果想计算真正的偶极子辐射方向图需要在远场近似下引入矢势和角分布那是另一个更深入的话题。6.3 把主脚本封装成可复用函数当模拟跑顺了我建议你把算场的部分抽成函数方便后续调用function [X, Y, V, Ex, Ey, E] dipoleField(q, d, xRange, yRange, step) % 输入电荷量 q偶极子间距 d范围 xRange、yRange网格步长 step % 输出网格坐标 X、Y电位 V电场分量 Ex、Ey模长 E x xRange(1):step:xRange(2); y yRange(1):step:yRange(2); [X, Y] meshgrid(x, y); r1 [ d/2, 0]; r2 [-d/2, 0]; R1 sqrt((X - r1(1)).^2 (Y - r1(2)).^2 1e-10); R2 sqrt((X - r2(1)).^2 (Y - r2(2)).^2 1e-10); k 8.99e9; V k * q ./ R1 - k * q ./ R2; Ex k * q .* (X - r1(1)) ./ R1.^3 - k * q .* (X - r2(1)) ./ R2.^3; Ey k * q .* (Y - r1(2)) ./ R1.^3 - k * q .* (Y - r2(2)) ./ R2.^3; E sqrt(Ex.^2 Ey.^2); end封装完之后你的主脚本只需要传入不同参数就能快速出图[X, Y, V, Ex, Ey, E] dipoleField(1e-9, 0.1, [-1 1], [-1 1], 0.02);这样无论是做参数扫描、写作业还是给组会临时画个示意图都能一行命令搞定。我后来写很多场相关的报告都直接复用这个函数改改电荷排布省了大量重复代码。从电偶极子出发你已经拥有了一个可以扩展到任意点电荷组合的模拟框架。换个电荷位置、多叠几组电荷就能得到完全不同物理内容。这正是我特别推荐从这个小项目入手的原因它小到你能完全掌控又大到能带你走完场模拟的完整流程。