matlab在电磁学中的应用

§5-3 带电粒子在电场和磁场中的运动

带电粒子在电场和磁场中的运动及其规律具有重要的应用价值，本节利用MATLAB讨论几个具体实例，读者可以参考这些例子来学习和掌握有关内容。

5.3.1 带电粒子在电场中的运动

电量为q的粒子在电场强度为E的静电场中所受的电场力为

FqE

该力将使质量为m的带电粒子产生一加速度

aF m若带电粒子的初速度为v0，在加速电压U作用下，其动能变化为

Ek1212mvmv0qU 22式中，v为被加速后粒子的末速度。

● 题目(ex5311)

在示波器的竖直偏转系统中加电压于两极板，在两极板之间产生均匀电场E，设电子质量为m，电荷为 －e，它以速度v0射进电场中，v0与E垂直，试讨论电子运动的轨迹。

● 解题分析

电子在两极板间电场中的运动和物体在地球重力场中的平抛运动相似。作用在电子上的电场力为F= －eE，电子的偏转方向与E相反（设为负y方向）。电子在垂直方向的加速度为 aeE。在水平方向和垂直方向电子的运动方程分别为 m11eE2xv0t； yat2t

22m为了讨论电子运动轨迹与初速度及电场的关系，使用了input函数供读者输入

E和v0，以观察不同电场和初速度情况下电子的运动轨迹。

● 程序(ex5311) clear,clf,

238 MATLAB及其在大学物理课程中的应用

E=input('E=','s'); %输入电场强度与时间的函数关系 e=1.6e-19; m=9.1e-31; %给定电子电荷和质量的数值 v0=input('v0='); %输入电子的水平初速度 t=0:0.01:10; %给定时间数组 x=v0.*t;

E1=eval(E); %运算输入的字符串E y=-1./2.*e.*E1.*t.^2./m; plot(x,y,x,0,'r-'),grid on, hold on

运行该程序，在提示后键入E的表达式。例如，取E=100, 或E =100 cos(5t)，便分别得到图 5-3-1和图 5-3-2。从图中可以看出，通过改变极板间的电场可控制带电粒子的运动轨迹。如果在水平方向再加一对极板，则可实现水平和垂直两个方向上的控制，此即示波器的原理。

图 5-3-1恒定电场下的电子运动轨迹 图 5-3-2 变化电场下的电子运动轨迹 ● 题目(ex5312) 一电子二极管由半径r1=0.50 mm的圆柱形阴极K和套在阴极外面的同轴圆筒状阳极A构成。阳极半径R=0.45 cm，阳极电势比阴极高300 V。设电子从阴极发射出来时的速度很小（可忽略不计），试求：

① 电子从K向A走过2.0 mm时的速度； ② 电子到达A时的速度。

● 解题分析 阴极与阳极之间的电位差为

r1r1UAUKEdlR2dr 40rR已知UA－UK＝300 V, 可由此解得单位长度上所带电量λ。

第5章 MATLAB在电磁学中的应用 239

设距阴极2 mm的点为B点，用积分法可求出B点和K点之间的电位差UB-Uk（用UBK表示）。然后，利用

1212e(UAUK) mvBe(UBUK) 和 mvA22求得vB和vA。

● 程序(ex5312)

% (ex53121)计算UA-UK和UB-UK，分别用UAK和UBK表示 syms lambda C0 r R r1; E=-2.*lambda.*C0./r;

UAK=simplify(int(E,'r',R,r1)) UBK=simplify(int(E,'r',2.5,0.5)) 运行结果：

UAK =-2*lambda*C0*(log(r1)-log(R)) UBK =2*lambda*C0*log(5)

% (ex53122)计算单位长度上所带电量lambda,B点的速度VB和到达A点时的速度VA syms VB VA m e lambda R r1;

R=0.0045;r1=0.5;m=9.1e-31;e=1.6e-19;C0=9e9;

lambda=solve('-2*lambda*C0*(log(r1)-log(R))=300','lambda'); lambda1=vpa(subs(lambda),3)

UBK=-2*lambda*C0*log(5); UBK1=vpa(subs(UBK),3);

