带电粒子在电磁场中运动的MATLAB仿真

目 录

第1章 概述......................................................... 1 第2章MATLAB的基础知识............................................. 2

2.1 MATLAB使用介绍.............................................. 2 2.2 MATLAB的基本知识............................................ 2 2.3利用MATLAB作图.............................................. 5 第3章 实验原理及仿真结果分析....................................... 7

3.1 带电粒子在电磁场中运动的原理 ................................ 7 3.2 质量较大的带电微粒在复合场中的运动 .......................... 7 3.3 带电粒子垂直射入E和B正交的叠加场的运动分析 ................ 8 3.4 实验内容 .................................................... 8 第4章 MATLAB仿真的应用 ........................................... 11

4.1 用MATLAB仿真带电粒子在磁场中运动的优点 .................... 11 4.2 用MATLAB仿真在物理实验中的应用 ............................ 11 第5章 总结........................................................ 12 参考文献........................................................... 13 附 录.............................................................. 14

第1章 概述

计算机数值模拟的研究方法已成为继实验研究和理论分析之外的第三种研究手段，在基础物理学习中适当引入计算机数值方法，有助于将一些高深的物理知识深入浅出、生动形象地学习。随着计算机的普及，MATLAB在基础物理中的应用日益广泛。MATLAB是当今最优秀的科技应用软件之一，它以强大的科学计算与可视化功能、简单易用、开放式可扩展环境，特别是所附带的30 多种面向不同领域的工具箱支持，使得它在许多科学领域中成为计算机辅助设计和分析、算法研究和应用开发的基本工具和首选平台。MATLAB具有其他高级语言难以比拟的一些优点，如编写简单、编程效率高、易学易懂等，因此MATLAB 语言也被通俗地称为演算纸式科学算法语言。在控制、通信、信号处理及科学计算等领域中，MATLAB 都被广泛地应用，已经被认为能够有效提高工作效率、改善设计手段的工具软件，掌握了MATLAB 好比掌握了开启这些专业领域大门的钥匙[1]。

带电体在复合场中运动的基本分析:这里所讲的复合场指电场、磁场和重力场并存, 或其中某两场并存, 或分区域存在, 带电体连续运动时, 一般须同时考虑电场力、洛仑兹力和重力的作用。

在不计粒子所受的重力的情况下，带电粒子只受电场和洛仑兹力的作用，粒子所受的合外力就是这两种力的合力，其运动加速度遵从牛顿第二定律。在相互垂直的匀强电场与匀强磁场构成的复合场中，如果粒子所受的电场力与洛仑兹力平衡，粒子将做匀速直线运动；如果所受的电场力与洛仑兹力不平衡，粒子将做一般曲线运动，而不可能做匀速圆周运动，也不可能做与抛体运动类似的运动。在相互垂直的点电荷产生的平面电场与匀强磁场垂直的复合场中，带电粒子有可能绕场电荷做匀速圆周运动。

1

第2章MATLAB的基础知识

2.1 MATLAB使用介绍

在Windows窗口中用鼠标双击Matlab图标即可进入Matlab的工作窗口（Command Window），如图2-11所示。没有图标可利用Matlab\\bin目录下的Matlab.exe文件在桌面上建立一个快捷方式。

退出Matlab的方法有三种：单击工作窗口右上角的关闭按钮；用菜单File→Exit Matlab命令；或者直接在工作窗口中输入quit后回车。

工作窗口是标准的Windows窗口形式，用户在命令窗口中输入各种指令，进行运算；在左侧的变量窗口中监控当前所创立的所有变量[2]。

当前工作路径 命令窗口 变量窗口

图2-1 Matlab工作窗口

2.2 MATLAB的基本知识

2.2.1 MATLAB基本命令

MATLAB基本命令如下表所示

2

表2-1 MATLAB基本命令

主题词 format who whos clear save load linspace length size max min sum 含义 设置数据显示格式 显示变量名 显示变量信息 清除内存变量 保存工作变量到文件 从文件装载变量 区间等分 获取数组长度 矩阵大小 最大值 最小值 求和 主题词 feval input disp tic toc help lookfor type which double str2num num2str 含义 函数求值 提示输入 输出 启动秒表 时间读数（秒） 帮助 查找 列程序清单 查找文件目录 双精度 字符串转化为数值 数值转化为字符串 2.2.2 数据和变量

1.表达式

在命令窗口作一些简单的计算，就如同使用一个功能强大的计算器，使用变量无须预先定义类型。

