电磁场的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 o

3、nplot(-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为电荷到测量点的距离。线电荷所产生的电位可用积分或叠加的方法来求。为此

4、把线电荷分为n段，每段长为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=(

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

6、ic); rc = input(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(

7、#%gn,ic); fprintf(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 delta=d;