UAK=300; %已知条件，阳极电势比阴极高300V VB=(2*e*UBK1/m)^(1/2); VB1=vpa(subs(VB),3) VA=(2*e*UAK/m)^(1/2); VA1=vpa(subs(VA),3) 运行结果：

lambda1 =-.692e-8 VB1 =.841e7 VA1 =.103e8

即，电子从K向A走过2.0 mm时的速度vB 以及电子到达A时的速度vA分别为 8.4×106和1.03×107 m / s。

5.3.2带电粒子在电磁场中的运动

如果一个带电粒子在既有电场又有磁场的区域里运动，则其受到的电磁力为

Fq(EvB)

240 MATLAB及其在大学物理课程中的应用

这里，我们进一步讨论既有电场又有磁场情况下带电粒子的运动情况。

● 题目(ex5321)

设质量为m，带电量为q的粒子在磁感应强度为B，电场强度为E的电磁场中运动，建立和求解该问题的微分方程，并分E≠0，B≠0；E=0，B≠0和E≠0，B=0三种情况画出带电粒子在电磁场中的运动轨迹。

● 解题分析

质量为m，带电量为q的粒子在电磁场中的运动微分方程为

d2rm2qEqvB dt选场中某点为原点，以E为O-y方向，B沿O-z方向，建立坐标系Oxyz。令上式的投影方程为

qB，md2xdydt2dtd2yqEdx dt2mdtd2z0dt2令y1=x, y2=dx/dt, y3=y, y4=dy/dt, y5=z, y6=dz/dt, 上述方程可改写为下列一阶微分方程组：

dydy1dyy2;2y4;3y4dtdtdt

dydydy4qEy2;5y6;60dtmdtdt● 程序(ex5321)

%%① (ex53211) 符号法求粒子运动微分方程的特解并绘图 clear

syms w x y z t B E m q;

E=input('E=');B=input('B='); %输入E和B值

[x,y,z]=dsolve('D2x=q*B/m*Dy','D2y=q*E/m-q*B/m*Dx','D2z=0','x(0)=0','y(0)=0','z(0)=0','Dx(0)=0.01','Dy(0)=6','Dz(0)=0.01') ;

第5章 MATLAB在电磁学中的应用 241

%初始条件取x(0)=y(0)=z(0)=0,Dx(0)=0.01,Dy(0)=6,Dz(0)=0.01 q=1.6e-2; m=0.02;

X=subs([x y z]); x=X(1),y=X(2),z=X(3), ezplot3(X(1),X(2),X(3))

运行上述程序，例如，取E=4, B=8可得下列特解并给出图5-3-4（a）。 x =-15/16*cos(32/5*t)-49/640*sin(32/5*t)+1/2*t+15/16 y =15/16*sin(32/5*t)-49/640*cos(32/5*t)+49/640 z =1/100*t (a) E=4, B=8 (b) E=0.01, B=8 (c) E=8, B=0 图 5-3-3 带电粒子在电磁场中的运动 下面我们给出一段用数值方法求解该问题的程序，以便读者比较和练习。 % ② ex5322采用数值方法求解并绘制粒子运动轨迹

q=1.6e-2; m=0.02; B=[2;2;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,y]=ode23('ex5322f',[0:0.1:20],[0,0.01,0,6,0,0.01],[],q,m,B(i),E(i)); axes('unit','normalized','position',[0.0293+(i-1)*0.325 0.062 0.28 0.658]); plot3(y(:,1),y(:,3),y(:,5),'linewidth',2); grid on

title(strd{i},'fontsize',12,'fontweight','demi'); view([-51,18]);

242 MATLAB及其在大学物理课程中的应用

end

函数文件是一个独立的文件，文件名为ex5322f.m function ydot=ex5322f(t,y,flag,q,m,b,e)

ydot=[y(2);q*b*y(4)/m;y(4);q*e/m-q*b*y(2)/m;y(6);0]; 运行该程序,可得到与图5-3-3相同的结果。

§5-4 电流和磁场

如前所述，运动电荷将在其周围产生电磁场，而导体中的电流是由自由电子的定向运动形成的，因此，在载流导体周围也必定有相应的电磁场存在。与电场类似，磁场也遵从叠加原理，所以我们也可将任意形状的载流导体划分成许多电流元；而整个载流导体所产生的磁场就是这些电流元所产生的磁场的叠加。毕奥和萨伐尔在实验的基础上导出了电流元产生的磁感应强度的表达式，称为毕奥－萨伐尔定律。

