电磁场的Matlab仿真

1、Matlab与电磁场模拟一 单电荷的场分布：单电荷的外部电位计算公式： 等位线就是连接距离电荷等距离的点，在图上表示就是一圈一圈的圆，而电力线就是由点向外辐射的线。MATLAB程序：theta=0:.01:2*pi;r=0:10;x=sin(theta)*r;y=cos(theta)*r;plot(x,y,b)x=linspace(-5,5,100);for theta=-pi/4 0 pi/4y=x*tan(theta);hold on;plot(x,y); endgrid on单电荷的等位线和电力线分布图：二 多个点电荷的电场情况：模拟一对同号点电荷的静电场 设有两个同号点电荷,其带电量分

2、别为 +Q1和+Q2(Q1、Q20 )距离为 2a则两电荷在点P(x, y)处产生的电势为: 由电场强度可得E = -U,在xOy平面上,电场强度的公式为:为了简单起见,对电势U做如下变换:。Matlab程序：q=1;xm=2.5;ym=2;x=linspace(-xm,xm);y=linspace(-ym,ym);X,Y=meshgrid(x,y);R1=sqrt(X+1).2+Y.2);R2=sqrt(X-1).2+Y.2);U=1./R1+q./R2;u=1:0.5:4;figurecontour(X,Y,U,u)grid onlegend(num2str(u)hold onplot(-

3、xm;xm,0;0)plot(0;0,-ym;ym)plot(-1,0,o,MarkerSize,12)plot(1,0,o,MarkerSize,12)DX,DY = gradient(U);quiver(X,Y,-DX,-DY);surf(X,Y,U);同号电荷的静电场图像为：同理，将程序稍作修改，便可以得到异号电荷的静电场图像：三、线电荷产生的电位：设电荷均匀分布在从z=-L到z=L,通过原点的线段上，其密度为q(单位C/m),求在xy平面上的电位分布。 点电荷产生的电位可表示为 是一个标量。其中r为电荷到测量点的距离。线电荷所产生的电位可用积分或叠加的方法来求。为此把线电荷分为N段，每

4、段长为dL。每段上电荷为q*dL,看作集中在中点的点电荷，它产生的电位为 然后对全部电荷求和即可。 把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

5、*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);endmax(V),min(V)plot(R,V),grid线电荷产生的静电位分布图：四 计算平面上N个电荷之间的库伦引力1 建模： 由库仑定律： 先输入电荷的数目，各电荷的坐标及电荷量，再选一个电荷，求其它电荷对它的作用力，叠加求合力。再选下一个电荷，依次类推。Matlab程序：clear all;N = input(N=:);for ic = 1:N fprintf(-/n#%gn,ic); rc = input(

6、x,y:); x(ic) = rc(1); y(ic) = rc(2); 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:N if(ic = jc) xij = x(ic)-x(jc);yij = y(ic)-y(jc); Rij = sqrt(xij2+yij2); Fx=Fx+C0*q(ic)*q(jc)*xij/Rij3; Fy=Fy+C0*q(ic)*q(jc)*yij/Rij3; end end fprintf(#%gn,ic); fprint

7、f(x-:%en,Fx); fprintf(y-:%en,Fy);end五 有限差分法处理电磁场问题Matlab程序:m=40for k=1:m for j=1:m if k=1 V(j,k)=1; elseif(j=1)|(j=m)|(k=m) V(j,k)=0; else V(j,k)=0.5; end endend cha=0.01;delta=0;n=0;while(1) n=n+1; for k=2:m-1 for j=2:m-1 Vnew(j,k)=1/4*(V(j+1,k)+V(j-1,k)+V(j,k+1)+V(j,k-1); d=abs(Vnew(j,k)-V(j,k)/V(j,k); if ddelta d