简介这份docx文档围绕带电粒子在电磁场中的运动系统讲解如何用MATLAB进行仿真与分析适合物理专业学生、科研人员以及希望用数值计算理解电磁场问题的MATLAB初学者。文档先介绍MATLAB的基础操作、数据表示、矩阵运算与作图方法再结合牛顿第二定律、电场力、洛伦兹力等知识点详细推导带电粒子在复合场中的运动原理。其中涵盖质量较大微粒在复合场中的运动、带电粒子垂直射入正交电场磁场等典型场景并给出完整的仿真结果分析与参数影响讨论帮助读者通过改变初始速度、电荷量等参数直观观察粒子轨迹变化。资源包为单个docx文件大小约608KB内部按概述、MATLAB基础、实验原理与仿真结果、应用及总结等章节组织结构清晰便于查阅。目前已有二百六十六人学习无论用于课程设计、物理实验竞赛还是科研预研都能借助其中的仿真思路和代码框架快速上手让抽象电磁场理论变得可视化、可验证。1. 带电粒子在电磁场中的运动仿真的难点从来不是画螺旋线把“带电粒子在电磁场中运动的MATLAB仿真”布置成课程设计题目网上能搜到一批大同小异的作业代码。真正动手会发现让粒子转起来很容易难的是让轨迹转得可信回旋半径对不对、漂移方向对不对、长时间积分后会不会“仿真发散”。这个题目拆开有三层把洛伦兹力的二阶运动方程改写成求解器能接住的一阶常微分方程组选定一个在回旋频率尺度下不膨胀、不吞误差的积分方案再把三维轨迹里的漂移趋势用人眼一秒看穿。整套流程里涉及的数值积分、量纲检查和单位制换算往后做等离子体诊断、加速器束流模拟或者电磁场对带电尘埃轨迹的影响分析都会以类似面貌再次出现。适合动手的人群是物理、电子工程方向高年级学生以及刚接触MATLAB、想用仿真验证电磁场理论的工程师。有一个反直觉的注意点无电场时动能严格守恒而大多数写坏的仿真发散本质上是数值方法在悄悄给粒子注入能量。2. 洛伦兹力方程的数值解从二阶ODE到六维状态向量2.1 运动方程与六维状态MATLAB求解器要的是一阶形式带电粒子在外场中的运动由洛伦兹力给出m dv/dt q(E v×B)dr/dt v这是典型的二阶向量常微分方程而MATLAB的ode45、ode23t等求解器只认dy/dt f(t, y)的一阶形式。常见做法是把位置r和速度v拼成一个6维状态向量y [x, y, z, vx, vy, vz]^T再把原方程拆成两个一阶块前三个分量是速度后三个分量是加速度。这样一来求轨迹就是求这个6维向量随时间的变化。后面任何操作——plot3画线、检查动能、统计漂移速度——都从这6个分量里取数。求解器函数有一个必须注意的维度细节。ode45传入的状态向量是6×1列向量而磁场、电场通常是1×3行向量直接用cross(v, B)可能报维度不一致的错。下面这个函数写法在本科生的作业里最常见也最不容易错function dydt lorentz(t, y, q, m, E, B) r y(1:3).; % 提取位置分量转成1×3行向量 v y(4:6).; % 提取速度分量 a (q/m) * (E cross(v, B)); % 洛伦兹力给出的加速度 dydt [v, a].; % 拼接后转置保证输出是6×1列向量 end这段函数里有三个细节值得说明。第一形参t虽然在计算中没有用到但ode45回调函数签名要求有它不能省。第二r向量在这里没有参与运算但保留提取操作能让后面扩展位置相关场时少改一行代码。第三dydt必须显式转成列向量否则MATLAB会提示输出维度与初始状态不一致。这样写下来函数是自包含的E、B、q、m全部从外部传入换参数时不需要改动核心方程。2.2 数值积分方法与求解器选型手写欧拉法为什么会让半径膨胀很多第一次做这个仿真的人会选择自己写两层循环先由加速度更新速度再由速度更新位置。这就是前向欧拉法。在只有均匀磁场、没有电场的情况下粒子在每个时间步都沿圆弧切线方向走一小段速度方向转过一个角度但轨道半径会一点点向外“吹气球”。原因是欧拉法对旋转运动有一个内在的放大效应离散每一步都会让速度模长略微增长几千步之后轨迹半径就肉眼可见地变粗。这个现象常被误判成“仿真发散”其实只是积分器没有守住物理上的能量守恒。下表是几种常用离散方式的对比按“步长与回旋周期T2πm/(|q|B)的关系”来估算离散方式单步误差典型参数取法轨迹行为显式欧拉O(dt)dt T/20半径明显膨胀动能线性爬升显式欧拉O(dt)dt T/500短时间肉眼勉强可接受经典RK4O(dt^4)dt T/50轨迹与解析解偏差很小ode45自适应4/5阶默认容差即可常规均匀场下足够可靠这里的经验值是“做这个标题最常见”的做法与其纠结欧拉法步长不如直接用ode45。ode45在每个步长内部做两次不同阶数的估计用它们的差控制下一步长既能保证精度又不需要人为估计回旋周期的分点数。只有当E特别强、粒子在单个时间步内被剧烈加速或者回旋频率与其他运动尺度相差几个数量级时才需要考虑ode23t这类面向刚性问题的求解器。2.3 先算特征周期再定tspan仿真时长的标定顺序拿到一组物理参数第一件事不是写tspan[0,1]而是先算出两个特征量回旋周期T2πm/(|q|B)和回旋半径r_Lmv_perp/(|q|B)。例如一个质子质量m1.673×10^-27 kg电荷q1.602×10^-19 C磁感应强度B0.01 T它的回旋角频率约9.58×10^5 rad/s周期约6.55 μs。如果初速度的垂直分量是1×10^4 m/s回旋半径约0.0104 m也就是1.04 cm。看到这两组数量级就能立刻判断tspan该怎么设。仿真八到十个回旋周期时间范围应该是[0, 5.24×10^-5] s而不是[0,1]。如果tspan取到1秒相当于让粒子转了十五万圈ode45的输出点会挤满内存画出来的图也不是螺旋线而是一根实心“管子”。很多同学在这个步骤上栽跟头其实先算一次T所有问题都消失了。这个思想也贯穿后面的参数标定任何仿真都要先找到问题自带的特征时间尺度和特征长度尺度再决定时间窗口和解算精度。3. 用ode45在MATLAB里跑通第一个交叉电磁场轨迹3.1 场景设定E沿y、B沿z预期会出现E×B漂移最常见的教学场景是均匀磁场沿z轴正方向均匀电场沿y轴正方向粒子初速度沿x轴。这样的构型在等离子体物理里有明确解析解粒子在xy平面做回旋同时整体沿x方向漂移漂移速度为v_E E×B/B²。取E(0,100,0) V/m、B(0,0,0.01) T则v_E10000 m/s方向沿x。这个数字特意选成与初速度v0(10000,0,0)相同轨迹会是一条干净的摆线肉眼能直接验证漂移量。物理量取值说明电荷q1.602×10^-19 C质子改成负值可对比正负电荷质量m1.673×10^-27 kg质子质量磁场B00.01 T沿z电场E0100 V/m沿y初始位置(0,0,0)初始速度(1×10^4, 0, 0) m/s沿x回旋周期T6.55×10^-6 s由qB/m算出回旋半径1.04×10^-2 mv_perp/ω漂移速度1×10^4 m/sE0/B0这个表格本身就是一张“参数标定表”跑完仿真后逐项核对轨迹半径对不对总漂移距离是否约等于v_E乘以仿真时长。如果对不上不是代码问题就是物理参数写错。3.2 核心脚本matlab代码与参数说明下面是完整脚本直接复制到MATLAB编辑器保存运行即可% charged_particle_em.m % 带电粒子在交叉电磁场中的运动仿真SI单位制 clear; clc; close all; % ---------- 物理参数 ---------- q 1.602e-19; % 电荷量单位 C m 1.673e-27; % 质量单位 kg B0 0.01; % 磁感应强度单位 T E0 100; % 电场强度单位 V/m B [0, 0, B0]; % 磁场方向z E [0, E0, 0]; % 电场方向y % ---------- 初始条件 ---------- r0 [0, 0, 0]; % 初始位置 v0 [1e4, 0, 0]; % 初始速度 Y0 [r0, v0].; % 6维状态向量列向量 % ---------- 时间范围先算回旋周期 ---------- omega abs(q)*B0/m; % 回旋角频率 T_cyc 2*pi/omega; % 回旋周期 tspan [0, 8*T_cyc]; % 仿真8个回旋周期 % ---------- 求解 ---------- [t, Y] ode45((t,y) lorentz(t,y,q,m,E,B), tspan, Y0); % ---------- 从结果矩阵提取分量 ---------- x Y(:,1); yc Y(:,2); z Y(:,3); vx Y(:,4); vy Y(:,5); vz Y(:,6); % ---------- 三维轨迹图 ---------- figure(Color,w); plot3(x/1e-2, yc/1e-2, z/1e-2, b-, LineWidth, 1.2); xlabel(x (cm)); ylabel(y (cm)); zlabel(z (cm)); grid on; axis equal; title(质子E沿yB沿z轨迹为E×B漂移下的摆线); % ---------- 运动方程函数放在脚本末尾---------- function dydt lorentz(t, y, q, m, E, B) r y(1:3).; v y(4:6).; a (q/m) * (E cross(v, B)); dydt [v, a].; end几个参数的逻辑说明如下。tspan只指定了起止时间没有指定步长ode45会自动在内部调整步长返回的t向量是不等间隔的。轨迹的采样点足够密绘图时不会出现明显折角。plot3里的x/1e-2把单位从米换算成厘米是因为回旋半径只有1厘米量级直接画“米”会导致三条轴的数字都带着科学计数法横竖不直观。axis equal保证三个坐标轴比例一致不会因为窗口长宽比把圆形轨迹压成椭圆。如果运行后z坐标始终为0这是正常的因为初速度没有z分量且B垂直于xy平面粒子不会沿z方向运动想看螺旋线把初始速度改成v0[1e4, 0, 2e4]即可粒子会一边沿z轴匀速前进一边在xy平面回旋。3.3 在轨迹图上叠加E、B方向让漂移方向可解释纯轨迹图只能让读者看到“在动”但解释不了“为什么往这个方向漂”。在图上叠加两个箭头信息量立刻不同hold on; % 磁场方向用红色箭头长度只作示意 quiver3(0,0,0, 0,0,0.04, r, LineWidth, 1.5, MaxHeadSize, 2); % 电场方向用绿色箭头 quiver3(0,0,0, 0,0.04,0, g, LineWidth, 1.5, MaxHeadSize, 2); legend(轨迹, B (z), E (y), Location, best);箭头长度并非物理尺度的真实比例这里只是把“方向”画出来红色指z绿色指y。把这个方向图与漂移公式对照就能确认轨迹整体前进方向是x而不是y。这个习惯很重要真正做科研绘图时把场方向画进图里会让审稿人少问一个问题。MATLAB的quiver3的第三个三元组是箭头向量不需要刻意缩放但要保证起点在轨迹起点附近。4. 均匀场仿真里最容易把轨迹弄“飘”的三个设置4.1 单位制陷阱高斯单位制里漏掉光速结果直接飞出台面电磁学文献里经常混用两种单位制。SI制下运动方程是m dv/dt q(E v×B)而高斯单位制CGS写成m dv/dt q(E v×B/c)其中E和B的量纲也完全不同。如果从某本书上抄到CGS公式又直接用SI的“伏特/米”和“特斯拉”代入计算粒子会在几个时间步内获得比物理真实值大3×10^8倍的加速度直观表现就是轨迹瞬间冲出坐标范围被当成“仿真发散”。判断方法很简单关掉电场只保留均匀磁场跑一个纯回旋把模拟出的回旋半径与r_L m v_perp/(|q|B)对比。如果两者一致单位制没问题如果半径差了好几个数量级先检查公式里是不是混进了光速c。这套自检流程比盯着报错信息有效得多因为单位制错误通常不产生运行时报错只是结果离谱。4.2 步长与求解器报告“发散”先看动能曲线回旋频率特别高时定步长欧拉法很容易“飘”。这里给的不是具体步长阈值而是一条检查路径如果仿真结果出现半径增长、轨迹抖动或报错“步长在端点趋于零”先画出动能曲线。无电场情形下动能应当是一条水平直线如果你看到动能随时间线性上升说明积分器在向系统注入数值能量。改用下面这组更严格的容差设置能解决大部分非刚性问题opts odeset(RelTol, 1e-8, AbsTol, 1e-10); [t, Y] ode45((t,y) lorentz(t,y,q,m,E,B), tspan, Y0, opts);RelTol控制的是相对误差对轨迹形态影响最大AbsTol兜底防止某个分量绝对值接近零时误差被无限放大。注意这两个容差不要同时设得太苛刻否则ode45会为了满足精度把步长压得极小计算时间陡增。一般先用默认值跑通再逐级收紧RelTol观察轨迹是否变化若轨迹几乎不变就说明默认容差已够用。4.3 E×B漂移的方向检查叉积顺序错了整体会反着走漂移速度的公式是v_E E×B/B²注意叉积顺序不能反。以第3章场景为例B沿zE沿yE×B的结果是沿x所以粒子整体向x正方向漂移。如果跑出来的轨迹整体沿-x走第一个要怀疑的就是E或B某一个分量符号写反了。用一段差分代码能从数值结果里反推漂移速度% 取轨迹末尾段做差分避免初始回旋相位干扰 idx end-200:end; drift_sim mean(diff(x(idx)) ./ diff(t(idx))); drift_theory E0 / B0; fprintf(仿真漂移速度%.3e m/s\n, drift_sim); fprintf(理论漂移速度%.3e m/s\n, drift_theory);为什么要取末尾段而不是整条轨迹因为粒子初速度里叠加了回旋分量瞬时速度一直在振荡只有对足够多个周期做平均回旋部分才会消掉剩下纯漂移速度。这个平均时间的选取也反过来印证了“先算回旋周期再定窗口”的价值末尾选200个点跨越多少个周期取决于ode45输出的点数密度。5. 把验证写进可视化animatedline追踪粒子与动能守恒检查5.1 用animatedline实现粒子轨迹的动态演示静态plot3适合出论文插图但课程答辩或汇报里动态轨迹更能说明问题。MATLAB里最省事的做法是用animatedline先创建一根动画线再循环添加新点figure(Color,w); h animatedline(Color,b, LineWidth, 1.2, MaximumNumPoints, 300); xlabel(x (cm)); ylabel(y (cm)); zlabel(z (cm)); grid on; axis equal; % 每隔50个点添加一个数据点画得太密会拖慢帧率 for k 1:50:length(t) addpoints(h, x(k)/1e-2, yc(k)/1e-2, z(k)/1e-2); drawnow limitrate; endMaximumNumPoints设为300的含义是动画线只保留最近300个点粒子走远后轨迹尾巴会“断掉”从而清晰看出整体漂移趋势。如果想去掉这个限制、显示完整轨迹直接删掉这个参数即可。drawnow limitrate限制了刷新频率避免循环被绘图拖慢。这里的采样间隔50是经验值取决于ode45输出的总点数如果动画跳帧严重就把50调大。5.2 用动能曲线和相平面验证一套不需要解析解的校对方法在没有解析解的非均匀场场景里能量守恒依然是最通用的验证工具。把E0改成0重跑一次然后画动能KE 0.5*m*(vx.^2 vy.^2 vz.^2); % 动能单位 J figure; plot(t/1e-6, KE/1.6e-13, b-); xlabel(t (μs)); ylabel(动能 (keV));这里除以1.6×10^-13是把焦耳换算成keV方便读图。理想的曲线是一条水平直线如果曲线缓慢上翘继续收紧RelTol重新求解。另一条校验路径是画vx-vy相平面磁场均匀时相轨迹应当近似是一个圆圆心在(v_E, 0)半径由初速度与漂移速度的差决定。圆不闭合说明能量有损耗或增益圆整体移动说明存在持续的净加速这些信息能直接定位到错误来源。相平面检查还有一个附带好处它能直观展示“回旋运动”与“漂移运动”的分离。带电粒子在电磁场中运动的MATLAB仿真做到这一步已经不再只是交作业而是一套可以复用到非均匀场、时变场和带电粒子束传输问题的标准验证流程。下次换一组参数先画动能曲线再画相平面两步下来基本能确定仿真结果敢不敢拿出去用。本文还有配套的精品资源点击获取 带电粒子在电磁场中运动的MATLAB仿真：从洛伦兹力到可视化验证
简介这份docx文档围绕带电粒子在电磁场中的运动系统讲解如何用MATLAB进行仿真与分析适合物理专业学生、科研人员以及希望用数值计算理解电磁场问题的MATLAB初学者。文档先介绍MATLAB的基础操作、数据表示、矩阵运算与作图方法再结合牛顿第二定律、电场力、洛伦兹力等知识点详细推导带电粒子在复合场中的运动原理。其中涵盖质量较大微粒在复合场中的运动、带电粒子垂直射入正交电场磁场等典型场景并给出完整的仿真结果分析与参数影响讨论帮助读者通过改变初始速度、电荷量等参数直观观察粒子轨迹变化。资源包为单个docx文件大小约608KB内部按概述、MATLAB基础、实验原理与仿真结果、应用及总结等章节组织结构清晰便于查阅。目前已有二百六十六人学习无论用于课程设计、物理实验竞赛还是科研预研都能借助其中的仿真思路和代码框架快速上手让抽象电磁场理论变得可视化、可验证。1. 带电粒子在电磁场中的运动仿真的难点从来不是画螺旋线把“带电粒子在电磁场中运动的MATLAB仿真”布置成课程设计题目网上能搜到一批大同小异的作业代码。真正动手会发现让粒子转起来很容易难的是让轨迹转得可信回旋半径对不对、漂移方向对不对、长时间积分后会不会“仿真发散”。这个题目拆开有三层把洛伦兹力的二阶运动方程改写成求解器能接住的一阶常微分方程组选定一个在回旋频率尺度下不膨胀、不吞误差的积分方案再把三维轨迹里的漂移趋势用人眼一秒看穿。整套流程里涉及的数值积分、量纲检查和单位制换算往后做等离子体诊断、加速器束流模拟或者电磁场对带电尘埃轨迹的影响分析都会以类似面貌再次出现。适合动手的人群是物理、电子工程方向高年级学生以及刚接触MATLAB、想用仿真验证电磁场理论的工程师。有一个反直觉的注意点无电场时动能严格守恒而大多数写坏的仿真发散本质上是数值方法在悄悄给粒子注入能量。2. 洛伦兹力方程的数值解从二阶ODE到六维状态向量2.1 运动方程与六维状态MATLAB求解器要的是一阶形式带电粒子在外场中的运动由洛伦兹力给出m dv/dt q(E v×B)dr/dt v这是典型的二阶向量常微分方程而MATLAB的ode45、ode23t等求解器只认dy/dt f(t, y)的一阶形式。常见做法是把位置r和速度v拼成一个6维状态向量y [x, y, z, vx, vy, vz]^T再把原方程拆成两个一阶块前三个分量是速度后三个分量是加速度。这样一来求轨迹就是求这个6维向量随时间的变化。后面任何操作——plot3画线、检查动能、统计漂移速度——都从这6个分量里取数。求解器函数有一个必须注意的维度细节。ode45传入的状态向量是6×1列向量而磁场、电场通常是1×3行向量直接用cross(v, B)可能报维度不一致的错。下面这个函数写法在本科生的作业里最常见也最不容易错function dydt lorentz(t, y, q, m, E, B) r y(1:3).; % 提取位置分量转成1×3行向量 v y(4:6).; % 提取速度分量 a (q/m) * (E cross(v, B)); % 洛伦兹力给出的加速度 dydt [v, a].; % 拼接后转置保证输出是6×1列向量 end这段函数里有三个细节值得说明。第一形参t虽然在计算中没有用到但ode45回调函数签名要求有它不能省。第二r向量在这里没有参与运算但保留提取操作能让后面扩展位置相关场时少改一行代码。第三dydt必须显式转成列向量否则MATLAB会提示输出维度与初始状态不一致。这样写下来函数是自包含的E、B、q、m全部从外部传入换参数时不需要改动核心方程。2.2 数值积分方法与求解器选型手写欧拉法为什么会让半径膨胀很多第一次做这个仿真的人会选择自己写两层循环先由加速度更新速度再由速度更新位置。这就是前向欧拉法。在只有均匀磁场、没有电场的情况下粒子在每个时间步都沿圆弧切线方向走一小段速度方向转过一个角度但轨道半径会一点点向外“吹气球”。原因是欧拉法对旋转运动有一个内在的放大效应离散每一步都会让速度模长略微增长几千步之后轨迹半径就肉眼可见地变粗。这个现象常被误判成“仿真发散”其实只是积分器没有守住物理上的能量守恒。下表是几种常用离散方式的对比按“步长与回旋周期T2πm/(|q|B)的关系”来估算离散方式单步误差典型参数取法轨迹行为显式欧拉O(dt)dt T/20半径明显膨胀动能线性爬升显式欧拉O(dt)dt T/500短时间肉眼勉强可接受经典RK4O(dt^4)dt T/50轨迹与解析解偏差很小ode45自适应4/5阶默认容差即可常规均匀场下足够可靠这里的经验值是“做这个标题最常见”的做法与其纠结欧拉法步长不如直接用ode45。ode45在每个步长内部做两次不同阶数的估计用它们的差控制下一步长既能保证精度又不需要人为估计回旋周期的分点数。只有当E特别强、粒子在单个时间步内被剧烈加速或者回旋频率与其他运动尺度相差几个数量级时才需要考虑ode23t这类面向刚性问题的求解器。2.3 先算特征周期再定tspan仿真时长的标定顺序拿到一组物理参数第一件事不是写tspan[0,1]而是先算出两个特征量回旋周期T2πm/(|q|B)和回旋半径r_Lmv_perp/(|q|B)。例如一个质子质量m1.673×10^-27 kg电荷q1.602×10^-19 C磁感应强度B0.01 T它的回旋角频率约9.58×10^5 rad/s周期约6.55 μs。如果初速度的垂直分量是1×10^4 m/s回旋半径约0.0104 m也就是1.04 cm。看到这两组数量级就能立刻判断tspan该怎么设。仿真八到十个回旋周期时间范围应该是[0, 5.24×10^-5] s而不是[0,1]。如果tspan取到1秒相当于让粒子转了十五万圈ode45的输出点会挤满内存画出来的图也不是螺旋线而是一根实心“管子”。很多同学在这个步骤上栽跟头其实先算一次T所有问题都消失了。这个思想也贯穿后面的参数标定任何仿真都要先找到问题自带的特征时间尺度和特征长度尺度再决定时间窗口和解算精度。3. 用ode45在MATLAB里跑通第一个交叉电磁场轨迹3.1 场景设定E沿y、B沿z预期会出现E×B漂移最常见的教学场景是均匀磁场沿z轴正方向均匀电场沿y轴正方向粒子初速度沿x轴。这样的构型在等离子体物理里有明确解析解粒子在xy平面做回旋同时整体沿x方向漂移漂移速度为v_E E×B/B²。取E(0,100,0) V/m、B(0,0,0.01) T则v_E10000 m/s方向沿x。这个数字特意选成与初速度v0(10000,0,0)相同轨迹会是一条干净的摆线肉眼能直接验证漂移量。物理量取值说明电荷q1.602×10^-19 C质子改成负值可对比正负电荷质量m1.673×10^-27 kg质子质量磁场B00.01 T沿z电场E0100 V/m沿y初始位置(0,0,0)初始速度(1×10^4, 0, 0) m/s沿x回旋周期T6.55×10^-6 s由qB/m算出回旋半径1.04×10^-2 mv_perp/ω漂移速度1×10^4 m/sE0/B0这个表格本身就是一张“参数标定表”跑完仿真后逐项核对轨迹半径对不对总漂移距离是否约等于v_E乘以仿真时长。如果对不上不是代码问题就是物理参数写错。3.2 核心脚本matlab代码与参数说明下面是完整脚本直接复制到MATLAB编辑器保存运行即可% charged_particle_em.m % 带电粒子在交叉电磁场中的运动仿真SI单位制 clear; clc; close all; % ---------- 物理参数 ---------- q 1.602e-19; % 电荷量单位 C m 1.673e-27; % 质量单位 kg B0 0.01; % 磁感应强度单位 T E0 100; % 电场强度单位 V/m B [0, 0, B0]; % 磁场方向z E [0, E0, 0]; % 电场方向y % ---------- 初始条件 ---------- r0 [0, 0, 0]; % 初始位置 v0 [1e4, 0, 0]; % 初始速度 Y0 [r0, v0].; % 6维状态向量列向量 % ---------- 时间范围先算回旋周期 ---------- omega abs(q)*B0/m; % 回旋角频率 T_cyc 2*pi/omega; % 回旋周期 tspan [0, 8*T_cyc]; % 仿真8个回旋周期 % ---------- 求解 ---------- [t, Y] ode45((t,y) lorentz(t,y,q,m,E,B), tspan, Y0); % ---------- 从结果矩阵提取分量 ---------- x Y(:,1); yc Y(:,2); z Y(:,3); vx Y(:,4); vy Y(:,5); vz Y(:,6); % ---------- 三维轨迹图 ---------- figure(Color,w); plot3(x/1e-2, yc/1e-2, z/1e-2, b-, LineWidth, 1.2); xlabel(x (cm)); ylabel(y (cm)); zlabel(z (cm)); grid on; axis equal; title(质子E沿yB沿z轨迹为E×B漂移下的摆线); % ---------- 运动方程函数放在脚本末尾---------- function dydt lorentz(t, y, q, m, E, B) r y(1:3).; v y(4:6).; a (q/m) * (E cross(v, B)); dydt [v, a].; end几个参数的逻辑说明如下。tspan只指定了起止时间没有指定步长ode45会自动在内部调整步长返回的t向量是不等间隔的。轨迹的采样点足够密绘图时不会出现明显折角。plot3里的x/1e-2把单位从米换算成厘米是因为回旋半径只有1厘米量级直接画“米”会导致三条轴的数字都带着科学计数法横竖不直观。axis equal保证三个坐标轴比例一致不会因为窗口长宽比把圆形轨迹压成椭圆。如果运行后z坐标始终为0这是正常的因为初速度没有z分量且B垂直于xy平面粒子不会沿z方向运动想看螺旋线把初始速度改成v0[1e4, 0, 2e4]即可粒子会一边沿z轴匀速前进一边在xy平面回旋。3.3 在轨迹图上叠加E、B方向让漂移方向可解释纯轨迹图只能让读者看到“在动”但解释不了“为什么往这个方向漂”。在图上叠加两个箭头信息量立刻不同hold on; % 磁场方向用红色箭头长度只作示意 quiver3(0,0,0, 0,0,0.04, r, LineWidth, 1.5, MaxHeadSize, 2); % 电场方向用绿色箭头 quiver3(0,0,0, 0,0.04,0, g, LineWidth, 1.5, MaxHeadSize, 2); legend(轨迹, B (z), E (y), Location, best);箭头长度并非物理尺度的真实比例这里只是把“方向”画出来红色指z绿色指y。把这个方向图与漂移公式对照就能确认轨迹整体前进方向是x而不是y。这个习惯很重要真正做科研绘图时把场方向画进图里会让审稿人少问一个问题。MATLAB的quiver3的第三个三元组是箭头向量不需要刻意缩放但要保证起点在轨迹起点附近。4. 均匀场仿真里最容易把轨迹弄“飘”的三个设置4.1 单位制陷阱高斯单位制里漏掉光速结果直接飞出台面电磁学文献里经常混用两种单位制。SI制下运动方程是m dv/dt q(E v×B)而高斯单位制CGS写成m dv/dt q(E v×B/c)其中E和B的量纲也完全不同。如果从某本书上抄到CGS公式又直接用SI的“伏特/米”和“特斯拉”代入计算粒子会在几个时间步内获得比物理真实值大3×10^8倍的加速度直观表现就是轨迹瞬间冲出坐标范围被当成“仿真发散”。判断方法很简单关掉电场只保留均匀磁场跑一个纯回旋把模拟出的回旋半径与r_L m v_perp/(|q|B)对比。如果两者一致单位制没问题如果半径差了好几个数量级先检查公式里是不是混进了光速c。这套自检流程比盯着报错信息有效得多因为单位制错误通常不产生运行时报错只是结果离谱。4.2 步长与求解器报告“发散”先看动能曲线回旋频率特别高时定步长欧拉法很容易“飘”。这里给的不是具体步长阈值而是一条检查路径如果仿真结果出现半径增长、轨迹抖动或报错“步长在端点趋于零”先画出动能曲线。无电场情形下动能应当是一条水平直线如果你看到动能随时间线性上升说明积分器在向系统注入数值能量。改用下面这组更严格的容差设置能解决大部分非刚性问题opts odeset(RelTol, 1e-8, AbsTol, 1e-10); [t, Y] ode45((t,y) lorentz(t,y,q,m,E,B), tspan, Y0, opts);RelTol控制的是相对误差对轨迹形态影响最大AbsTol兜底防止某个分量绝对值接近零时误差被无限放大。注意这两个容差不要同时设得太苛刻否则ode45会为了满足精度把步长压得极小计算时间陡增。一般先用默认值跑通再逐级收紧RelTol观察轨迹是否变化若轨迹几乎不变就说明默认容差已够用。4.3 E×B漂移的方向检查叉积顺序错了整体会反着走漂移速度的公式是v_E E×B/B²注意叉积顺序不能反。以第3章场景为例B沿zE沿yE×B的结果是沿x所以粒子整体向x正方向漂移。如果跑出来的轨迹整体沿-x走第一个要怀疑的就是E或B某一个分量符号写反了。用一段差分代码能从数值结果里反推漂移速度% 取轨迹末尾段做差分避免初始回旋相位干扰 idx end-200:end; drift_sim mean(diff(x(idx)) ./ diff(t(idx))); drift_theory E0 / B0; fprintf(仿真漂移速度%.3e m/s\n, drift_sim); fprintf(理论漂移速度%.3e m/s\n, drift_theory);为什么要取末尾段而不是整条轨迹因为粒子初速度里叠加了回旋分量瞬时速度一直在振荡只有对足够多个周期做平均回旋部分才会消掉剩下纯漂移速度。这个平均时间的选取也反过来印证了“先算回旋周期再定窗口”的价值末尾选200个点跨越多少个周期取决于ode45输出的点数密度。5. 把验证写进可视化animatedline追踪粒子与动能守恒检查5.1 用animatedline实现粒子轨迹的动态演示静态plot3适合出论文插图但课程答辩或汇报里动态轨迹更能说明问题。MATLAB里最省事的做法是用animatedline先创建一根动画线再循环添加新点figure(Color,w); h animatedline(Color,b, LineWidth, 1.2, MaximumNumPoints, 300); xlabel(x (cm)); ylabel(y (cm)); zlabel(z (cm)); grid on; axis equal; % 每隔50个点添加一个数据点画得太密会拖慢帧率 for k 1:50:length(t) addpoints(h, x(k)/1e-2, yc(k)/1e-2, z(k)/1e-2); drawnow limitrate; endMaximumNumPoints设为300的含义是动画线只保留最近300个点粒子走远后轨迹尾巴会“断掉”从而清晰看出整体漂移趋势。如果想去掉这个限制、显示完整轨迹直接删掉这个参数即可。drawnow limitrate限制了刷新频率避免循环被绘图拖慢。这里的采样间隔50是经验值取决于ode45输出的总点数如果动画跳帧严重就把50调大。5.2 用动能曲线和相平面验证一套不需要解析解的校对方法在没有解析解的非均匀场场景里能量守恒依然是最通用的验证工具。把E0改成0重跑一次然后画动能KE 0.5*m*(vx.^2 vy.^2 vz.^2); % 动能单位 J figure; plot(t/1e-6, KE/1.6e-13, b-); xlabel(t (μs)); ylabel(动能 (keV));这里除以1.6×10^-13是把焦耳换算成keV方便读图。理想的曲线是一条水平直线如果曲线缓慢上翘继续收紧RelTol重新求解。另一条校验路径是画vx-vy相平面磁场均匀时相轨迹应当近似是一个圆圆心在(v_E, 0)半径由初速度与漂移速度的差决定。圆不闭合说明能量有损耗或增益圆整体移动说明存在持续的净加速这些信息能直接定位到错误来源。相平面检查还有一个附带好处它能直观展示“回旋运动”与“漂移运动”的分离。带电粒子在电磁场中运动的MATLAB仿真做到这一步已经不再只是交作业而是一套可以复用到非均匀场、时变场和带电粒子束传输问题的标准验证流程。下次换一组参数先画动能曲线再画相平面两步下来基本能确定仿真结果敢不敢拿出去用。本文还有配套的精品资源点击获取