Matlab 在电磁场中的应用专业: 电气信息与自动化班级：2012级自动化3班学号：***********学院：物电学院指导老师：**完成日期：2013年12月15日Matlab 在电磁场中的应用摘要Matlab是美国Mathworks公司于80年代推出的大型数学软件，通过多年的升级换代，现在已发展成为集数值计算、符号计算、可视化功能以及诸多的工具箱为一体的大型科学计算软件，它已广泛应用于科研院所、工程技术等各个部门，并成为大学生、研究生必备的工具软件。

电磁学是物理学的一个分支，是研究电场和电磁的相互作用现象。

电磁学从原来互相独立的两门科学(电学、磁学)发展成为物理学中一个完整的分支学科，主要是基于电流的磁效应和变化的磁场的电效应的发现。

这两个实验现象，加上麦克斯韦关于变化电场产生磁场的假设，奠定了电磁学的整个理论体系，发展了对现代文明起重大影响的电工和电子技术。

针对电磁场学习理论性强、概念抽象等特点，利用Matlab强大的数值计算和图形技术，通过具体实例进行仿真，绘制相应的图形，使其形象化，便于对其的理解和掌握。

将Matlab引入电磁学中，利用其可视化功能对电磁学实验现象进行计算机模拟，可以提高学习效率于学习积极性，使学习效果明显。

本文通过Matlab软件工具，对点电荷电场、线电荷产生的电位、平面上N 个电荷之间的库仑引力、仿真电荷在变化磁场中的运动等问题分别给出了直观形象的的仿真图，形实现了可视化学习，丰富了学习内容，提高了对电磁场理论知识的兴趣。

关键词：Matlab 电磁学仿真计算机模拟一、点电荷电场问题描述：真空中，两个带正电的点电荷，在电量相同和电量不同情况下的电场分布。

根据电学知识，若电荷在空间激发的电势分布为V，则电场强度等于电势梯度的负值，即：根据题意，真空中若以无穷远为电势零点，则在两个点电荷的电场中，空间的电势分布为：程序实现：clear allep0=8.85*1e-12;c0=1/(4*pi*ep0);e=1.60e-10;h=0.018;x=-0.5:h:0.5;y=-0.5:h:0.5;[X,Y]=meshgrid(x,y);q=[e;1.9*e];for i=1:2V=c0*e./sqrt((X+0.2).^2+Y .^2)+c0.*q(i)./sqrt((X-0.2).^2+Y .^2);[Ex,Ey]=gradient(-V ,h);figure(i)contour(X(:,:,1),Y(:,:,1),V ,...[20,-20,19,-19,18,-18,17,-17,...16,-16,15,-15,14,-14,13,-13,...12,-12,11,-11,10,-10]);axis([-0.38,0.38,-0.28,0.28])hold onphi=0:pi/17:2*pi;sx1=0.2+0.01*cos(phi);sy1=0.01*sin(phi);streamline(X(:,:,1),Y(:,:,1),Ex,Ey,sx1,sy1);hold onsx2=-0.2+0.01*cos(phi);sy2=0.01*sin(phi);streamline(X(:,:,1),Y(:,:,1),Ex,Ey,sx2,sy2);title(str{i})text(-0.212,0,'+','fontsize',20);text(0.187,0,'+','fontsize',20);endE V =-∇1212010244q q V V V r r πεπε=+=+图1-1 两个同号等量电荷的电场分布 图1-2 两个同号不等量电荷的电场分布二、线电荷产生的电位设电荷均匀分布在从z=-L 到z=L,通过原点的线段上，其密度为q(单位C/m),求在xy 平面上的电位分布。

点电荷产生的电位可表示为 0/4V Q r πε= 是一个标量。

其中r 为电荷到测量点的距离。

线电荷所产生的电位可用积分或叠加的方法来求。

为此把线电荷分为N 段，每段长为dL 。

每段上电荷为q*dL,看作集中在中点的点电荷，它产生的电位为04qdL dV r πε=然后对全部电荷求和即可。

把xy 平面分成网格，因为xy 平面上的电位仅取决于离原点的垂直距离R ，所以可以省略一维，只取R 为自变量。

把R 从0到10米分成Nr+1点，对每一点计算其电位。

matlab 程序clear all;L=input(‘线电荷长度L ＝：’);N=input(‘分段数N ＝：’);Nr=input(‘分段数Nr ＝：’);q=input(‘电荷密度q=:’);E0=8.85e-12;C0=1/4/pi/E0;L0=linspace(-L,L,N+1);L1=L0(1:N);L2=L0(2:N+1);Lm=(L1+L2)/2;dL=2*L/N;R=linspace(0,10,Nr+1);for k=1:Nr+1Rk=sqrt(Lm.^2+R(k)^2);Vk=C0*dL*q./Rk;V(k)=sum(Vk);end[max(V),min(V)]plot(R,V),grad输入:线电荷长度L＝：5分段数N＝：50分段数Nr＝：50电荷密度q=:1可得最大值和最小值为:ans =1.0e+010 *[9.3199 0.8654]图（2-1）线电荷产生的静电位分布图三、平面上N个电荷之间的库仑引力建模：由库仑定律：3120/4F q q r πε=其分量的公式可以写成：312210312210()/4()/4x y F q q x x r F q q y y r r πεπε=-=-=编写程序时，先输入电荷的数目，各电荷的坐标及电荷量，再选一个电荷，求其它电荷对它的作用力，叠加求合力。

再选下一个电荷，依次类推。