几个表达式可以写在一行，用分号（；）或逗号（，）分割，用分号（；）使该表达式运算结果不显示，而逗号（，）则显示结果。也可以将一个长表达式分在几行上写，用三点（…）续行。

若需要修改已执行过的命令行，可以在命令历史中找到该命令行复制，再粘贴至命令窗口修改。也可以直接使用键盘↑↓调出已执行过的命令行修改。

2.数据显示格式

MATLAB默认的数据显示格式为短格式（short）：当结果为整数，就作为整数显示；当结果是实数，以小数点后四位的长度显示。若结果的有效数字超出一定范围，以科学计数法显示（如3.2000e-006表示3.2106）。数据显示格式可使用命令Format改变。

3.预定义变量

MATLAB有一些预定义变量（表2-2），启动时就已赋值，可以直接使用，如前我们使用的圆周率pi和虚数单位i。

表2-2 常用预定义变量

变量名 i或j

说 明 虚数单位1 3

pi eps realmin realmax inf NaN 圆周率3.14159… 浮点数识别精度2^(-52)=2.22041016 最小正实数2.225110308 最大正实数1.797710308 无穷大 没有意义的数 4.用户变量

MATLAB变量名总以字母开头，以字母、数字或下划线组成，区分大小写，有效字符长度为63个。如A，a，a1，a_b都是合法的，且a与A表示不同变量。在Command Window中使用的变量一旦被赋值，就会携带这个值存在于工作空间，直到被清除或被赋予新的值。

2.2.3数组和矩阵运算

MATLAB基本数据单元是无需指定维数的数组。数组运算是MATLAB最鲜明的特点，一方面可以使得计算程序简明易读，另一方面可以提高计算速度。

1.数组的输入

最常用的数组是双精度数值数组（double array）。一维数组相当于向量，二维数组相当于矩阵，一维数组可以视为二维数组的特例。二维数组的第一维称为“行”，第二维称为“列”。MATLAB数组无需预先定义维数。直接输入数组的元素，用中括号（[]）表示一个数组，同行元素间用空格或逗号分隔，不同行间用分号或回车分隔。

2.数组运算

数组运算是指数组对应元素之间的运算，也称点运算。矩阵的乘法、乘方和除法有特殊的数学含义，并不是数组对应元素的运算，所数组乘法、乘方和除法的运算符前特别加了一个点。特别要区分数组运算在乘法、乘方和除法上的意义和表示上与矩阵运算的不同[3]。

表2-3 数组运算符

运算 数组加与减 数乘数组 数与数组加减 数组乘数组 数组乘方 数除以数组 数组除法 符号 A+B与A-B k*A或A*k k+A或k-A A.*B A.^k,k.^A k./A 左除A.\\B,右除B./A 说明 对应元素之间加减 k乘A的每个元素 k加（减）A的每个元素 点运算只有点乘、点乘方、点除三个，表示对应元素之间的运算；（.*）是一个整体，点(.)不能漏掉，（.）和（*）之间也不能有空格 4

3.数学函数

数组的数学函数也是按每个元素的运算，使用通常的函数符号，常用数学函数见表2-4。

表2-4 数学函数

函数 sin cos tan cot asin acos sqrt 意义 正弦 余弦 正切 余切 反正弦 反余弦 开方 函数 fix mod rem abs exp log log10 意义 向0取整 模余 除法余数 绝对值（模） 指数函数 自然对数 以10为底的对数 2.3利用MATLAB作图

1.常用作图命令和函数（表2-5）。

表2-5 常用作图命令和函数

主题词 plot fplot ezplot plot3 meshgrid mesh surf contour contour3 title xlabel 含 义 基本二维图形 一元函数图像 画二维曲线的符号命令 空间曲线 网格数据生成 网面图 曲面图 等高线图 三维等高线图 标题 x轴说明 主题词 clabel grid hold axis view subplot figure clf close ylabel zlabel 含 义 等高线高度标志 格栅 图形保持 定制坐标轴 改变视点 子图 新图形窗口 清除图形 关闭图形窗口 y轴说明 z轴说明 2.曲线图

plot(x,y) 作出以数据(x(i),y(i))为节点的折线图，其中x,y为同长度的向量 fplot('fun',[a,b]) 作出函数fun在区间[a,b]上的函数图，fun可以是M函数主名，

也可是字符串

ezplot(fun,lims) 绘制字符串fun（可是显函数、隐函数或参数方程）指定的函数 plot3(x,y,z) 空间曲线图，其中x，y，z为同长度向量

图形显示在图形窗口。在图形窗口可以使用File菜单保存(Save)为M文件，

5

