MATLAB带电粒子混合电磁场轨迹仿真与ODE求解器实战

简介：一份基于MATLAB的带电粒子在混合场运动仿真模拟实验源码，针对均匀电场、均匀磁场及其叠加等不同混合场情景，能够准确计算并绘制带电粒子的运动轨迹，帮助学习者直观理解电磁场对粒子运动的影响。压缩包内共7个文件，其中两个fig文件用于图形界面布局，两个m文件负责运动方程求解与轨迹可视化，另含一个zip附件、一个txt运行说明及一个md项目介绍，整体仅62KB，轻量便携。目前已有97人学习下载，适合MATLAB初学者、物理或电子信息专业学生及毕业设计人员参考。代码经测试可稳定运行，提供完整源码与资料，既可直接用于课程设计、项目演示，也能在此基础上修改电场磁场参数，扩展研究不同混合场组合下的粒子动力学行为。

1. 带电粒子在混合场中的运动方程与仿真思路

一个常被忽略的事实：带电粒子在真空电磁场中的轨迹，绝大多数没有闭式解。匀强磁场单独作用还能得到圆周运动，叠一个与磁场垂直的匀强电场后，解立刻变成摆线与漂移的组合；场一随空间变化，教科书公式基本失效。MATLAB 仿真实验的价值正是补上这段空白：把运动方程写成六维一阶 ODE，交给求解器数值积分，只要场函数定义准确，不同混合场情景下的轨迹就能可靠重现。这个方向主要服务三类人：做大学物理仿真实验的本科生、准备等离子体或加速器课题的研究生，以及需要快速估算粒子偏转路径的工程师。下面这套方案从洛伦兹力函数写起，覆盖求解器参数、三种混合场建模、轨迹绘图与回归验证，照着敲完就能在本地跑出可对照解析解的轨迹图。

2. 用 MATLAB 把混合场写成洛伦兹力函数，再交给 ODE 求解器

2.1 六维状态向量与比荷归一化：代码的第一步是整理方程

带电粒子在混合场中同时受到电场力 qE 和洛伦兹力 q(v×B)，牛顿-洛伦兹方程是一个二阶常微分方程组：

m dv/dt = q(E + v×B)，同时满足 dx/dt = v。

MATLAB 的 ODE 求解器只接受一阶显式系统，所以先把状态定义为六维列向量 y = [x; y; z; vx; vy; vz]，原来的二阶方程组拆成两个三组一阶方程，右侧恰好是 [v; a]。这种写法不只是迁就 API，它把“场如何随位置和时间变化”隔离到独立的场函数里，后续切换匀强场、梯度场、时变场逻辑非常自然。

E 和 B 在这个框架里建议写成函数句柄而不是普通变量。混合场场景中，电场可能只存在于某个区域，磁场可能沿坐标变化，用函数句柄传参比在 ODE 右端函数内部做 switch 分支干净得多。比荷 qm = q/m 放在外面预先算好，归一化时先令 qm=1 验证几何形状，再换真实粒子的值。

function dydt = lorentz_rhs(t, y, qm, E, B) % 状态向量 y = [x; y; z; vx; vy; vz] r = y(1:3); v = y(4:6); % E、B 是函数句柄，入参为位置 r 和时间 t Eval = E(r, t); Bval = B(r, t); % 洛伦兹加速度 a = (q/m) * (E + v x B) a = qm * (Eval + cross(v, Bval)); dydt = [v; a]; end

这段代码把 r 和 v 拆成 3×1 列向量，cross 要求两个向量尺寸一致，返回的 a 也是 3×1。qm 是预先算好的标量，避免每个时刻重复做除法和类型转换。函数签名里保留 t 这个占位参数，是因为 ode45 会按固定签名调用该函数，即使当前是静态场也不能省略。E、B 两个函数句柄在主脚本里定义，例如 B = @(r,t) [0;0;1] 就是沿 z 轴正向的匀强磁场。

最小可跑通的主脚本如下：

qm = 1.0; % 归一化比荷，先验证轨迹形状 E = @(r, t) [0; 0; 0]; % 零电场 B = @(r, t) [0; 0; 1]; % 匀强磁场，归一化单位 tspan = [0 40]; % 覆盖约 6 个回旋周期 y0 = [0; 0; 0; 1; 0; 0]; % 初始位置原点，速度沿 x，与 B 垂直 [t, y] = ode45(@(t, y) lorentz_rhs(t, y, qm, E, B), ... tspan, y0, odeset('RelTol', 1e-9, 'AbsTol', 1e-11)); plot3(y(:,1), y(:,2), y(:,3));