Matlab 程序：clear all;N = input('输入电荷数目N=:');for ic = 1:N %输入给定条件fprintf('----/n 对电荷#%g\n',ic);rc = input('输入电荷位置[x,y]（米）:');x(ic) = rc(1); %电荷ic 的x 坐标y(ic) = rc(2); %电荷ic 的y 坐标q(ic) = input('输入电荷量（库仑）：');endE0 = 8.85e-12; %真空中的常数C0 = 1/(4*pi*E0); %合并常数for ic = 1:N %循环计每个电荷所受的力Fx = 0.0;Fy = 0.0;for jc = 1:Nif(ic ~= jc)xij = x(ic)-x(jc);yij = y(ic)-y(jc);Rij = sqrt(xij^2+yij^2);Fx = Fx+C0*q(ic)*q(jc)*xij/Rij^3;Fy = Fy+C0*q(ic)*q(jc)*yij/Rij^3;endendfprintf('其它电荷作用在电荷#%g 上的合力为：\n',ic);fprintf('x-分量:%gN\n',Fx);fprintf('y-分量:%gN\n',Fy);end本程序注意学会循环提示并输入参数的方法，以及用双循环解决较复杂的计算过程的编程问题。

输入已知条件:输入电荷数目N=3-------对电荷#1输入电荷位置[x,y](m):[1 2]输入电荷量（库仑）:2-------对电荷#2输入电荷位置[x,y](m):[1 1]输入电荷量（库仑）:1-------对电荷#3输入电荷位置[x,y](m):[3 3]输入电荷量（库仑）:3计算结果:其它电荷作用在# 1 上的合力为:X-分量为:-9.65102e+009NY-分量为1.31581e+010其它电荷作用在# 2 上的合力为:X-分量为:-2.38431e+009NY-分量为-2.03679e+010其它电荷作用在# 3 上的合力为:X-分量为:1.20353e+010NY-分量为7.20982e+009利用matlab软件仿真电荷在变化磁场中的运动程序一%电荷在非均匀磁场中的运动v=10;sita=pi/6; %设定带电粒子的初速度及入射角v=v*cos(sita);u=v*sin(sita); %计算x，y方向的初速度w=0;[t,y] = ode23('yy',[0:0.002:2],[0,v,0,u,0,w]); %求解名为“yy”的微分方程组figure %描绘运动轨迹plot(t,y(:,1)); %绘制一般二维曲线%comet(t,y(:,1)); %绘制二维动态曲线xlabel('t');ylabel('x');figureplot(t,y(:,3));%comet(t,y(:,3));xlabel('t');ylabel('y');figureplot(t,y(:,5));%comet(t,y(:,5));xlabel('t');ylabel('z');figureplot(y(:,3),y(:,5));%comet(y(:,3),y(:,5));xlabel('y');ylabel('z');figureplot3(y(:,1),y(:,3),y(:,5)) %绘制一般三维曲线图%comet3(y(:,1),y(:,3),y(:,5)) %绘制三维动态轨迹xlabel('x');ylabe('y');zlabel('z');%电荷在非均匀磁场中运动的微分方程function f=yy(t,y);global A; %定义全局变量A=100; %设定qB0/mf=[y(2);0;y(4);A*y(6)*y(1);y(6);-A*y(4)*y(1)]; %写入微分方程截图图（4-1）电荷在x轴上运动轨迹图（4-2）电荷在y轴上的运动轨迹图（4-3）电荷在z轴上的运动轨迹图（4-4）电荷在yz平面上的运动轨迹图（4-5）电荷在三维空间中的运动轨接着讨论尖端放电现象function pdemodel[pde_fig,ax]=pdeinit;pdetool('appl_cb',1);set(ax,'DataAspectRatio',[21.103448275862068 15.416666666666664 1]);set(ax,'PlotBoxAspectRatio',[1 1 1]);set(ax,'XLim',[-20.793103448275865 21.41379310344827]);set(ax,'YLim',[-16.5277777777778 14.305555555555529]);set(ax,'XTickMode','auto');set(ax,'YTickMode','auto');% Geometry description:pdecirc(0,0,50,'C1');pdepoly([ -0.36641221374044619,56.061068702290072, 56.610687022900777,... ],[ 1.0992366412213741,1.0992366412213741,-8.5190839694656475,],...'P1');set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','C1-P1')% Boundary conditions:pdetool('changemode',0)pdesetbd(6,'dir',1,'1',...'0')pdesetbd(5,'dir',1,'1',...'0')pdesetbd(4,'dir',...1,'1','0')pdesetbd(3,'dir',...1,'1','0')pdesetbd(2,'dir',1,'1','100')pdesetbd(1,'dir',1,'1','100')% Mesh generation:setappdata(pde_fig,'Hgrad',1.3);setappdata(pde_fig,'refinemethod','regular');pdetool('initmesh')pdetool('refine')pdetool('refine')pdetool('jiggle')pdetool('refine')pdetool('refine')% PDE coefficients:pdeseteq(1,'1.0','0.0','0','1.0','0:10','0.0','0.0','[0 100]') setappdata(pde_fig,'currparam',['1.0';'0.0';'0 ';'1.0'])% Solve parameters:setappdata(pde_fig,'solveparam',...str2mat('0','95232','10','pdeadworst',...'0.5','longest','0','1E-4','','fixed','Inf'))% Plotflags and user data strings:setappdata(pde_fig,'plotflags',[4 1 1 2 1 1 1 1 0 0 0 1 1 1 0 1 0 1]); setappdata(pde_fig,'colstring','u');setappdata(pde_fig,'arrowstring','');setappdata(pde_fig,'deformstring','');setappdata(pde_fig,'heightstring','');% Solve PDE:pdetool('solve')劈尖带电50V ，由图可见在尖端出的场强明现比别出大。