导出(Export)为图形文件。也可利用图形窗口Edit菜单Copy figure作为图片复制到剪贴板，从而进一步粘贴到Word或其他应用程序中。

图形的线型、标记、颜色均可根据要求设定。常用的见表2-6。

表2-6 图形元素设定

颜色 b g r m y k 蓝（默认） 绿 红 洋红 青 黑 - ： -. -- 线型 实线（默认） 虚线 点划线 划线 . o x + * 标记 无标记（默认） 点 圈 叉 十字 星 3.曲面图

[x,y]=meshgrid(xa,ya) 当xa,ya分别为m维和n维行向量，得到x和y均为n

行m列矩阵。meshgrid常用于生成X-Y平面上的网格数据

mesh(x,y,z) 绘制网面图，是最基本的曲面图形命令，其中x,y,z是

同阶矩阵，表示曲面三维数据

surf(x,y,z) 绘制曲面图，与mesh用法类似 contour(x,y,z) 绘制等高线图，与mesh用法类似 contour3(x,y,z) 绘制三维等高线图，与mesh用法类似

4.图形说明和定制

title('字符串') 图形标题说明

xlabel,ylabel,zlabel 用法类似于title，分别说明坐标轴x,y,z grid on/off 显示/不显示格栅 box on/off 使用/不使用坐标框

axis([a,b,c,d]) 定制二维坐标轴范围 alegend(str1,str2,……) 图例，字符串str1,str2,……依次为各图形对象说明

6

第3章 实验原理及仿真结果分析

3.1 带电粒子在电磁场中运动的原理

带电粒子在磁场中运动时会受到洛伦兹力的作用，且随着初始运动方向和磁

场分布的不同，其运动轨迹会发生不同的变化。由洛伦兹力的推导公式可知，它垂直于粒子的运动速度，不对运动粒子作功，只改变其运动方向，其大小为：

FqBsin (3—1)

因此，综合牛顿运动定律就可以精确确定带电粒子在磁场中的运动轨迹。

3.2 质量较大的带电微粒在复合场中的运动

这里我们只研究垂直射入磁场的带电微粒在垂直磁场的平面内的运动，并分几种情况进行讨论。

（1）只受重力和洛仑兹力：此种情况下，要使微粒在垂直磁场的平面内运动，磁场方向必须是水平的。微粒所受的合外力就是重力与洛仑兹力的合力。在此合力作用下，微粒不可能再做匀速圆周运动，也不可能做与抛体运动类似的运动。在合外力不等于零的情况下微粒将做一般曲线运动，其运动加速度遵从牛顿第二定律；在合外力等于零的情况下，微粒将做匀速直线运动。

无论微粒在垂直匀强磁场的平面内如何运动，由于洛仑兹力不做功，只有重力做功，因此微粒的机械能守恒，即：

1212 （3—2） mghamvamghbmvb22 （2）微粒受有重力、电场力和洛仑兹力：此种情况下。要使微粒在垂直磁场的平面内运动，匀强磁场若沿水平方向，则所加的匀强电场必须与磁场方向垂直。

在上述复合场中，带电微粒受重力、电场力和洛仑兹力。这三种力的矢量和即是微粒所受的合外力，其运动加速度遵从牛顿第二定律。如果微粒所受的重力与电场力相抵消，微粒相当于只受洛仑兹力，微粒将以洛仑兹力为向心力，以射入时的速率做匀速圆周运动。若重力与电场力不相抵，微粒不可能再做匀速圆周运动，也不可能做与抛体运动类似的运动，而只能做一般曲线运动。如果微粒所受的合外力为零，即所受的三种力平衡，微粒将做匀速直线运动。

无论微粒在复合场中如何运动，洛仑兹力对微粒不做功。若只有重力对微粒做功，则微粒的机械能守恒；若只有电场力对微粒做功，则微粒的电势能和动能

7

的总和守恒；若重力和电场力都对微粒做功，则微粒的电势能与机械能的总和守恒，用公式表示为：

1212 （3—3） qUamghamvaqUbmghbmvb223.3 带电粒子垂直射入E和B正交的叠加场的运动分析

带电粒子垂直射入E和B正交的叠加场常见的模型就是速度旋转器。如图3-1（a）所示，粒子受力特点——电场力F与洛仑兹力f方向相反。

粒子匀速通过速度选择器的条件——带电粒子从小孔S1水平射入，匀速通过叠加场，并从小孔S2水平射出，从不同角度看有三种等效条件：从力的角度—电场力与洛仑兹力平衡，即：qEBq0；从速度角度——v0的大小等于E与B的比值，即：v0E；从功的角度——电场力对粒子不做功，即：WFqEd0； BE，粒子将因侧BE时：粒子向fB两个重要的功能关系——当粒子进入速度选择器时速度v0移而不能通过选择器。