这段脚本跑完就能看到 x-y 平面内的闭合圆周。tspan 写成两个元素时，求解器按内部误差控制决定输出点密度，不需要外部指定步长。初始速度沿 x、磁场沿 z，洛伦兹力只提供向心力，符合匀强磁场圆周运动的前提。RelTol 放到 1e-9 而不是默认的 1e-3，是因为轨迹多圈后对相位误差敏感，默认容差下弧线不能闭合，看起来像缓慢进动。

2.2 ode45 不是唯一选择：刚性问题与容差参数怎么定

求解器	方法性质	适合的场景	使用提示
ode45	显式 4/5 阶 Runge-Kutta	场函数平滑、时间尺度差异不大	默认首选
ode23	显式 2/3 阶	快速预览、允许较大误差	精度下限较高
ode15s	隐式变阶多步	刚性问题	不受高频回旋制约步长
ode23t	隐式梯形	中等刚性、希望能量追踪较好	稳定但单步开销大

绝大多数混合场仿真用 ode45 就够，尤其是电场磁场只随位置缓慢变化的情况。什么时候换刚性求解器？回旋频率远高于粒子通过整个场区所需时间的时候，例如 ωc = qB/m 达到 1e10 rad/s，而研究对象是秒级宏观漂移。ode45 为了解析每个回旋圈被迫采用极小步长，整体慢到不可接受。这种场景一个选项是改用引导中心近似；坚持做全轨道模拟的话，再尝试 ode15s 或 ode23t。判断是否刚性的快速办法是把 RelTol 降低一个数量级，看求解器耗时是否翻倍以上，如果是，说明步长受限于局部振荡。

容差参数方面，RelTol 是相对误差，AbsTol 是绝对误差下限。轨迹仿真里 RelTol 一般取 1e-8 到 1e-10，AbsTol 再低一到两个数量级。位置和速度分量数值相差很大时，把 AbsTol 写成六维向量分别约束：

opt = odeset('RelTol', 1e-9, ... 'AbsTol', [1e-12; 1e-12; 1e-12; 1e-8; 1e-8; 1e-8]);

六个数字分别对应 y 中的 x、y、z、vx、vy、vz。位置分量通常在零点附近，给太松会导致转向点识别模糊；速度分量可能很大，给太紧会让求解器在每个回旋圈上浪费计算量。这个技巧在质谱仪、偏转磁铁这类对轨迹细节敏感的模拟里很实用。

提示：默认 RelTol 是 1e-3，直接跑轨迹只能演示“有东西在动”。要画出能和解析解对比的闭合圆周、螺旋或摆线，至少把 RelTol 降到 1e-6 以下。

2.3 时间跨度、输出点数与 deval：采样策略决定后续绘图是否卡顿

tspan 写成 [0 40] 时，ode45 返回的点数由误差控制自动决定，常规回旋运动大约几千到几万点，绘图和动画都能承受。如果为了“平滑”把 tspan 写成 0:0.001:40 这种固定向量，实质是强制求解器在每个网格点输出，内存和绘图负担成倍增加，内部积分精度并不会因此提高。我一般在需要固定密度输出时改用结构体输出加 deval 插值：

sol = ode45(@(t, y) lorentz_rhs(t, y, qm, E, B), ... [0 40], y0, odeset('RelTol', 1e-9, 'AbsTol', 1e-11)); ts = linspace(0, 40, 10000); ys = deval(sol, ts).'; % 转置成 10000×6

deval 是在已通过误差检验的解上做插值，而不是重新积分，几乎不增加计算时间。ys 的列顺序与 y 一致，第一到三列是位置，第四到六列是速度。这个输出直接喂给 animatedline 或 comet3 都很合适。

如果模拟粒子轰击平面、强聚焦磁铁出口这类带边界条件的场景，可以用事件函数让积分在指定边界上停止，避免粒子越过边界后场函数返回无效值导致 NaN。事件函数需要返回检测值、是否终止和方向三个量，检测 z 平面的写法如下：

function [value, isterminal, direction] = stop_at_plane(t, y) value = y(3); % 当前 z 坐标 isterminal = 1; % 触发后终止积分 direction = -1; % 只捕获 z 从正到负的穿越 end

