简介这份基于 Matlab 的仿真实验报告面向电子信息、电气工程等专业学生及电磁场初学者围绕线电荷周围电场与电位分布的建模与可视化展开可直接用于课程实验、仿真练习与实验报告写作参考。内容以库仑定律、电场强度与势函数为理论主线先讲清将线电荷视为点电荷集合的微积分思路再给出完整的 Matlab 实现设置 50 个线电荷、用 meshgrid 构建坐标网点、循环累加总电势并调用 surf、contour 与 quiver 分别呈现电势三维曲面、等位线和归一化场强矢量图末尾附实验结论以及线电荷密度选取、矩阵维度一致性等排错经验。成品为 1 个 docx 文档压缩包约 227KB涵盖实验目的、原理、步骤、代码与结果图表共 5 页结构完整可直接套用。已有 3041 人学习适合希望快速入门电磁场数值仿真流程或需要现成实验报告框架的读者参考。1. 图跑得挺漂亮数量级却是错的matlab 线电荷电场仿真的第一课第一次照着讲义把这段 matlab 代码敲进去surf 出来的三维电势面漂亮得让我当场截图——等位线齐整quiver 的绿箭头铺满整个网格看上去一切正常。直到某天顺手用两点电势差除以距离去反推场强得到的数跟图上箭头方向对得上、大小却差了几十倍才发现这套代码从常数写法到梯度算子至少有五处是「凑巧跑通」。这篇不写实验报告模板而是把线电荷电场仿真从公式到出图这段路上真正绊人的地方一条条拆开参数怎么定、三张图各管什么、nr 取多少算够、哪里会静默出错。刚上手 matlab 做电磁场仿真、既要交报告又想讲清原理的人照着改参数就能跑做过一遍但没验证过数值的人大概能在里面找到自己踩过的那两个坑。2. 从库仑定律到能跑的循环把物理公式翻译成矩阵运算2.1 先算电势再求梯度比直接合成场强省一半事点电荷的电场是矢量要算空间任一点的 E得把每个电荷产生的场按 x、y 分量分别投影再叠加代码里就是两套循环、两套符号判断。而电势是标量直接把 51 个点电荷的贡献加起来就行等算完了再用 E -∇U 反推场强一次 gradient 调用拿到 Ex、Ey 两个分量。这条路径是这类仿真里最常见的做法不是说它精度最高而是它把「矢量叠加」这件容易写错的事交给了一个 matlab 内置函数去兜底。代价也很明确gradient 是数值差分精度天然低于解析求导而且在网格边界只能做单侧差分边界那一圈算出来的场强方向不可信。所以看 quiver 图的时候我对最外两圈箭头基本不做判断只关心中间区域。想验证某一点的场强准不准正确做法是拿解析式 q/(4πε₀r²) 单独算一个点跟它比而不是盯着图看箭头长得顺不顺眼。电势叠加的写法本身没有歧义U Σ q/(4πε₀·rₖ)rₖ 是第 k 个电荷到目标点的距离。这里有个前提容易被忽略——点电荷电势取无穷远为零势点所以你算出来的 U 是相对值绝对值大小完全取决于 q 和 r 的单位制。如果坐标用米、电量用库仑整张图的 U 落在 1e-9 量级matlab 的 double 有 15 到 16 位有效数字完全撑得住但如果你把 q 写成 1 库仑U 会冲到 1e10surf 的默认坐标轴刻度会开始丢精度图上看着没问题数值其实已经在抖。2.2 meshgrid 出来的 X、Y 到底进了什么[X,Y]meshgrid(-40:0.3:40,-40:0.3:40)这一行生成两个 267×267 的矩阵X 的每一行都一样、沿列方向递增Y 的每一列都一样、沿行方向递增两个矩阵同一个下标位置拼起来就是平面上一个点的坐标。后面所有运算都是逐元素矩阵运算.^2、./这些点运算符不能省省了就变成矩阵乘除轻则报错重则算出一堆毫无意义但能出图的数。区间 -40 到 40、步长 0.3实际点数不是 267 整——(40-(-40))/0.3 除不尽matlab 会取到 39.8 为止共 267 个点最后一行/列并不在 40 上。这个细节在需要跟解析解逐点对比的时候会咬人你以为取的是边缘点实际差了 0.2。我一般会改成linspace(-40,40,267)明确点数或者干脆用能整除的步长让网格端点落在整数上。循环体里 X、Y 从头到尾不变变的只有每次算出来的 Rk 和 Uk。Uk 是一个和 X 同尺寸的矩阵每个元素对应一个网格点的电势贡献累加到 U 上51 次循环之后 U 就是完整的电势分布。整个过程内存占用很小267×267 的 double 矩阵不到 600KB就算把步长压到 0.01、网格放大到 8000×8000单个矩阵也就 500MB 出头普通机器扛得住但循环次数乘以网格点数之后跑一次要等的时间就不是「顺手看一眼」的节奏了。2.3 一张参数表把全局定死这类仿真最容易出的问题不是公式写错而是参数之间互相不匹配——电荷的位置范围跟网格范围对不上、电荷间距跟步长撞在一起、电荷个数改了但每个电荷的电量没跟着改。我习惯在动手写循环之前先把参数列成一张表每个数都写清楚它影响什么改的时候才知道会连带改什么。参数原报告取值物理含义改动后的连带影响q1.6e-18 C原文写作 1.6*10e-19单个点电荷电量直接按比例缩放整张电势图的量级e01e-9/(36π) ≈ 8.842e-12真空电容率与真值 8.854e-12 差 0.14%做解析对照时会体现出来nr50线电荷离散段数决定近场能否分辨单个电荷也决定总电荷量如果 q 不跟着改网格步长0.3空间采样精度太小则计算量暴涨太大则近场直接失真坐标范围-40 到 40观察窗口必须完整包住电荷分布否则图会缺一块电荷位置x k-25k1..51 → -24 到 26线电荷的几何位置和 nr 强耦合改 nr 必须同步改这里这张表里最值得盯的是最后两行。原代码的电荷位置是X-k25里直接写死的k 从 1 走到 nr1位置就从 -24 走到 26跨度 50跟网格的 -40 到 40 只对上了一部分——左边空 16右边空 14图看起来大致居中其实是偏的。而定义电荷分布的那行Rlinspace(0,10,nr1)从头到尾没被引用过属于典型的死变量。这两处放在一起意味着你改 nr 的时候电荷的位置分布根本不会跟着变只是循环次数变了图当然会变但变的原因跟你以为的不是一回事。3. 三张图三种职责surf、contour、quiver 各自的调参与判断标准3.1 surf 看的是全局形状和对称性surf(X,Y,U)把电势当作高度画成曲面它最大的价值不是好看而是让你一眼看出对称性对不对。电荷沿 y0 这条线排布整个电势场关于 y0 严格上下对称如果你的曲面看起来一边高一边低那一定是 Rk 里某一项写反了符号或者网格的 Y 没有以 0 为中心。这种错误在 contour 图上不一定明显但在三维曲面上会立刻暴露。surf 有个固有的视觉缺陷近场电势按 1/r 发散靠近电荷的地方能到 1e-7 量级远处只有 1e-9两者差两个数量级默认的 z 轴自动缩放会把远处压成一张平板看上去「远处没有电场」实际上那里还有起伏只是被近场的尖峰吃掉了。想看清远处的结构要么用zlim手动截断要么对 U 取对数再画——后者在物理上不太严谨但作为观察手段非常好用等位线的疏密关系会变得直观。真正判断 surf 画得对不对我只看三件事曲面是否关于 y0 对称、峰值是否出现在电荷位置上、远处是否平滑单调衰减。这三条都过了说明从公式到循环这一段没写错剩下的问题都在参数取值上。3.2 contour 的等位线数量决定图能不能读原代码用CVlinspace(Vmin,Vmax,30)生成 30 个等位值问题在于 Vmin 和 Vmax 的跨度是两个数量级以上而且是线性均分结果就是绝大部分等位线挤在电荷附近那几个像素里远处一根线都没有。这种图看起来「有内容」实际上中场和远场的梯度信息全丢了而这恰恰是判断场强衰减规律最需要看的部分。常见的改法有两种。一种是把等位值改成对数间隔CV logspace(log10(abs(Vmax)*1e-3), log10(abs(Vmax)), 30)让远场也分到线另一种是先算一遍 U用prctile或者手动排序砍掉最靠前的 1% 极值再线性均分。我倾向前者因为对数间隔和 1/r 的衰减规律天然匹配画出来的等位线间距在视觉上更均匀读场强大小也更直接——等位线越密场强越大这个关系在任何一张等位线图上都成立。还有一个细节contour 默认会画在当前的坐标轴范围里如果前面画过 surffigure 没有切干净等位线会叠在三维图上。原代码里figure(2)是显式指定的这一点做得对只是后面的hold on在循环里每画一个点调用一次虽然不报错但逻辑上是冗余的正确写法是进循环前hold on一次循环结束hold off一次。3.3 quiver 之前必须归一化不归一化会看到什么quiver 画的是矢量场每个网格点上一个箭头长度代表场强大小。如果直接把gradient算出来的 Ex、Ey 丢进去matlab 的自动缩放会按全场最大值来定比例——而近场因为 1/r² 的关系场强能比远场大四五个数量级自动缩放之后远场的箭头全部缩成一个点整张图只剩下电荷附近一团乱麻。这就是必须先把场强归一化的原因。归一化的写法是AEsqrt(Ex.^2Ey.^2); ExEx./AE; EyEy./AE;做完之后所有箭头长度一致只剩方向信息。代价是你看不出场强大小了所以归一化的 quiver 必须和 contour 叠在一起看箭头告诉你方向等位线密度告诉你大小两者配合才是一张完整的场图。这里有个静默出错的点AE 在某些位置可能为 0除零之后 Ex、Ey 变成 NaNquiver 遇到 NaN 会直接跳过不画图上就会出现莫名其妙的空白而且不报任何错。稳妥的写法是AE max(AE, eps)或者干脆把 AE 小于某个阈值的点单独标记出来。另外gradient的第二个、第三个参数是网格间距默认值 1而这里的实际步长是 0.3不写这两个参数的话算出来的场强数值会整体偏小 0.3 倍——归一化之后看不出来但只要你想拿它跟解析解对比这个偏差就会跳出来。3.4 叠加顺序和坐标比例两个决定图能不能用的细节画图顺序我固定成先 contour 画等位线再 quiver 叠箭头最后 plot 把电荷位置标上去全程 hold on。这样等位线在底层不会被箭头盖住电荷标记在最上层最醒目。电荷标记用plot(k-25, 0, ro, MarkerSize, 5, MarkerFaceColor, r)一行一个点比原代码那种一次传两组坐标加两个样式的写法清楚得多出问题也好定位。比顺序更容易被忽略的是axis equal。matlab 默认会按窗口比例拉伸坐标轴x 方向 80 个单位、y 方向也是 80 个单位看起来像是正方形但窗口一旦被拖动过比例就变了。这时候等位线图上的圆会被拉成椭圆你可能会误以为线电荷的场不是轴对称的然后回头去查公式——查半天发现公式没错是坐标轴没设等比例。这类问题属于纯粹的血泪经验不加axis equal的图看上去只是「稍微有点扁」很难第一时间联想到是显示设置。4. 线电荷的离散化取舍nr 从 5 到 1000 到底改了什么4.1 离散化的本质用有限个点电荷代替连续分布真实的线电荷是连续分布单位长度带电量 λ空间任一点的电势要用电场积分算U (λ/4πε₀)·ln(...)是个解析式。matlab 里用循环叠加点电荷本质是把这条线切成 nr1 段每段用一个集中点电荷代替段数越多越接近连续分布。这里的关键是「每段电量怎么取」——如果保持 λ 不变每段电量应该是 λ·L/nr总数 nr1 个点电荷加起来总电量是 λ·L·(nr1)/nrnr 越大越接近 λ·L如果像原代码那样每个点电荷固定 q那总电量就是 (nr1)·q改 nr 等于改总电量物理模型直接变了。这是我认为整个实验里最隐蔽的一个坑。你把 nr 从 50 改成 100本来只想提高分辨率结果总电荷翻了一倍整张图的电势量级整体抬高而你还在纳闷「为什么分辨率提高之后数值也变了」。正确做法是把 λ 和 L 作为输入q 由它们反算q lambda*L/nr这样 nr 只影响离散精度不影响物理量。离散化误差的量级也值得说清楚。对于远处的点整条线电荷可以近似看成一个总电量为 Q 的点电荷离散成 nr 段带来的相对误差是 O((Δ/r)²) 量级Δ 是段长、r 是观察距离所以远处几乎看不出差别。近场就不一样了靠近某个具体点电荷的位置电势主要由最近的那一段贡献你会明显看到一个个独立的球面等位线这正是「离散化」的痕迹。想要近场也平滑只能靠加密 nr或者老老实实用解析式。4.2 不同 nr 的实际表现对比同一套参数、同一个网格步长只改 nr跑出来的东西差别很大。下面这张表是我反复调过之后总结的定性对比数值部分是循环次数的硬计算现象部分是可以直接对着图核对的。nr点电荷个数段长 Δx循环次数近场形态中场与远处56106单个电荷的球面等位线清晰可辨整条线呈明显的「糖葫芦」状已退化为单根线电荷的场与解析解偏差可忽略1011511仍能分辨单个电荷但球面之间开始连成带状同上5051151近场仍有个体特征中场已经连成连续等位线原报告取值与解析解偏差在 1% 量级2002010.25201近场细腻个体特征只在最近一两个网格内可见循环耗时明显上升但仍在几秒内100010010.051001近场已接近连续线电荷但奇异点更尖锐surf 的尖峰更细计算量和内存都上一个台阶收益递减怎么选取决于这个实验要说明什么。如果目的是「演示叠加原理、看清多个点电荷如何合成一根线的场」nr 取 50 是合理的因为近场还保留着个体特征图上能看出「这是很多个点电荷叠出来的」讲原理的时候有东西可指。如果目的是「逼近理想线电荷、跟解析解做定量对比」50 远远不够至少要 500 以上而且必须把近场的奇异点处理掉否则误差主要来自那几个 Inf 点而不是离散化本身。我见过有人为了追求「好看」直接把 nr 拉到 5000结果内存没爆、时间也没爆但 surf 图上出现一根细得几乎看不见的尖刺直冲天花板其余部分全被压平——这时候问题不在 nr在奇异点没处理。加密 nr 会让最近的那个电荷离某个网格点更近Rk 更接近 0Uk 更接近 Inf视觉上就是尖刺越来越尖。4.3 拿解析解做一次对照别只信图好看有限长均匀线电荷在中垂面上距离 r 处的电场强度有闭式解E Q/(4πε₀·r·√(r²L²/4))其中 Q 是总电量、L 是线长。当 r 远大于 L/2 时根号项趋近 r公式退化成点电荷的 Q/(4πε₀r²)。这两条给了两个可以独立验证的锚点。第一个锚点是远场退化。取 r 200、L 50修正因子 1/√(1(L/2r)²) ≈ 1-0.0078也就是说数值解应该比纯点电荷公式小 0.8% 左右。验证的时候要把网格范围单独放大到能盖住 r200 的位置比如 -400:4:400 跑一次在 (0,200) 点取 U跟 Q/(4πε₀·200) 比差在 1% 以内就算通过。这个测试能同时验证总电量算得对不对、距离算得对不对、e0 取值合不合理。第二个锚点是收敛性。同一网格点nr 取 50 和 nr 取 500 算出来的 U 应该非常接近相对差在 2% 以内。如果不收敛而是随 nr 单调变化八成是总电量没保持恒定——也就是 4.1 里说的那个坑。这个测试不需要任何解析式纯靠两次运行对比非常适合放在报告里当「结果可信度」的证据。5. 踩坑排查五处让我重跑实验的细节5.1 常数写成 10e-19量级差了整整十倍现象图能跑出来形态也对但把某个网格点的电势手算一遍发现程序算出来的是手算的 10 倍一开始还以为是 e0 取值的锅换了几种写法都对不上。原因1.6*10e-19在 matlab 里不是 1.6×10⁻¹⁹而是 1.6×10⁻¹⁸。科学计数法的尾数部分本身就带指数10e-19等于 10×10⁻¹⁹也就是 1e-18。这个写法在纸上推导时很容易顺理成章敲进代码就变成另一回事而且 matlab 不会给任何提示图照样出。解决写成q 1.6e-19并且养成习惯——所有量级敏感的常数定义完之后紧跟一行fprintf(q %e\n, q)把实际值打出来看一眼。我现在的做法是常数区统一用1.6e-19这种紧凑写法绝不在尾数里放 10 的倍数。顺带说 e0。1e-9/(36*pi)算出来是 8.842e-12真值 8.854e-12差 0.14%。做教学演示完全够用但如果你要跟解析解做定量对照这 0.14% 会实打实地进入误差项。我的处理是直接用 8.854e-12一行代码的事省得后面解释误差来源。5.2 gradient 默认间距是 1不是你的网格步长现象quiver 图的方向完全正确箭头分布也符合预期但只要把 Ex、Ey 的数值拿出来跟解析场强比就发现整体偏小比例恰好是 0.3 倍左右。原因gradient(F)在不给间距参数时默认按 1 计算差分而实际网格步长是 0.3。差分公式是 (F(i1)-F(i-1))/(2h)h 用错结果就整体差一个系数。解决[Ex,Ey] gradient(-U, 0.3, 0.3)第二个参数是 x 方向间距第三个是 y 方向。这一步做完再做归一化方向和大小才都对得上。这个坑危险的地方在于它被归一化掩盖了。只要做了ExEx./AE系数就被除掉图上一点问题看不出来可能一整年都不会发现数值是错的。我的对策是在报告里保留一小段不归一化的输出单独打印某个点的场强数值跟解析解并排放在一起。图上看着对不等于数算得对这一条我记了很久。5.3 U[0] 能跑通靠的是隐式扩展现象在新版 matlab 上跑得好好的换一台装旧版本的机器直接报维度不匹配另有一次把 U 的初值改成别的形状之后整张图的数值全乱了。原因U[0]是 1×1 标量矩阵UUk里 Uk 是 267×267靠 matlab 的隐式扩展广播成 267×267。第一次循环之后 U 变成了 267×267后面就正常了。这套逻辑「碰巧对」但只要运行环境不支持广播、或者你在循环前不小心把 U 赋成了别的尺寸就会出问题。解决U zeros(size(X))把初值设成和网格同尺寸的全零矩阵循环里U U Uk就是纯粹的逐元素累加不依赖任何扩展规则。这行改动不影响结果但让代码的意图变得明确——U 从一开始就是一张和网格一样大的电势图只是初始值全为零。还有一个衍生问题如果一次脚本里要算多组参数比如扫描不同的 nr忘了在每组开始前把 U 清零电势会一路累加下去图越来越亮而且不报错。我现在的习惯是把整个计算段封成一个函数U 在函数内部初始化调用多少次都互不影响。5.4 网格步长撞上电荷位置surf 直接出一根尖刺现象把步长从 0.3 改成 1本来只是想加快速度结果 surf 图出现一根直冲天花板的白线其余部分被压成一张平板等位线图也变成一堆同心圆挤在一起。原因步长 1 的时候网格点正好包含整数坐标而电荷位置 x k-25 也是整数两者重合的位置 Rk 0Uk Inf。surf 遇到 Inf 会把整个 z 轴范围撑到无穷视觉上就是一根刺加一片平。解决三条路。一是步长取非整数比如 0.3用「错开」的方式回避二是给距离加下限Rk max(Rk, 1e-3)把近场电势截断在一个有限值三是把 Inf 点置为 NaNUk(~isfinite(Uk)) NaNsurf 遇到 NaN 会自动挖洞图上留下一个小坑而不是一根刺。我一般用第二种因为物理上你本来就在意近场的精确值直接挖洞等于放弃那片区域。但加了下限之后要在报告里说明近场那几个点的电势是被截断的不能拿去跟解析解比比较只在中远场做。5.5 改了 nr 电荷位置却没动图缺了一半现象把 nr 从 50 改成 20循环次数确实少了但图明显不对称——左边一堆等位线右边一大片空白看上去像是电荷只分布了半边。原因电荷位置写在X-k25里由 k 直接决定跟 nr 没有函数关系。k 只走到 nr121位置就到 -2420 -4 为止右边 -4 到 40 整片区域没有电荷。而Rlinspace(0,10,nr1)这行定义了电荷分布范围却从头到尾没被引用。解决把位置生成和场计算解耦。先算位置数组xs linspace(-L/2, L/2, nr1)再算每段电量qs lambda*L/nr循环里用Rk sqrt((X-xs(k)).^2 Y.^2)。这样 nr 只控制分辨率位置范围和总电量都由 L、lambda 决定改 nr 不会连带改物理模型。顺手把没用到的R那行删掉交报告前少一个被追问的点。6. 把参数定死一份能自检的线电荷电场仿真配置6.1 可以直接抄的核心循环把前面所有结论收拢核心计算段大概长这样。参数区、计算区、绘图层三段分开改参数只动第一段。% ---- 参数区 ---- L 50; % 线电荷总长度 lambda 1e-9; % 线电荷密度 C/m nr 50; % 离散段数 q0 1.6e-19; % 单个电子电量用于对照量级 e0 8.854e-12; % 真空电容率别用 1e-9/(36*pi) 凑 h 0.3; % 网格步长取非整数以避开电荷位置 % ---- 计算区 ---- xs linspace(-L/2, L/2, nr1); % 电荷位置与 nr 解耦 qs lambda*L/nr; % 每段电量保持总电量恒定 [X,Y] meshgrid(-40:h:40, -40:h:40); U zeros(size(X)); % 显式初始化不用 U[0] for k 1:nr1 Rk sqrt((X-xs(k)).^2 Y.^2); Rk max(Rk, 1e-3); % 近场截断防 Inf U U qs/(4*pi*e0)./Rk; % 标量电势逐点累加 end [Ex,Ey] gradient(-U, h, h); % 间距必须与 h 一致 AE max(sqrt(Ex.^2Ey.^2), eps); % 防除零 Ex Ex./AE; Ey Ey./AE;参数区里lambda和L是物理输入nr和h是数值精度开关q0单独列出来只是方便对照量级不参与计算。计算区第一件事就是由lambda、L、nr反算每段电量这一步做完改 nr 就纯粹是改精度。max(Rk,1e-3)那一行把近场截断在 1 毫米代价是最近的那几个点数值不可信收益是整张图不会因为一个 Inf 崩掉。gradient后面两个h不能省省了数值就错。6.2 两个自检跑完图先别急着截图第一个自检是 y 轴对称性。电荷全在 y0 上U 必须关于 y0 严格对称一行代码就能验err max(abs(U(:) - flipud(U)(:)))结果应该是 0 或者机器精度量级。如果不是说明网格的 Y 没有以 0 为中心或者 Rk 里某一项符号写反了。第二个自检是 nr 收敛性。同样参数跑 nr50 和 nr500 两遍在中场随便挑几个网格点比较 U 的相对差应该小于 2%。如果差很多先回去查qs有没有跟着 nr 变再查xs的范围有没有变。这两个自检加起来不到十行代码但能挡住我见过的绝大多数「图好看但数不对」的情况。做完再截图心里踏实。6.3 想再往前走一步的话离散求和这套写法适合讲原理不适合做精密计算。真要算线电荷的场解析式更省事无限长线电荷 E λ/(2πε₀r)有限长的用 4.3 里那个闭式解一行就能算还不用担心奇异点。离散求和的价值在于它把「叠加原理」这件事变得可见——你能亲眼看到 51 个点电荷的球面等位线怎么连成一条线这个视觉过程是解析式给不了的。回过头看我第一次跑这个仿真的时候最该做的不是把图调好看而是先拿一个点手算一遍跟程序对一下。那次之后我养成了一个习惯任何涉及物理常数的仿真先挑一个解析可算的点做单点验证通过了再画整张图。图是给人看的数是给自己兜底的。这份报告里的代码、参数表和几张对照图我整理在一起了改nr和h就能跑出自己的版本希望帮到你。本文还有配套的精品资源点击获取