5.4.1 毕奥－萨伐尔定律

毕奥－萨伐尔定律可表述为：载流回路的任一电流元Idl，在空间任一点P处所产生的磁感应强度dB可表示为

dB0Idlr4r3

其中，r是电流元Idl到场点P的径矢，I为电流。可以看出，dB的方向垂直于Idl与r所在的平面，其指向遵守右手螺旋法则。dB的大小为 dB0Idlsin4r2

利用叠加原理，对上式积分，便可求得任意形状的载流导线所产生的磁感应强度，即 BLdBIdlr

4Lr30 ● 题目(ex5411)

如图5-4-1所示，求垂直于无限长载流直导线的平面内磁感应强度的分布。 ● 解题分析

设场点P的位置为x0iy0jz0k，电流元位置为xiyjzk，电流元矢量

第5章 MATLAB在电磁学中的应用 243

为IdlI(dxidyjdzk)。由此，场点P相对于电流元的位置矢量为

r(x0x)i(y0y)j(z0z)k

i利用行列式计算Idl×r ，可写为 dlrjdyy0ykdzz0z

dxx0x也可利用MATLAB中的det 命令函数来求该行列式，程序如下：

syms dx dy dz x0 x y0 y z0 z; dl=[dx,dy,dz];

r=[x0-x,y0-y, z0-z];

d1cr=cross(dl,r) %求dl×r的积 运行结果为

d1cr = [ dy*(z0-z)-dz*(y0-y), dz*(x0-x)-dx*(z0-z), dx*(y0-y)-dy*(x0-x)] 即 dlr(z0z)dy(y0y)dzi(x0x)dz(z0z)dxj(y0y)dx(x0x)dyk 又，r的大小为

r(x0x)2(y0y)2(z0z)2 设载流导体通过坐标原点垂直于x-y平面放置，电流元Idl沿z轴正向，场点P位于x-y平面上。对本题目而言，dx=dy=0,x=y=0, z0=0, 矢量叉乘积为

dlry0dzix0dzj r的大小为

rx02y02z2 由毕奥－萨伐尔定律

dB有

Idlr4Lr30图5-4-1 题目ex5411示意图 244 MATLAB及其在大学物理课程中的应用

dBx0Iy0dz0Ix0dzdB； y2223/22223/24(x0y0z)4(x0y0z)

Bx0I4y0dz0I； By22(x0y0z2)3/24x0dz 2223/2(x0y0z)BBx2By2

● 程序 (ex5411)

%① 用符号运算求B的表达式(ex54111) syms C0 I z x y r r0;

Bx=C0.*I.*int(-y./(x.^2+y.^2+z.^2).^(3/2),'z',-inf,inf) By=C0.*I.*int(x./(x.^2+y.^2+z.^2).^(3/2),'z',-inf,inf) B=(Bx.^2+By.^2).^0.5 运行结果：

Bx =-2*C0*I*y/(x^2+y^2)^(3/2)/(1/(x^2+y^2))^(1/2) By =2*C0*I*x/(x^2+y^2)^(3/2)/(1/(x^2+y^2))^(1/2)

B =(4*C0^2*I^2*y^2/(x^2+y^2)^2+4*C0^2*I^2*x^2/(x^2+y^2)^2)^(1/2) 即 Bx0I0I2xdz2ydz02IB；； By222222224(xy)4(xy)4xy%② 绘制磁场大小分布图和矢量场图（ex54112）

x=-0.5:0.05:0.5;y=x;

I=input('请输入电流I=');

mu0=4*pi*1e-7; C0=mu0/(4*pi); [X,Y]=meshgrid(x,y);

Bx =-2.*C0.*I.*Y./(X.^2+Y.^2).^(3./2.)./(1./(X.^2+Y.^2)).^(1./2); By =2.*C0.*I.*X./(X.^2+Y.^2).^(3./2)./(1./(X.^2+Y.^2)).^(1./2);