用 odeset('Events', @stop_at_plane) 注册后，ode45 会在 z 从正到负穿越 0 的时刻精确结束。这个机制比事后在轨迹数据里找最近点更可靠，也能避免粒子进入没有定义的场区域后轨迹发散。

3. 三种常见混合场情景的参数设计与 MATLAB 建模

3.1 匀强电场与匀强磁场垂直：E×B 漂移先于轨迹验证

E 和 B 垂直时，运动方程不再有简单圆周解，稳态解是 E×B 漂移，漂移速度与电荷符号无关。初始速度为零时，粒子先被电场加速、被磁场偏转，形成摆线轨迹并整体漂移。这个场景适合验证场函数与求解器是否配对正确，因为漂移速度可以直接用叉积公式算出来。

qm = 1.0; % 用结构体保存场函数和初始条件，方便一键切换场景 scene = 'crossed'; % uniformB / crossed / mirror switch scene case 'uniformB' s.E = @(r, t) [0; 0; 0]; s.B = @(r, t) [0; 0; 1]; s.y0 = [0; 0; 0; 1; 0; 0]; % 纯圆周运动 case 'crossed' s.E = @(r, t) [0.2; 0; 0]; % 电场沿 x s.B = @(r, t) [0; 0; 1]; % 磁场沿 z s.y0 = [0; 0; 0; 0; 0; 0]; % 初速为零，观察摆线 case 'mirror' s.E = @(r, t) [0; 0; 0]; s.B = @(r, t) [0; 0; 1 + 0.2 * r(3).^2]; s.y0 = [0; 0; -6; 1; 0; 4]; % 沿 z 方向进入磁镜 end [t, y] = ode45(@(t, y) lorentz_rhs(t, y, qm, s.E, s.B), ... [0 60], s.y0, odeset('RelTol', 1e-9, 'AbsTol', 1e-11));

switch 分支只适合测试阶段，后期建议把场景写进配置文件或 GUI 的下拉框。crossed 场景中 E=[0.2;0;0]、B=[0;0;1]，理论漂移速度为 E×B/|B|² = [0; -0.2; 0]。数值上可以从速度分量做时间平均得到：

v_d_sim = mean(y(end-500:end, 4:6), 1); % 取末段平均 Bval = s.B(zeros(3,1), 0); Eval = s.E(zeros(3,1), 0); v_d_theory = cross(Eval, Bval) / sum(Bval.^2);

为什么取最后 500 个点而不是全程平均？因为粒子的回旋运动贡献了一个大幅振荡项，全程平均会引入相位项；末段已包含数百个回旋周期，振荡项衰减到很小。对比可以发现两者误差通常在 1e-4 量级以下。

3.2 梯度磁场与磁镜效应：只改 Bz(z) 就能演示捕获

带电粒子进入磁场增强区域时垂直回旋动能增大，按绝热不变量 μ = v⊥²/(2B) 守恒，平行动能相应减少，到磁镜点 v∥=0 后反向，形成磁镜之间的往复。这是磁约束装置里最常被仿真的场型之一，也是混合场仿真里最容易出现意外结果的场景，因为初速在 z 方向稍大一点，反射点就会外移很远。

上面的 mirror 分支只让 Bz 随 z 变化，属于一维磁镜模型。这个模型的好处是参数直观、计算快，能演示捕获和反射；代价是单分量磁场不满足三维散度约束，不能直接用于等离子体边界分析。工程上遇到真实磁镜位形时，常见做法是用有限元场数据插值替代这个函数句柄，运动方程部分不用改。

粒子在 z=-6 处以 v=[1;0;4] 出发，B(z)=1+0.2z²，z=0 处磁场最小。绝热不变量给出反射点磁场 B_m = B(z_start) × |v|² / v⊥0²。这里 B(z_start)=8.2，v⊥0²=1，|v|²=17，所以 B_m≈139.4，对应反射点 z≈±26.3。粒子会在约 ±26 的纵向范围内往复。用下面的代码检查 μ 沿轨迹的波动：

vper = hypot(y(:,4), y(:,5)); % B 沿 z 时，垂直速度在 x-y 平面 mu = vper.^2 ./ (2 * (1 + 0.2 * y(:,3).^2)); plot(t, mu);

μ 不是严格常数，绝热不变量只在磁场尺度远大于回旋半径时近似守恒。曲线若出现明显斜坡，通常是 v⊥ 或磁场公式写错；若只在振荡中夹杂缓慢漂移，属于可以接受的绝热误差。把 z 方向初速从 4 改成 10，反射点会外移到约 ±64，纵向范围显著扩大；若想仿真真正的逃逸，需要把磁场改成有限区间剖面，让梯度只出现在局部区域。