如图3-1（b），设在电场方向侧移d后粒子速度为v，当v方向侧移，F做负功—粒子动能减少电势能增加，有：

121mv0qEdmv2； 22E当<时，粒子向F方向侧移，F做正功——粒子动能增加，电势能减少，

B112有： mv0qEdmv2

22

(a) （b）

图3-1 速度选择器原理图

3.4 实验内容

1、1．用Matlab 数值模拟的方法模拟带电粒子在恒定磁场中的螺旋运动，即

8

带电粒子进入磁场的方向与磁场方向的角度（0 90）。

（1）带电粒子质量为m, 电量为q( q>0)的带电粒子在均匀稳定的电磁场中的 运动微分方程为：

md2rqEqB （3-4）

以电磁场中某点为原点，以E为Oy方向，B为Oz方向建立坐标系O-xyz。由于 qBm，则上式的投影方程为：

d2xdy2dtdtd2yqEdx （3-5） 2mdtdtd2z20dt将其转换为一阶微分线性方程组，以便用差分替代微分作数值计算，令

w1x，w2dydx，w3y，w4dydt，w5z，w6dzdt，则上式化简为：

dw1dtdw2dtdw3dtdw4dtdw5dtdw6dtw2w4w4qEw2m （3-6）

w60下面使用MATLAB R2009b 编写程序，分三种情况考虑: ① 电场强度和磁场强度都不为零； ② 电场强度为零, 磁场强度不为零； ③ 电场强度不为零, 磁场强度为零。 程序见附录1。

运行程序得到该带电粒子在均匀稳定的电磁场中的运动轨迹如图3-2所示：

9

图 3-2 带电粒子在均匀稳定的电磁场中的运动轨迹

sin（即x0，（2）当v与B有一定夹角时，可将v分解为：cos，y0，z0）。若只有分量，带电粒子将在垂直的平面内作匀速圆周运动，若只有分量，粒子将沿磁场B方向作匀速直线运动，当两分量同时存在时，带

电粒子的运动轨迹将是一条螺旋线。用MATLAB R2009b编写程序，源程序见附录2.

运行程序，得到带电粒子的运动轨迹如图3-3所示：

图 3-3 带电粒子的运动轨迹（螺旋线）

10

第4章 MATLAB仿真的应用

4.1 用MATLAB仿真带电粒子在磁场中运动的优点

由于电磁场看不见，摸不着，它不像普通的“三态”物质那样由原子、分子 构成，也没有可见的形态，但其具有可以被检测的运动速度、能量和动量，占有 空间，是一种真实的客观存在。实验中通过MATLAB 数值模拟模拟带电粒子在恒定磁场中的螺旋运动和磁聚焦现象，让我们对带电粒子在磁场中的运动有了更深的感性认识，用MATLAB 画出的立体图也更有利于对此的理解，对于对应知识的理解和吸收有很大的帮助。运用数值模拟方法作出的螺旋线运动和磁聚焦现象的轨迹不仅生动形象,而且结果是建立在定量计算的基础上的, 不是简单的定性的模拟, 因此更具有科学性与说服力。

4.2 用MATLAB仿真在物理实验中的应用

MATLAB现已成为集数值分析、符号运算、图形处理、系统仿真等功能于一体的科学与工程计算软件平台。MATLAB的最大特点在于功能强大而界面友好,其丰富的库函数和各种专用工具箱,将使用者从繁琐的底层编程中解放出来,使他们有更多的时间和精力去探究科学问题本身；它对科学计算结果迅捷而准确的可视化能力,有助于使用者化抽象思维为形象思维,从而更好地洞察含义、理解概念、发现规律。目前,MATLAB已成为工程界最流行的软件工具,在大学理工实验课程教学中的应用亦渐成热点[4]。

在构建物理模型的前提下,利用MATLAB的计算、声音、图像、动画等功能来模拟物理现象,诠释物理规律,演示运动过程,并由此构建学生自主探索问题和解决问题的平台。利用它开展大学物理计算机模拟实验,一方面使抽象的问题形象化,加深学生对物理现象和规律的理解,提高教学效率,另一方面则有利于激发学的学习兴趣,发展学生的想象力,开发学生的创造力[5]。在物理实验教学过程中,教师可利用MATLAB 绘图和动画功能,选择其中合适的部分用在自己的讲解和演示文稿及多媒体课件中,以便说明讲解的结构,形象地演示物理教学中某些难以理解的内容;或用模拟的图表、动画等展示动态的变化过程和理论模型等,帮助学生了解和理解所学的知识。形象而生动,非常便于学生接受。在信息技术、计算机辅助教学迅速发展的时代,上述工作显示了比较大的发展空间和潜力,已经广泛地引起了学生的各方面兴趣,正在成为学生自我学习的平台,已显示出在提高学生素质、提高学生学习质量以及在教学方法和教学内容改革等方面所具有的优势。