B=(4.*C0.^2.*I.^2.*Y.^2./(X.^2+Y.^2).^2+4.*C0.^2.*I.^2.*X.^2./(X.^2+Y.^2).^2).^(1./2); subplot(1,2,1)

quiver(X,Y,Bx,By,2), axis([-0.5,0.5,-0.5,0.5]), axis('square'), subplot(1,2,2) mesh(X,Y,B)

第5章 MATLAB在电磁学中的应用 245

运行该程序，在命令窗中的提示后键入I0值 (例如，取I0＝100A)，便得到图5-4-2所示图形。

图5-4-2 长载流导线的磁场在x-y平面上的分布

● 题目(ex5412)

设圆线圈的中心为O，半径为R，放置于y-z平面，线圈通过的电流为I0，如图5-4-3所示。用毕奥－萨伐尔定律计算载流圆线圈在z=0处x-y平面上的磁场分布。

● 解题分析 ① 方法1

根据毕奥－萨伐尔定律，

Idlr dB3L4r0 图5-4-3 题目ex5412示意图 线圈上任一点处的电流元在x-y平面上一点

P产生的元磁场为dB。在编制程序时，将电流环分为N段，每一小段视为一电流元，然后求出每一电流元在观察点处的磁场分量，求出总磁场，最后叠加。 ● 程序 (ex54121)

clear all R=input('请输入圆环半径,R='); I0=input('请输入电流,I0=');

246 MATLAB及其在大学物理课程中的应用

mu0=4*pi*1e-7; C0=mu0/(4*pi); %归并常数 N=20; %电流环分段数 x=linspace(-3,3,N); y=x; %确定观测点范围 theta0=linspace(0,2*pi,N+1); %环的圆周角分段 theta1=theta0(1:N);

y1=R*cos(theta1); z1=R*sin(theta1); %环各段矢量的起始坐标y1,z1 theta2=theta0(2:N+1);

y2=R*cos(theta2); z2=R*sin(theta2); %环各段矢量的终点坐标y2,z2

xc=0; yc=(y2+y1)./2; zc=(z2+z1)./2; %计算环各段矢量中点的三个坐标分量xc,yc,zc dlx=0;dly=y2-y1;dlz=z2-z1; %计算环各段矢量dl的三个长度分量，其中x1=x2=0。 NGx=N; NGy=NGx; %网格线数

for i=1:NGy %循环计算各网点上的B（x,y）值 for j=1:NGx

rx=x(j)-xc; ry=y(i)-yc; rz=0-zc; %计算径矢r的3个长度分量，r在z=0的面上。 r3=sqrt(rx.^2+ry.^2+rz.^2).^3; %计算r3

dlXr_x=dly.*rz-dlz.*ry; %计算叉乘dl×r的x和y分量，z分量为0 dlXr_y=dlz.*rx-dlx.*rz;

Bx(i,j)=sum(C0*I0.*dlXr_x./r3); %把环各段产生的磁场分量累加

By(i,j)=sum(C0*I0.*dlXr_y./r3);

B=(Bx.^2+By.^2).^0.5; %计算B的大小 end end

subplot(1,2,1), quiver(x,y,Bx,By), %画矢量场图 hold on

plot(0,1,'ro',0,-1,'bo'),

xlabel('x'),ylabel('y'), %修饰图形，标注坐标轴 axis([-3,3,-3,3]), subplot(1,2,2)

mesh(x,y,B);axis([-3,3,-3,3,0,1e-4]) %画磁场大小分布图 xlabel('x'),ylabel('y'),zlabel('B')

运行该程序，例如，在命令窗中的R和I0的提示后分别键入1和100，所得结果如图5-4-4 所示。

第5章 MATLAB在电磁学中的应用 247

图5-4-4 载流圆线圈的磁场分布 ② 方法2

下面给出用数值积分计算本题的程序，读者可以比较二者的不同，从中体会编程的方法和技巧。

利用题目ex5411中叉乘积的计算表达式，先求出dl×r。

设载流线圈通过坐标原点垂直于x-y平面放置，电流元Idl沿z轴正向，场点P位于x-y平面上。对本题目而言，x= 0，y= Rsinθ，z=Rcosθ，z0=0，dx=0， dy=R cosθdθ，dz= -Rsinθdθ。可以编制下列程序，很方便地求出dl×r。 %求dl×r