3.3 时变电场与回旋共振：频率比扫参看能量增长

混合场不限于静态场，回旋共振加速就是射频电场与匀强磁场的典型组合。沿 x 方向加频率接近回旋频率的电场，粒子每个回旋周期都在同一相位获得能量，轨迹半径逐圈增大；频率失谐时能量增长明显下降。

wc = qm * 1; % qm=1, B=1 时回旋频率 s.E = @(r, t) [0.5 * sin(0.95 * wc * t); 0; 0]; s.B = @(r, t) [0; 0; 1]; s.y0 = [0; 0; 0; 1; 0; 0]; [t, y] = ode45(@(t, y) lorentz_rhs(t, y, qm, s.E, s.B), ... [0 400], s.y0, odeset('RelTol', 1e-9, 'AbsTol', 1e-10)); Ek = 0.5 * (y(:,4).^2 + y(:,5).^2 + y(:,6).^2); semilogy(t, Ek);

0.95 是扫参的起点，改成 1.0 就是精确共振。把 0.95 换成 0.8:0.01:1.2 的循环，记录每个频率比对应的最终能量，得到共振峰曲线。能量曲线上出现周期性小突起不是数值 bug，那是失谐电场产生的边带调制。要注意时变场函数必须足够光滑，方波或三角波会在场函数一阶导不连续处触发误差估计失败，这种场景应拆分时间区间分段积分。

4. 轨迹可视化与诊断：plot3、quiver3、animatedline 的正确打开方式

4.1 画出真实比例：axis equal 与轨迹完整性检查

先画 3D 轨迹再设置坐标轴比例，是很多仿真实验漏掉的一步。MATLAB 默认会拉伸坐标轴，让本来闭合的圆周看起来像椭圆，磁场垂直平面里的摆线也会被压成不可辨认的波形。基础绘图代码如下：

plot3(y(:,1), y(:,2), y(:,3), '.-', 'MarkerSize', 4); hold on; scatter3(y0(1), y0(2), y0(3), 40, 'r', 'filled'); xlabel('x'); ylabel('y'); zlabel('z'); axis equal; grid on; view(3);

axis equal 会在 x、y、z 三个方向使用相同的数据单位长度，view(3) 把视角放到默认三维视角。轨迹完整性检查的一个习惯是看首尾是否闭合：匀强磁场无电场时，初始点和自身重合的程度直接反映数值误差；使用默认 RelTol 时，多圈轨迹往往出现明显偏出圆周的尾巴，此时不要调整绘图属性，先回头改误差容限。

4.2 quiver3 叠加场矢量：画箭头之前先想清楚稀疏度

单独一条轨迹看不出场的分布，叠加场矢量能让轨迹与场结构对应起来。quiver3 的问题是箭头密度容易失控，网格点数上到 20×20 之后整个图变成刺猬，轨迹反而不清楚。我一般把网格控制在 5×5×4 左右，用数据范围做线性采样：

xg = linspace(min(y(:,1)), max(y(:,1)), 5); yg = linspace(min(y(:,2)), max(y(:,2)), 5); zg = linspace(min(y(:,3)), max(y(:,3)), 5); [xx, yy, zz] = meshgrid(xg, yg, zg); % 匀强磁场沿 z 时的场矢量；空间变化场建议写采样函数 quiver3(xx, yy, zz, ... zeros(size(xx)), zeros(size(yy)), ones(size(zz)), ... 0.4, 'Color', [0.6 0.6 0.6], 'LineWidth', 0.8);

quiver3 前三个参数是箭头起点，后三个参数是场矢量分量；最后的 0.4 是全局缩放因子，控制箭头长度占数据范围的比重。对 mirror 这类梯度场，把 ones(size(zz)) 换成 1+0.2*zz.^2，并在调用前把场函数整理成可接受矩阵、返回 Bx、By、Bz 的函数，避免把 ODE 里的函数句柄逻辑重复写一份。

4.3 animatedline 做长轨迹动画：从逐帧重绘改成增量绘图

动画的关键是还原时间维。comet3 函数只需要一行代码，但轨迹长时彗尾重绘开销大，在实时脚本里拖动参数滑块会明显卡顿。animatedline 是增量式绘制，适合几万点以上的轨迹动画。常用写法是：