11

第5章 总结

基于MATLAB计算机物理模拟方法实现的物理图像和物理过程,教学意图明确,操作简单,人机交互性强,交互参数调节方便,响应速度快,具有很高的容错能力。在物理实验教学中,以MATLAB为实验演示平台,对物理实验现象进行模拟,作为演示配合物理理论的讲授,能使学生加深对物理现象和规律的理解,提高教学效率。MATLAB绘图和动画功能,能将计算结果或外部文件数据以二维、三维图形呈现而且提供了方便的图形导出功能,其图形可以jpg、tif 等多种格式导出,在物理实验教学中开发利用,可作为物理实验CAI 课件的素材图[6]。

通过这次试验，进一步学习和掌握了MATLAB的使用，为今后的学习和应用MATLAB打下了更好的基础。另外通过这次使用MATLAB仿真带电粒子在不同条件下在磁场中的运动，深入且形象的认识了带电粒子在磁场中的运动方式，而不在限于之前的理论认识。

同时，在实验时也遇到了不少的问题。1、还是像上次实验一样不能熟练使用MATLAB软件，导致花费时间较多。2、模拟效果感觉不是很好。希望通过这次的学习及今后的进一步学习使用MATLAB，能更好的使用MATLAB解决更多的问题。

12

参考文献

[1] Steven E. Koonin，秦克诚译.计算物理学.北京：高等教育出版社，1993. [2] 马文淦等。计算物理学.合肥：中国科学技术大学出版社，1992. [3] 张志涌。精通MATLAB6.5.北京：北京航空航天大学出版社，2003. [4] 罗志恒,陆建隆.MATLAB 在普通物理教学中的应用[J].物理通报,2003. [5] 梁齐,许晓琳、王明虎.光干涉实验的计算机模拟[J].大学物理验,2003. [6] 张星辉,钟晓春,张世昌.大学物理教学中引用Matlab图像和动画的实践与思

考[J].贵州大报,2004.

13

附 录

附录1

%用Matlab数值模拟的方法模拟带电粒子在恒定磁场中的螺旋运动

function lx

global q m b e q=1.6e-2; m=0.02; b=[2;1;0]; e=[1;0;1]; figure

strd{1}='e\\neq 0,b\\neq 0'; strd{2}='e=0,b\\neq 0'; strd{3}='e\\neq 0, b=0'; for i=1:3

[t,w]=ode23(@ddlzfun,[0:0.01:20],[0,0.01,0,6,0,0.01],[ ],q,m,b(i),e(i));

axes ('unit','normalized','position', [0.045+(i-1)*0.35 0.062 0.2786 0.6583 ]);

plot3(w(:,1),w(:,3),w(:,5),'linewidth',2); grid on

title(strd{i},'fontsize',12,'fontweight','demi'); xlabel('x'); ylabel('y'); zlabel('z'); view([-51,18]); end

function wdot=ddlzfun(t,w,q,m,b,e) %该函数实现(3)式的差分迭代 syms q m b e q=1.6e-2; m=0.02; b=[2;1;0]; e=[1;0;1]; wdot=[w(2); q.*b.*w(4)./m;w(4); q.*e./m- q.*b.*w(2)./m; w(6); 0]

14

附录2：

%用Matlab数值模拟的方法模拟带电粒子在恒定磁场中的螺旋运动 function lx global q m B

[t,y]=ode45(@ddlzfun,[0:0.01:20],[0,0.01,0,6,0,0.01],[ ],q,m,B,0);

%用ode23解微分方程组，时间设为20s %指定初始条件，传递相关参数 comet3(y(:,1),y(:,3),y(:,5)); plot3(y(:,1),y(:,3),y(:,5)); grid on %开启坐标网格线

xlabel('x'); ylabel('y'); zlabel('z'); title('模拟带电粒子在恒定磁场中的螺旋运动');

function ydot=ddlzfun(t,y,q,m,B,E) %q,m,B,E为参量 syms q m B q=1.6e-2; B=2; m=0.02;

ydot=[y(2);q*B*y(4)/m;y(4);-q*B*y(2)/m;y(6);0]

15