简介这份仿真实验报告面向电子信息、电气工程及物理相关专业的学生与Matlab初学者围绕线电荷周围电场与电位分布的可视化展开可用于电磁场课程实验参考、课程设计入门或自学复现。文档完整梳理实验目的、库仑定律与场强/电势公式、网格点设置、电势累加流程及绘图实现并附可运行的Matlab代码便于对照理解数值计算与图形绘制。压缩包仅含1个docx文件大小约227KB结构紧凑包含实验原理、程序流程、代码、结果图与结论感悟等部分。读者可重点参考其中用surf绘制电势三维图、contour画等位线、quiver叠加电场矢量的做法以及矩阵维度一致、初始电势赋零、线电荷密度选取等排错与调参思路。目前已有3041人学习适合作为短周期实验报告模板和线电荷电场仿真入门材料。1. 一份能跑通的 matlab 线电荷场分布实验报告卡住人的从来不是公式电磁场课设里“线电荷周围的电场与电位分布”是个经典题目公式在纸上十分钟就能推完真正让人熬到半夜的是 matlab 里那三十几行代码电位算出来是一根通天柱子、等位线画成一团乱麻、quiver 箭头根根一样长。这份资源是一套完整的仿真实验报告实验目的、原理推导、程序流程、可运行代码、三张结果图和结论讨论都在里面核心代码只有三十多行 matlab。它适合三类人第一次做电磁场数值仿真的同学、要给课程设计找一份能改参数的模板的人以及想用最小代价验证“离散点电荷叠加能不能逼近连续线电荷”这件事的工程师。下面按“物理量对到变量名 → 画图分工 → 踩坑排查 → 拿解析解对账”的顺序拆。2. 从库仑定律到 51 次循环把物理量一个个对到变量名上2.1 为什么先叠电势再求梯度而不是直接叠电场点电荷的电势是标量$U_i q_i/(4\pi\varepsilon_0 r_i)$N 个电荷的总电势就是 N 个标量相加没有方向、没有分量符号写起来几乎不会错。电场是矢量直接叠加要拆成 $E_x$、$E_y$ 两个分量分别累加每个分量还得带上 $\Delta x/r^3$、$\Delta y/r^3$ 这样的方向因子符号一错整张图就反了。所以常见做法是先算总电势矩阵U再用中心差分E -∇U反推场强代价是损失一点精度后面第 5 章会量化这个损失。把代码变量和物理量对齐是读这段程序的第一步符号代码变量物理含义单位$q$q每个离散点电荷的电量C$\varepsilon_0$e0真空介电常数取 8.854e-12F/m$r_i$Rk场点到第 i 个电荷的距离m$U$U叠加后的总电位V$E_x,E_y$Ex,Ey场强两个分量V/m$1/(4\pi\varepsilon_0)$ 约等于 8.99e9原代码写成q/4/pi/e0注意这是从左到右的除法等价于 $q/(4\pi\varepsilon_0)$不是 $q/4 \cdot \pi/\varepsilon_0$别自己改写的时候改错结合顺序。线电荷的处理思路是把一段长度 $L$、线密度 $\lambda$ 的带电直线切成 N 段每段等效成一个点电荷电量 $q \lambda L / N$。N 越大越接近连续分布但计算量线性增长这是后面调参的主线。2.2 网格、预分配与内存U[0]为什么能跑却是个坏习惯先看坐标网格这一段的写法h 0.3; % 网格步长单位 m xm -40:h:40; % x 方向采样点共 267 个 [X, Y] meshgrid(xm, xm); % 生成 267×267 的坐标矩阵 U zeros(size(X)); % 按网格尺寸预分配别写成 U[0]meshgrid输出的是矩阵而不是向量X的每一行都是xm的复制Y的每一列都是xm的复制这样X(i,j)、Y(i,j)就直接对应平面上第 (i,j) 个点的坐标后面所有运算都是矩阵整体运算不用双重 for 循环。267×267 约等于 7.1 万个点一个 double 矩阵约 570 KBU、Ex、Ey、AE加起来不到 3 MB普通笔记本毫无压力。原代码里的U[0]第一次执行U U Uk时会被广播成 267×267结果没错但有两个隐患一是每次累加都要重新分配内存51 次循环白白多花时间二是如果循环体一次都不执行比如把nr设成 -1 或者提前 returnU还是个标量后面的min(min(U))、contour直接报维度错。养成zeros(size(X))的习惯这类玄学报错能少一半。2.3 循环体里的坐标平移电荷到底摆在哪原代码这一行是整段程序的核心Rk ((X - k 25).^2 Y.^2).^0.5; % 第 k 个电荷在 (k-25, 0) Uk q/4/pi/e0 ./ Rk; % 注意这里是 ./ 逐元素除 U U Uk;把X - k 25整理一下就是X - (k-25)所以第 k 个电荷的坐标是(k-25, 0)k 从 1 到nr1也就是 51x 坐标从 -24 到 26间距 1 m。这里有两个容易忽略的点线电荷的中心在 x1 而不是原点图像会整体右偏另外程序里那行R linspace(0,10,nr1)算完之后从来没被引用过是从别的版本残留下来的一句改参数的时候别被它误导。我一般会重写成下面这样中心在原点、长度和线密度都能显式控制N 51; % 离散段数越大越接近连续线电荷 L 50; % 线电荷长度m lam 1e-9; % 线电荷密度C/m q lam * L / N; % 每段等效点电荷电量 x0 linspace(-L/2, L/2, N); % 电荷位置中心在原点 rmin h/2; % 奇点截断半径避免 r0 for k 1:N r hypot(X - x0(k), Y); % hypot 比 sqrt(x^2y^2) 数值更稳 r max(r, rmin); % 截断防止出现 Inf U U q/(4*pi*e0) ./ r; endhypot内部做了溢出保护虽然这里数值不大用不上但换成大坐标范围时sqrt(x^2y^2)会有精度损失。rmin h/2的取法是有讲究的网格上离电荷最近的距离最大就是半个步长把小于h/2的距离统一截断相当于把点电荷“摊”到一个网格单元里物理上对应线电荷有有限半径做出来的图不会有一根刺。2.4 参数表改哪几个数能让图和报告对得上变量作用建议取值改大/改小的后果N线电荷离散段数21 / 51 / 101 三档对比太小图像呈锯齿太大循环变慢但结果趋稳h网格步长0.3~1.0减小一半点数变 4 倍内存和时间同步翻倍L线电荷长度20~200相对场点距离太短时退化不成无限长线电荷lam线电荷密度1e-9~1e-8 C/m只影响幅值不影响图形形状rmin奇点截断半径h/2过大削平近场过小出现 Inf参数扫描是写报告时最省事的一招固定其余参数只改N出三张图拼在一起讨论“离散化误差随 N 怎么收敛”比空写一段结论有说服力得多。2.5 量级对不上1.6*10e-19不是电子电量原代码写q1.6*10e-19。matlab 里10e-19等于 $10 \times 10^{-19} 10^{-18}$所以这一句算出来是1.6e-18比元电荷大了整整 10 倍。图形的形状完全不受影响整体等比缩放但报告里如果写“取单位电荷”数据就对不上了。正确的写法是1.6e-19指数位不要画蛇添足加个 10。更实际的问题是数量级观感。51 个 1.6e-18 C 的电荷铺在 50 m 上总电势大概在 1e-7 V 量级surf 出来的曲面高度全是 0.0000001 这种数字z 轴标签写U(k)更让人一头雾水。我的习惯是直接给定线密度lam 1e-9让电位落在个位到几十伏之间图上标U (V)报告里注明“线密度取 1e-9 C/m实际物理量按比例缩放即可”。这样读者一眼能看出近处电位高、远处低不用去数小数点后面有几个零。注意q/4/pi/e0里的./和/不能混。Uk q/4/pi/e0./Rk中只有最后一步对Rk用逐元素除前面全是标量运算顺序写反会变成矩阵除矩阵报错信息非常难懂。3. 三张图的分工surf 看趋势、contour 看等位、quiver 看方向3.1 surf三维电位面最容易变成一根针figure(1); surf(X, Y, U, EdgeColor, none); shading interp; % 去掉网格线曲面更平滑 colormap parula; view(35, 30); % 方位角、俯仰角比默认视角好看 xlabel(x (m)); ylabel(y (m)); zlabel(电位 U (V)); zlim([0 prctile(U(:), 99)]); % 关键截掉顶部 1% 的尖峰 grid on; box on;surf默认会给每个网格面画黑边7 万个面叠在一起就是一片黑EdgeColor,none或者shading interp必须加一个。真正让图“变成一根针”的是电荷附近电位按 $1/r$ 发散离电荷最近的那个网格点电位可能是远处的上千倍z 轴自动缩放之后其余部分全被压在底部。zlim用手工上限或者prctile(U(:), 99)截掉最高的 1%视觉上立刻正常。老版本 matlab 用caxis控制色标范围新版本改成了clim写报告时按自己装的版本选。3.2 contour等位线画多少条取决于电位怎么分布figure(2); Vmin min(U(:)); Vmax max(U(:)); CV linspace(Vmin, Vmax, 30); % 30 条等位线 contour(X, Y, U, CV, LineWidth, 0.8); hold on; plot(x0, zeros(size(x0)), r., MarkerSize, 10); % 标出电荷位置 hold off; axis equal; xlabel(x (m)); ylabel(y (m));原代码用linspace(Vmin,Vmax,30)均分 30 条问题是 $1/r$ 型电位的动态范围太大均分之后大部分等位线挤在电荷附近半径一两米内远处的等位线稀疏得看不出形状。两种改法一是手工指定电平比如CV [0.5 1 2 5 10 20 50]按对数分布取二是干脆先对U取log10再画等位线做出来的图更符合“电位随距离对数下降”的物理直觉。axis equal容易被忘掉不加的话等位线会被拉成椭圆看起来像各向异性介质其实是坐标轴比例不一致。3.3 gradient 求场强那个默认间距 1 是个陷阱[Ex, Ey] gradient(-U, h, h); % h 必须传进去 AE hypot(Ex, Ey); % 场强幅值gradient(F)不给间距时默认按 1 处理也就是把相邻网格点当成相距 1 m。但实际步长是 0.3 m算出来的场强整体偏大 $1/0.3 \approx 3.33$ 倍。方向不受影响正数缩放不改变向量方向所以归一化之后再画箭头肉眼完全看不出来只有和解析解一对账才露馅——这是最典型的“图看着对、数全错”。gradient(-U, h, h)里两个h分别对应 x 方向和 y 方向的间距。这里两个方向步长一样所以无所谓如果以后做非均匀网格比如 x 方向 0.1、y 方向 0.5两个参数写反了场强方向会明显歪掉。负号来自 $E -\nabla U$别漏。3.4 quiver 降采样7 万个箭头等于没画step 8; % 每 8 个点取一个 idx 1:step:size(X,1); Exs Ex(idx, idx); Eys Ey(idx, idx); AEs hypot(Exs, Eys); quiver(X(idx,idx), Y(idx,idx), Exs./AEs, Eys./AEs, 0.6, Color, [0 0.55 0]); hold off;267×267 个网格点上每个都画一个箭头屏幕上就是一团绿色噪点。按 8 取一个剩下 34×34 约 1156 个箭头方向和疏密都看得清。除以AEs是归一化只保留方向不保留大小好处是远近箭头一样长不会互相遮挡坏处是读者看不出哪里场强更大。要同时表达方向和大小时有两个办法一是去掉归一化直接传Exs, Eys让 quiver 自动缩放代价是远处箭头短到看不见二是把AEs映射到颜色上用quiver的Color参数配合colormap箭头方向和强弱一图全有。quiver第五个参数 0.6 是缩放系数设成 0 表示自动缩放网格范围一改就得重新试我一般直接写死一个 0.5~1.0 之间的数出图更稳定。3.5 出图分辨率不够报告里全是糊的exportgraphics(gcf, E_potential.png, Resolution, 300);exportgraphics从 R2020a 开始提供直接按分辨率导出比saveas清晰得多。老版本用print(gcf, -dpng, -r300, E_potential.png)。三张图分别导成 300 dpi 的 png插进报告里放大看箭头边缘还是干净的用默认的saveas存出来只有 72 dpi缩放到 A4 宽度就糊了。提示图窗的Position属性会直接影响导出比例set(gcf,Position,[100 100 800 600])固定成 4:3三张图拼在一起才不会一张扁一张长。4. 排查清单五个把结果画成“黑匣子”的坑4.1 图上一根通天柱、contour 直接报错现象surf出来只有一根尖刺其余部分贴地contour报Input must be finite或者画出来只有一条线。原因网格步长 0.3 是从 -40 开始按 0.3 递增的落在网格节点上的 x 坐标满足 $x \equiv 2 \pmod 3$。原代码里 51 个电荷的 x 坐标是 -24 到 26 的整数其中 -22、-19、…、26 一共 17 个正好压在网格节点上Rk算出来是 0Uk变成Inf整个U矩阵被污染。Vmax max(max(U))得到Inflinspace(Inf, Inf, 30)全是NaN。解决循环体里加一句r max(r, rmin)rmin取h/2。或者改网格起点让电荷永远不落在节点上比如xm -40.15:0.3:40但这个方法不通用换参数又会撞上。4.2 电荷标记只显示“”看不到“o”现象plot(k-25,0,ro,k-25,0,r)跑完图上只有一个加号圆圈去哪了原因这一句在同一个坐标上画了两组标记先画圆圈再画加号加号盖在圆圈正中间视觉上只剩下一堆加号。另外hold on写在循环体里第一次循环时前面还没有图形对象实际效果是等到第二次循环才生效。解决一个位置只画一种标记plot(x0, zeros(size(x0)), ro, MarkerFaceColor, r)一条语句传向量比循环里逐个 plot 快得多也干净得多。hold on提到循环外。4.3 箭头根根一样长看不出哪里场强高现象quiver 图上所有箭头长度一致近处远处毫无区别。原因为了美观做了Ex Ex./AE归一化幅值信息被整体除掉了AE算完就丢掉。解决要么保留幅值交给 quiver 自动缩放要么把AE映射成颜色。做报告时我更倾向后者——箭头方向清晰、颜色表达强弱一张图两个信息量比两张图并排更省版面。4.4 场强数值和手算差了 3.33 倍现象取 y 轴上一个点用 $E q/(4\pi\varepsilon_0 r^2)$ 手算得 8.99 V/m程序里读出来 29.9 V/m比值刚好 3.33。原因gradient(-U)没传网格间距默认按 1 m 做差分而实际步长是 0.3 m$1/0.3 3.33$。解决gradient(-U, h, h)。这个错误在验证环节才会暴露光看图永远发现不了所以第 5 章的解析解对账一定要做。4.5 改了nr图像就歪线电荷不在中间现象把nr从 50 改成 20线电荷整体右移等位线明显不对称。原因电荷坐标写死成k-25k 的范围是1:nr1nr一变电荷区间的起点和终点同时变化中心位置跟着漂和网格的固定范围 -40~40 对不上。解决电荷位置用linspace(-L/2, L/2, N)显式给出长度L单独一个变量改N只影响离散密度不影响几何中心。注意排查顺序建议从“数据有没有 Inf/NaN”开始再看“数量级对不对”最后才怀疑画图函数。绝大多数“图画不出来”其实是U矩阵本身已经坏了调contour的参数是白费功夫。5. 数值解准不准拿无限长线电荷的解析解当场对账5.1 解析解长什么样无限长线电荷线密度 $\lambda$距离轴线 $r$ 处的场强$$E_r \frac{\lambda}{2\pi\varepsilon_0 r}$$电位是 $U(r) -\frac{\lambda}{2\pi\varepsilon_0}\ln r C$随 r 对数发散所以电位只能取相对参考点这也是为什么仿真里电位图的绝对数值没有物理意义只有形状有意义。有限长线段长度 L在中垂面上距离轴线 r 处的场强是可以精确写出来的$$E_r \frac{\lambda L}{4\pi\varepsilon_0 r\sqrt{r^2 (L/2)^2}}$$当 $r \ll L/2$ 时根号里 $r^2$ 可以忽略式子退化成 $\lambda/(2\pi\varepsilon_0 r)$和无限长一致。这个条件就是判断“我的有限长模型能不能当无限长用”的判据。5.2 两条算场的路线从 U 求梯度还是直接叠 E路线一原代码走的先叠电势再gradient。优点是代码短缺点是误差来自两处——中心差分本身的 $O(h^2)$ 截断误差以及电荷附近的电位剧烈变化导致梯度失真。路线二直接对每个点电荷叠加库仑定律的矢量形式$E \sum k q_i (\vec{r} - \vec{r}_i)/|\vec{r}-\vec{r}_i|^3$。代码多十几行但没有差分误差近场精度好得多。% 沿 y 轴取剖面两条路线都对一遍 r [1 2 5 10 15]; k0 1/(4*pi*e0); Eana lam*L ./ (4*pi*e0 .* r .* sqrt(r.^2 (L/2).^2)); % 有限长解析解 Enum_g interp2(X, Y, hypot(Ex, Ey), 0, r); % 路线一梯度 Ex_d zeros(size(r)); Ey_d zeros(size(r)); % 路线二直接叠加 for k 1:N dx 0 - x0(k); dy r; d max(hypot(dx, dy), rmin); Ex_d Ex_d k0*q*dx./d.^3; Ey_d Ey_d k0*q*dy./d.^3; end Enum_d hypot(Ex_d, Ey_d); err_g abs(Enum_g - Eana)./Eana; err_d abs(Enum_d - Eana)./Eana; disp(table(r, Eana, Enum_g, Enum_d, err_g, err_d));interp2直接吃meshgrid产生的X、Y不用手动构造网格向量。取点选在 y 轴上是为了避开电荷本体r 从 1 起步也是这个原因——r 太小时有限长和无限长的差别、离散化的锯齿都会混进来误差分析反而看不清主因。5.3 误差从哪来往哪调误差来源随什么变化该往哪调典型现象离散化N越大越小$O(1/N^2)$ 量级加大N或缩短L近场呈波纹等位线有锯齿差分截断步长h越小越小减小h或用直接叠加梯度路线整体偏移形状对数值不准奇点截断rmin过大时近场被削平rmin h/2r 在 1 以内的点误差突然变大有限长近似当 $r$ 接近 $L/2$ 时误差上升加长L或换解析式远场比无限长解偏小做报告时把这张表和一两组实测误差一起放进去比写一段“仿真结果与理论基本吻合”有价值得多。5.4 加长 L 看收敛顺手验证“无限长”假设把L从 20 依次调到 50、100、200固定N101盯住同一个场点比如 r5的场强变化。r5 时 $L/2 10$比值 0.5 已经不算小有限长解比无限长解明显偏小L加到 200 之后 $r \ll L/2$两者差距进入个位数百分比。这个收敛过程可以直接画成一条曲线正好回答报告里“为什么取 L50”这个必答题。提示验证脚本和主程序用同一套X、Y、x0不要另起一套网格。两边参数不一致的时候误差里混进的是参数差不是算法差排查会走歪。6. 把脚本收成一个函数任意电荷分布都能套进去把这段代码重构成函数下次碰到点电荷对、环形电荷、平行板只换输入不换逻辑function [U, Ex, Ey] chargeField(x0, y0, q, X, Y, rmin) % x0,y0,q 为等长向量描述任意离散电荷分布 % X,Y 为 meshgrid 输出rmin 为奇点截断半径 e0 8.854e-12; k0 1/(4*pi*e0); U zeros(size(X)); Ex zeros(size(X)); Ey zeros(size(X)); for k 1:numel(q) dx X - x0(k); dy Y - y0(k); d max(hypot(dx, dy), rmin); U U k0*q(k) ./ d; Ex Ex k0*q(k)*dx ./ d.^3; Ey Ey k0*q(k)*dy ./ d.^3; end end函数里直接把E算出来不再依赖gradient省掉步长换算这个坑rmin作为入参显式传边界条件不再藏在循环体里。线电荷只要一句[U,Ex,Ey] chargeField(x0, zeros(size(x0)), q*ones(size(x0)), X, Y, h/2)就能调起来同轴线、电偶极子、平行板电容也照这个套路喂电荷坐标和电量。每次交图之前我有三个固定动作一个都不跳先any(~isfinite(U(:)))确认没有Inf/NaN再打印max(U(:))和min(U(:))看数量级是不是眼熟最后挑一个场点用解析解手算一遍比对。这三步加起来不到两分钟能挡掉前面提到的绝大多数返工。血泪经验是有一回赶着交报告跳过了第一步U里混进Infsurf画出来是一根针我在contour的参数上折腾了一个多小时才发现数据本身就是坏的。从那以后我每次出新图之前都强制走一遍这三步宁可多花两分钟也不想再体验一次对着屏幕猜黑匣子的感觉。希望帮到你。本文还有配套的精品资源点击获取