h = animatedline('Color', [0.1 0.3 0.9], 'LineWidth', 1.2); axis equal; grid on; view(3); for k = 1:10:length(t) addpoints(h, y(k,1), y(k,2), y(k,3)); if mod(k, 200) == 0 drawnow limitrate; % 限制刷新率，防止实时脚本卡死 end end drawnow;

1:10:length(t) 把绘制点稀疏到采样密度的十分之一，视觉上仍足够平滑，动画速度提升明显。mod(k,200)==0 每 200 帧才刷新一次。drawnow limitrate 是较新版本引入的功能，老版本直接换成 drawnow。做磁镜轨迹时，把 z 方向的往复和 v∥、v⊥ 的能量交换放在同一幅图不同子图里，能直观判断捕获点位置。

4.4 诊断：用能量不变量反查数值误差

无电场时粒子动能严格守恒，任何数值误差都会表现为能量曲线的波动或漂移。这个判据比肉眼看轨迹闭合可靠得多，因为轨迹误差在小范围内可能被对称性掩盖。代码如下：

Ek = 0.5 * (y(:,4).^2 + y(:,5).^2 + y(:,6).^2); Ek_rel = (max(Ek) - min(Ek)) / mean(Ek); fprintf('动能相对波动: %.3e\n', Ek_rel);

匀强磁场且 RelTol=1e-9 时，Ek_rel 通常在 1e-8 以下；如果出现 1e-4 量级的波动，先检查是否把电场错误地设成了非零值，再检查 tspan 的时间单位是否与场定义一致。时变电场存在时动能不再守恒，不要用这个判据，改用 3.2 节里的绝热不变量 μ 做诊断。很多大学物理可视化实验会把这段代码直接放进实时脚本，用滑条调 RelTol 观察能量波动，是个很直观的误差演示。

5. 回归验证与参数标定：让仿真轨迹经得起实验对照

5.1 三个回归用例

仿真代码改场函数时容易把正负号、比荷映射改错。固定三个零争议用例，每次修改先跑一遍：

用例	场设置	解析预期	数值验证方式
匀速直线	E=0，B=0	位移=v0·t	任意两点位移/时间差
回旋圆周	E=0，匀强 B	T=2π/(qm·B)	vx 峰值间隔
E×B 漂移	E⊥B 匀强	v_d=E×B/B²	末段速度平均

回旋周期的数值验证不依赖额外工具箱，vx 在匀强磁场中按正弦变化，局部峰值间隔就是一个回旋周期：

locs = find(diff(sign(diff(y(:,4)))) < 0) + 1; % vx 局部峰值 T_sim = mean(diff(t(locs))); T_theory = 2 * pi / (qm * 1); % 归一化 qm=1, B=1 rel_error = abs(T_sim - T_theory) / T_theory;

这里用二阶差分方向变化定位峰值，结果与 findpeaks 基本一致；峰值间隔取平均可抵消相邻周期之间的相位抖动。

5.2 从归一化参数换到真实物理单位

几何验证跑通后换真实单位，需要同步调整比荷、磁场量级和 AbsTol。以质子为例，qm=9.5788e7 C/kg，B=1 T 时回旋周期约 6.56e-8 s：

qm = 9.5788e7; % 质子比荷 T_cyc = 2 * pi / (qm * 1); tspan = [0 20 * T_cyc]; y0 = [0; 0; 0; 1e6; 0; 0]; % 1e6 m/s 垂直入射 opt = odeset('RelTol', 1e-9, ... 'AbsTol', [1e-8; 1e-8; 1e-8; 1; 1; 1]);

位置 AbsTol 取 1e-8 m，速度 AbsTol 取 1 m/s，适用毫米级偏转轨道；长距离漂移仿真可把位置容差放宽到 1e-6 m。

5.3 三个容易耽误事的坑

坐标轴比例未固定是轨迹图最常被挑剔的问题，圆变椭圆、摆线变形都源于默认 3D 坐标轴长宽比。输出点过多导致绘图内存占用高时，优先用 deval 插值而不是缩小 tspan。场函数返回 NaN 时不会立即报错，轨迹会在某点截断，排错统一用 any(isnan(y(:))) 定位。初始速度完全平行于磁场时洛伦兹力为零，轨迹是直线，这本身不是仿真错误；把它写进回归用例反而能验证坐标正负号处理。把三个用例包成一个独立函数，每次改完场建模先跑回归，再画轨迹，跟解析解对照的返工时间能省一大半。

本文还有配套的精品资源，点击获取