syms dx dy dz x0 x y0 y z0 z R theta;

z0=0;dx=0;x=0; y=R.*sin(theta); z=R.*cos(theta); dy=diff(y,theta); dz=diff(z,theta) dl=[dx dy dz]; r=[x0-x,y0-y, z0-z];

dlcr=cross(dl,r) %求矢量的叉乘积 运行后返回： dlcr=

[-R^2*cos(theta)^2+R*sin(theta)*(y0-R*sin(theta)), -R*sin(theta)*x0, -R*cos(theta)*x0] 即

dlrR(Ry0sin)diRx0sindjx0Rcosdk

248 MATLAB及其在大学物理课程中的应用

r的大小为 rx02(y0Rsin)2(Rcos)2 由毕奥－萨伐尔定律

dB有

Idlr 4Lr30dBx0IR(Ry0sin)d 24[x0(y0Rsin)2(Rcos)2]3/2x0Rsind 2223/24[x0(y0Rsin)(Rcos)]dBy0IdBz0Ix0Rcosd 24[x0(y0Rsin)2(Rcos)2]3/2%用数值积分计算本题目的程序(ex54122)

clear figure(1)

R=0.35;y=-1:0.04:1;x=-1:0.04:1; theta=0:pi/20:2*pi; [X,Y,Theta]=meshgrid(x,y,theta);

r=sqrt(X.^2+(Y-R*sin(Theta)).^2+(R.*cos(Theta)).^2); r3=r.^3;

dBx=-R.*(R-Y.*sin(Theta))./r3; dBy=-X.*R.*sin(Theta)./r3;

Bx=pi/40*trapz(dBx,3); By=pi/40*trapz(dBy,3); [BSX,BSY]=meshgrid(0,[0:0.05:0.2]);

h1=streamline(X(:,:,1),Y(:,:,1),Bx,By,BSX,BSY,[0.1,1000]); %绘制流线图 h2=copyobj(h1,gca);

rotate(h2,[1,0,0],180,[0,0,0]); h3=copyobj(allchild(gca),gca); rotate(h3,[0,1,0],180,[0,0,0]); for kk=1:4

[BSX,BSY]=meshgrid(0,0.2+kk*0.02);

streamline(X(:,:,1),Y(:,:,1),Bx,By,BSX,BSY,[0.02/(kk+1),4500]);

第5章 MATLAB在电磁学中的应用 249

streamline(X(:,:,1),-Y(:,:,1),Bx,-By,BSX,-BSY,[0.02/(kk+1),4500]); end

运行结果如图5-4-5所示。

图5-4-5 载流圆线圈的磁力线示意 图5-4-6 题目ex5413示意图

● 题目(ex5413)

一对相同的圆形线圈，彼此平行而共轴。设两线圈内的电流都是I，且回绕方向一致，线圈的半径为R，二者的间距为a（当a=R时，称为亥姆霍兹线圈），求轴线附近的磁场分布。

● 解题分析

本题与题目2类似，只是把观测区域取在两线圈之间的小范围内。如图5-4-6 所示，线圈B生成的左边的磁场等于线圈A的左边磁场。因 此，A、B两线圈在中间部分的合成磁场等于A线圈的右磁场与其左磁场平移R后的和。

● 程序(ex5413) clear all

mu0=4*pi*1e-7;C0=mu0/(4*pi); %归并常数

I0=5.0;R=1; NGx=21;NGy=21; %设定网格线数 x=linspace(-1,1,NGx); %确定观测点范围

y=linspace(-1,1,NGy); N=20; %电流环分段数 theta0=linspace(0,2*pi,N+1); %环的圆周角分段 theta1=theta0(1:N);

y1=R*cos(theta1); z1=R*sin(theta1); %环各段矢量的起始坐标y1,z1

250 MATLAB及其在大学物理课程中的应用

theta2=theta0(2:N+1);

y2=R*cos(theta2); z2=R*sin(theta2); %环各段矢量的终点坐标y2,z2 dlx=0;dly=y2-y1;dlz=z2-z1;

