简介这是一份基于MATLAB实现的带电粒子在混合场中运动的仿真模拟实验资源能够针对不同混合场配置准确绘制粒子运动轨迹。资源主要面向电磁场与MATLAB仿真方向的本科生也适用于毕业设计、课程设计或相关课题的前期演示适合有一定MATLAB基础并希望快速搭建完整仿真界面的学习者。包内共7个文件包括2个fig界面文件、2个m源码脚本、1个txt说明、1个md文档以及1个zip附带数据包整体压缩后仅62KB结构紧凑、便于下载运行。已有97人学习浏览源码经过功能测试可正常运行代码中保留了清晰的GUI启动入口读者可直接运行观察不同工况下的带电粒子轨迹也可结合自身课题修改参数或扩展功能。整体可用于理解电磁混合场运动规律、快速生成可视化实验结果是一份现成参考。1. 带电粒子在混合场中解的不是轨迹而是力仿真一个带电粒子在混合场里的运动轨迹最容易犯的错误不是选错积分器而是把电场当作速度的直接增量去累加。混合场和单一匀强磁场最大的区别在于洛伦兹力中的 v × B 项让加速度与当前速度方向纠缠在一起电场把能量送进系统磁场又把速度方向掰弯轨迹因此变成一类“看着长得很像但物理来源完全不同”的螺旋和漂移曲线。这篇内容围绕带电粒子在电场与磁场叠加、梯度磁场、时变电场三类混合情景讲清楚怎么把问题写成六维一阶 ODE、如何在 MATLAB 里用 ode45 稳定求解并画出带电粒子的运动轨迹以及怎么确认画出来的曲线不是数值误差伪造出来的。目标读者是做过基础 MATLAB 计算、想快速搭建带电粒子仿真模拟实验的人。2. 混合场的数学描述洛伦兹力与运动方程的归一化2.1 洛伦兹力是状态相关的力场物理学上带电粒子在电磁场中的受力由洛伦兹力给出F q(E v × B)磁场项 q(v × B) 永远垂直于速度本身这意味着磁场对粒子不做功电场项 qE 才负责改变粒子动能。混合场的含义是 E、B 同时存在并且各自可以随空间、时间变化。处理这类问题时第一步不是画图而是把受力写清楚。把 F ma 改写成加速度形式dv/dt (q/m)E (q/m)(v × B)这里最值得留意的是磁场项的系数同时也是回旋角频率的源头。定义 ω_c |q|B/m 后会发现整个运动方程里粒子身份只通过 q/m 这一个比值出现这就给下一步的归一化提供了条件。混合场按场函数性质分成三类建模手法各有侧重混合方式场函数形式仿真注意点匀强电场 匀强磁场E、B 为常矢量存在精确解析解适合验证积分器空间非均匀磁场B B(x)粒子可能被磁镜反射需要记录全程位置时变电场 稳恒磁场E E(t)非自治 ODE能量注入对步长敏感2.2 归一化到回旋单位系避开数值病态直接拿 SI 单位写代码会带来三个问题位置量级在 10⁻³速度量级在 10⁶磁场量级在 1 附近ode45 内部误差控制会偏向数值绝对值大的分量相对误差反而失真。我一般会在动手前把整套方程归一到回旋单位系。选取三个基本尺度参考磁场强度 B₀参考速度 v₀回旋角频率 ω_c qB₀/m由此导出两个派生尺度回旋半径拉莫尔半径 r_L v₀/ω_c电场尺度 E₀ v₀B₀变量替换关系如下物理量SI 写法归一化变量换算关系时间tt′t′ tω_c位置xx′x′ x/r_L速度vv′v′ v/v₀电场EE′E′ E/(v₀B₀)磁场BB′B′ B/B₀这个变换做一次之后运动方程变成看起来极其干净的形式dx′/dt′ v′ dv′/dt′ E′ v′ × B′q/m 在这个方程里完全消失。这样做的好处是同一个算例只要重新换算 ω_c 和 r_L就能从电子迁移到质子代码主体不需要改这也是带电粒子仿真模拟实验里标准的建模顺序。2.3 把二阶运动方程组装成六维一阶 ODEMATLAB 的 ode45 只处理一阶常微分方程组所以要把位置的二阶方程降阶。取状态向量 y [x, y, z, vx, vy, vz]ᵀ并让场函数以“位置 时间”为输入、矢量场为输出function dydt particle_motion(t, y, E_fun, B_fun) r y(1:3); % 当前三维位置 v y(4:6); % 当前三维速度 E E_fun(r, t); % 调用电场函数句柄 B B_fun(r, t); % 调用磁场函数句柄 a E cross(v, B); % 归一化后的洛伦兹加速度 dydt [v; a]; % 位置导数 速度速度导数 加速度 end这段代码有三个关键点。第一场函数统一走函数句柄 E_fun、B_fun之后的匀强场、非均匀场、时变场都通过替换函数实现运动方程文件不再改动。第二cross(v, B) 的参数顺序必须保持 v 在前调换方向会让整个轨迹反转。第三如果后续要在混合场里加入重力或阻尼只需要在加速度行追加项a E cross(v, B) g drag_coeff * v;运动方程的结构保持不变。这样组织以后仿真程序的主体就只做一件事给定场函数求六维状态随时间的变化。3. 用 MATLAB 跑通最小实现:场函数、ode45 与轨迹绘图3.1 用函数句柄定义一个匀强电磁混合场先搭一个能跑的最小案例。取归一化单位制设磁场沿 z 方向强度为 1电场沿 x 方向强度为 0.3粒子初始速度沿 y 方向。在 MATLAB 里按下面的方式定义场函数E_fun (x, t) [0.3; 0; 0]; % 恒定电场沿 x B_fun (x, t) [0; 0; 1]; % 恒定磁场沿 z注意 E_fun 的输入参数虽然没用到 x 和 t但必须保留这两个输入位因为 particle_motion 里调用方式是 E_fun(r, t)接口要保持一致。如果你把某个场函数写成只接收一个参数调用时会直接报错这是刚切到函数句柄方案时最常碰到的问题。3.2 配置 ode45 的容差与时间跨度主程序这样写t_span [0, 12*pi]; % 仿真 6 个回旋周期 y0 [0; 0; 0; 0; 1; 0]; % 位置在原点速度沿 y 方向 opts odeset(RelTol, 1e-8, AbsTol, 1e-10); [t, y] ode45((t, y) particle_motion(t, y, E_fun, B_fun), ... t_span, y0, opts);时间范围取 12π 而不是任意一个整数值原因是归一化后一个完整回旋周期恰好是 2π12π 覆盖 6 圈既能看清螺旋结构又不至于因为轨迹太长而让视图拥挤。RelTol 和 AbsTol 分别设到 1e-8 与 1e-10这是轨迹类问题比较保险的组合。MATLAB 默认的 1e-3 相对容差在长时域积分里会让粒子的相位逐渐漂移画出来螺旋会慢慢散开看起来像是有磁场不均匀性实际却是数值误差。3.3 用彩色散点而不是纯线条画运动轨迹很多初学的人直接用 plot3 画轨迹线这在 MATLAB 可视化大学物理轨迹题时够用但看不出粒子的速度变化。用 scatter3 把速度模映射到颜色上一张图就能同时呈现轨迹形态和能量分布speed sqrt(sum(y(:,4:6).^2, 2)); % 逐点计算速度模 figure(Color, w); scatter3(y(:,1), y(:,2), y(:,3), 6, speed, filled); xlabel(x / r_L); ylabel(y / r_L); zlabel(z / r_L); cb colorbar; ylabel(cb, |v| / v_0); axis equal;scatter3 的 6 是点大小speed 是颜色映射值filled 让空心圆圈变成实心点。axis equal 这一行容易漏但非常重要三个坐标轴刻度不一致时圆形回旋轨迹会被压成椭圆你会误判粒子运动形态。如果要做动画演示可以使用 comet3(y(:,1), y(:,2), y(:,3))但 comet3 每次播放都重新计算尾迹批量比较多个算例时不推荐适合单独展示时用。3.4 场函数复杂时切换到文件函数匿名函数句柄适合常数场和简单表达式。当混合场涉及分段定义、多重参数时文件函数可读性和可调试性更好function E electric_field(x, t) E [0.3; 0; 0]; end function B magnetic_field(x, t) r2 x(1)^2 x(2)^2 x(3)^2; B [0; 0; 1 0.05 * r2]; end调用方式相应改成[t, y] ode45((t, y) particle_motion(t, y, electric_field, magnetic_field), ... t_span, y0, opts);用函数句柄传场函数的最大优点是后续只要在 electric_field 或 magnetic_field 里改返回值就能切换任意混合场情景不会碰到运动方程那部分代码。调试时还可以单独调场函数、检查输出向量维度不需要每次跑完整积分。4. 三类混合场情景的参数设置与轨迹合理性判据4.1 匀强 B 匀强 E用 E×B 漂移做定量验证匀强电场叠加匀强磁场是唯一有精确解析解的混合场所以它是最合适的验证算例。物理上粒子在垂直交叉的电磁场中会获得一个持续漂移速度v_E (E × B) / B²取 E [1,0,0]ᵀB [0,0,2]ᵀ理论漂移速度就是沿 y 方向的 0.5。仿真设置如下E_fun (x,t) [1; 0; 0]; B_fun (x,t) [0; 0; 2]; y0 [0; 0; 0; 0; 0.5; 0]; % 初始速度沿 y与漂移方向一致 [t, y] ode45((t,y) particle_motion(t,y,E_fun,B_fun), ... [0, 40*pi], y0, odeset(RelTol, 1e-9, AbsTol, 1e-11));验证漂移速度时不要用整段轨迹去拟合取后半段避开初始回旋调整区idx t 20*pi; y_window y(idx, 2); t_window t(idx); slope diff(y_window([1, end])) / diff(t_window([1, end])); disp(slope); % 理论值应为 0.5这个算例的判断标准是数值漂移速度与理论值偏差在 0.5% 以内。偏差偏大时优先检查初始速度里是否携带了额外的直流分量。E×B 漂移方向只取决于 E 和 B 的叉积方向与粒子电荷正负完全无关所以如果你把粒子换成负电荷后漂移方向反了那一定是场函数写错。4.2 变梯度磁场磁镜反射的临界条件把磁场改成随空间变化的梯度场就出现一个在聚变装置和空间物理里都常见的现象——磁镜。取对称磁场B(z) 1 (z/L)²其中 L 是磁场梯度尺度。当粒子从低场区向高场区运动时磁矩 μ v⊥²/(2B) 近似守恒平行速度 v∥ 会逐渐减小减到零时粒子被反射回去。用归一化参数写场函数L 5; B_fun (x, t) [0; 0; 1 (x(3)/L)^2]; E_fun (x, t) [0; 0; 0];取两种初始条件对比y0_a [0; 0; -8; 1; 0; 0.5]; % v_perp1, v_par0.5会被反射 y0_b [0; 0; -8; 1; 0; 1.5]; % v_perp1, v_par1.5会穿出 opts odeset(RelTol, 1e-9, AbsTol, 1e-11); [tA, yA] ode45((t,y) particle_motion(t,y,E_fun,B_fun), ... [0, 200], y0_a, opts); [tB, yB] ode45((t,y) particle_motion(t,y,E_fun,B_fun), ... [0, 200], y0_b, opts);画出两个算例的 z 坐标随时间变化figure; plot(tA, yA(:,3), DisplayName, 会被反射); hold on; plot(tB, yB(:,3), DisplayName, 会穿出); xlabel(t / \omega_c^{-1}); ylabel(z / r_L); legend;反射的临界条件是 v∥/v⊥ sqrt((B_max−B_min)/B_min)。当前算例中 B_min 1 (8/5)² 3.56B_max 在 z8 处与 B_min 相等整体是对称磁镜所以判断标准变成粒子能否越过 z8 的对称点。y0_a 的平行速度在到达边界前降到零并反向y0_b 则越过边界继续前进。这个算例里仿真时长取 200 个回旋周期而不是固定时间因为反射粒子的来回周期会因初始条件变化时间跨度短了会误判成粒子被“困住”。4.3 时变电场 稳恒磁场共振加热的数值观察第三个场景把电场改成时间振荡一项模拟射频波与回旋运动的共振。取回旋频率 ω_c 1电场频率等于 ω_cwc 1; E_fun (x,t) [cos(wc*t); 0; 0]; % 电场频率等于回旋频率 B_fun (x,t) [0; 0; 1]; y0 [0; 0; 0; 0; 1; 0]; [t, y] ode45((t,y) particle_motion(t,y,E_fun,B_fun), ... [0, 30*pi], y0, odeset(RelTol, 1e-9, AbsTol, 1e-11));对比组把频率调离共振点E_off (x,t) [cos(0.8*wc*t); 0; 0]; [t2, y2] ode45((t,y) particle_motion(t,y,E_off,B_fun), ... [0, 30*pi], y0, odeset(RelTol, 1e-9, AbsTol, 1e-11));看速度模随时间的变化v_res sqrt(sum(y(:,4:6).^2, 2)); v_off sqrt(sum(y2(:,4:6).^2, 2)); plot(t, v_res, t2, v_off); xlabel(t / \omega_c^{-1}); ylabel(|v| / v_0);共振情况下粒子每转一圈都在电场的同一相位获得能量速度幅值近似随时间线性增长失谐时速度出现周期性起伏但整体没有持续增长。需要特别说明的是这里的共振加热不考虑空间边界。实际等离子体或加速器中粒子回旋半径增大后会跑出有限区域所以在做长时间模拟时应该给场函数加空间包络比如E_fun (x,t) [cos(wc*t) * exp(-(x(1)^2x(2)^2)/400); 0; 0];不加这个约束的算例只适合展示短时间共振行为。5. 混合场算例收尾前加一道能量残留检查混合场里最隐蔽的问题不是轨迹画不出来而是画出来的轨迹看着光滑、实则能量误差已经积累到不可接受。洛伦兹力中的磁场部分不做功所以只要电场项为零或只在有限时间作用粒子动能应该在无电场的时间段保持严格守恒。ode45 是自适应步长积分器但它不保证约束能量守恒检查能量残留是确认步长和容差设置是否合理的快速手段。把这一段写成独立的验证脚本% verify_energy.m % 纯磁场场景判断容差设置是否合理 E_zero (x,t) [0;0;0]; B_test (x,t) [0;0;10.2*cos(0.5*x(3))]; opts odeset(RelTol, 1e-8, AbsTol, 1e-10); [tV, yV] ode45((t,y) particle_motion(t,y,E_zero,B_test), ... [0, 100*pi], [0;0;0;1;0;0.5], opts); Ek 0.5 * sum(yV(:,4:6).^2, 2); resid (Ek - Ek(1)) / Ek(1); semilogy(tV, abs(resid)); xlabel(t / \omega_c^{-1}); ylabel(relative energy error);把容差从 1e-6 扫到 1e-10观察能量残留的变化tolerance_list [1e-6, 1e-8, 1e-10]; energy_residual zeros(1, 3); for k 1:3 rt tolerance_list(k); opts odeset(RelTol, rt, AbsTol, rt*1e-2); [~, yk] ode45((t,y) particle_motion(t,y,E_zero,B_test), ... [0, 100*pi], [0;0;0;1;0;0.5], opts); Ek_k 0.5 * sum(yk(:,4:6).^2, 2); energy_residual(k) max(abs((Ek_k - Ek_k(1)) / Ek_k(1))); end disp(energy_residual);一般数值结果会显示三档容差对应的能量漂移大约在 1e-4、1e-8、1e-10 这个量级。如果某档容差下能量残留比预期高两个数量级以上那就是场函数在空间局部变化太快ode45 没有在这些区域主动加密步长这时候要检查 AbsTol 是否设得比 RelTol 相对过松。如果能量残留出现单调上升的斜坡而不是随机波动说明积分时间跨度内误差系统性累积时间范围要拆成多段重新积分或者改用 ode15s 这一类适用于刚性的求解器。只有能量残留检查通过之后前面三种混合场算例的轨迹判据才真正有说服力把这段脚本放在与 particle_motion.m 同一目录每次换场函数前先运行一次能省下不少排查轨迹异常的时间。本文还有配套的精品资源点击获取