%计算环各段矢量dl的三个长度分量，其中x1=x2=0。 xc=0; yc=(y2+y1)/2; zc=(z2+z1)/2; for i=1:NGy for j=1:NGx

rx=x(j)-xc; ry=y(i)-yc; rz=0-zc; %计算径矢r的3个长度分量,r在z=0的面上。 r3=sqrt(rx.^2+ry.^2+rz.^2).^3;

dlXr_x=dly.*rz-dlz.*ry; %计算叉乘dl×r的x和y分量，z分量为0 dlXr_y=dly.*rx-dlx.*rz;

Bx(i,j)=sum(C0*I0*dlXr_x./r3); %把环各段产生的磁场分量累加 By(i,j)=sum(C0*I0*dlXr_y./r3);

end end

Bax=Bx(:,11:21)+Bx(:,1:11); %把x<0区域内的磁场平移，叠加到x>0区域 Bay=By(:,11:21)+By(:,1:11);

subplot(1,2,1)mesh(x(11:21),y,Bax);xlabel('x');ylabel('y');zlabel('B'); subplot(1,2,2),quiver(x,y,Bx,By,1.5), axis('square'),axis([-1,1,-1,1]), xlabel('x');ylabel('y');

运行结果如图5-4-7 所示。可以看出，在轴线附近磁场大小均匀且沿x方向。

图5-4-7 亥姆霍兹线圈轴线附近Bx在x-y平面上的分布及矢量场 第5章 MATLAB在电磁学中的应用 251

5.4.2 安培定律

如前所述，磁场的基本性质是对处于磁场中的运动电荷施以作用力。因此，如将载流导线置于磁场中，导线中定向运动的自由电子将会受到洛仑兹力的作用，这在宏观上就表现导线为受到了磁场的作用力。

1820年，安培发现了电流之间的相互作用力，得出了电流元之间的的相互作用定律，即安培定律（Ampere law）。该定律的表述为：在磁场中的任一点P处的电流元Idl所受到的磁场作用力dF可表达为

dFIdlB

其中，B是场点P处的磁感应强度，dF也称为安培力，它是洛仑兹力的结果。

对于闭合回路，可将其划分为许多电流元，整个闭合回路所受到的安培力就是各电流元所受到的安培力的矢量和，即 FdFIdlB

● 题目(ex5421)

设两无限长载流导线之间的垂直距离为a，导线中的电流分别为I1和I2。将导线1放置在坐标原点处，电流沿z轴正向。导线2也垂直于x-y平面放置，讨论二者电流的大小和方向与相互作用力的关系。

● 解题分析 由5.4.2 (1)题目1知，导线1在与其垂直距离为a的磁感应强度为

x2y2的导线2处产生

B02I4x2y24a02I

由安培力公式，导线2的一段dl2所受力的大小为

F12I2dl2B10I1I22adl2

同样，导线2 产生的磁场作用在导线1一段dl1上力的大小为

F21I1dl1B20I1I22adl1

由此，在单位长度导线上的作用力的大小为

252 MATLAB及其在大学物理课程中的应用

fF12F210I1I2 dl2dl12a● 程序 (ex5421)

mu0=4*pi*1e-7;

I1=input('I1=');I2=input('I2='); a=0:0.001:1.2;

f=mu0.*I1.*I2./(2.*pi.*a); plot(a,f,a,0,'r-'), hold on, grid on,

axis([0,1.2,-1.5e-6,1.5e-6]) 分两次运行该程序，分别在提示后键入1和－1，得到如图5-4-8 所示的结果。 从图中可以看出，当两导线中电流流向相反时，导线之间表现为排斥力，流向相同时，则表现为吸引力。若两电流大小相等，取a=1米，I=1安培，依次点击图形窗口中的Tools\\Data Cursor，将十字线放在a=1处，对应

的曲线的数据是：x=1, 图5-4-8 两无限长载流导线f与a的关系 y=2e-007。即，此时二导线 －7

间的相互作用力为2×10牛顿／米。反过来，我们用此法可定义电流的单位。事实上，国际计量委员会的正式文件中对电流强度单位“安培”的定义就是：“载有

－

等量电流、相距1米的两个无限长直导线，每米长度上的作用力为2×107 N时，每根导线中的电流强度。”