CENTRAL SOUTH UNIVERSITY题目利用Matlab模拟点电荷电场的分布姓名xxxx学号xxxxxxxxxx班级电气xxxx班任课老师xxxx实验日期2010-10电磁场理论 实验一——利用Matlab 模拟点电荷电场的分布一．实验目的：1．熟悉单个点电荷及一对点电荷的电场分布情况； 2．学会使用Matlab 进行数值计算，并绘出相应的图形；二．实验原理：根据库伦定律：在真空中，两个静止点电荷之间的作用力与这两个电荷的电量乘积成正比，与它们之间距离的平方成反比，作用力的方向在两个电荷的连线上，两电荷同号为斥力，异号为吸力，它们之间的力F 满足：R R Q Q k F ˆ212= (式1)由电场强度E 的定义可知：R R kQ E ˆ2= (式2)对于点电荷，根据场论基础中的定义，有势场E 的势函数为R kQU = (式3)而 U E -∇= (式4)在Matlab 中，由以上公式算出各点的电势U ，电场强度E 后，可以用Matlab 自带的库函数绘出相应电荷的电场分布情况。
三．实验内容：1. 单个点电荷点电荷的平面电力线和等势线真空中点电荷的场强大小是E=kq /r^2 ,其中k 为静电力恒量, q 为电量, r 为点电荷到场点P(x,y)的距离。
电场呈球对称分布, 取电量q> 0, 电力线是以电荷为起点的射线簇。
以无穷远处为零势点, 点电荷的电势为U=kq /r,当U 取常数时, 此式就是等势面方程.等势面是以电荷为中心以r 为半径的球面。
平面电力线的画法在平面上, 电力线是等角分布的射线簇, 用MATLAB 画射线簇很简单。
取射线的半径为( 都取国际制单位) r0=, 不同的角度用向量表示( 单位为弧度)th=linspace(0,2*pi,13)。
射线簇的终点的直角坐标为:[x,y]=pol2cart(th,r0)。
插入x 的起始坐标x=[x; *x].同样插入y 的起始坐标, y=[y; *y], x 和y 都是二维数组, 每一列是一条射线的起始和终止坐标。
用二维画线命令plot(x,y)就画出所有电力线。
平面等势线的画法在过电荷的截面上, 等势线就是以电荷为中心的圆簇, 用MATLAB 画等势线更加简单。
静电力常量为k=9e9, 电量可取为q=1e- 9; 最大的等势线的半径应该比射线的半径小一点 r0=。
其电势为u0=k8q /r0。
如果从外到里取7 条等势线, 最里面的等势线的电势是最外面的3 倍, 那么各条线的电势用向量表示为: u=linspace(1,3,7)*u0。
从- r0 到r0 取偶数个点, 例如100 个点, 使最中心点的坐标绕过0, 各点的坐标可用向量表示: x=linspace(- r0,r0,100), 在直角坐标系中可形成网格坐标: [X,Y]=meshgrid(x)。
各点到原点的距离为: r=sqrt(X.^2+Y.^2), 在乘方时, 乘方号前面要加点, 表示对变量中的元素进行乘方计算。
各点的电势为U=k8q. /r, 在进行除法运算时, 除号前面也要加点, 同样表示对变量中的元素进行除法运算。
用等高线命令即可画出等势线contour(X,Y,U,u), 在画等势线后一般会把电力线擦除, 在画等势线之前插入如下命令hold on 就行了。
平面电力线和等势线如图1, 其中插入了标题等等。
越靠近点电荷的中心, 电势越高, 电场强度越大, 电力线和等势线也越密。
xy单个点电荷的电场线与等势线图1源程序：%点电荷的平面电力线和等势线 %平面电力线的画法 q=1e-9; r0=;th=linspace(0,2*pi,13); [x,y]=pol2cart(th,r0); x=[x;*x]; y=[y;*y]; plot(x,y); grid on hold onplot(0,0,'o','MarkerSize',12) xlabel('x','fontsize',16) ylabel('y','fontsize',16)title('单个点电荷的电场线与等势线','fontsize',20)%平面等势线的画法 k=9e9; r0=; u0=k*q/r0;u=linspace(1,3,7)*u0;x=linspace(-r0,r0,100);[X,Y]=meshgrid(x);r=sqrt(X.^2+Y.^2);U=k*q./r;hold on;contour(X,Y,U,u)clear;点电荷的立体电力线和等势面立体电力线的画法先形成三维单位球面坐标, 绕z 轴一周有8 条电力线[X,Y,Z]=sphere(8),每维都是9×9 的网格矩阵, 将X 化为行向量, 就形成各条电力线的终点x 坐标x=r=X(:)′, 其他两个坐标也可同样形成终点坐标y=r+Y(:)' , z=r+Z(:)' 。
对x坐标插入原点x=[x(zeros(size(x))], 其他两个坐标如下形成y=[y(zeros(size(y))], z=[z(zeros(size(z))], 用三维画线命令plot3(x,y,z), 就画出所有电力线。
立体等势面的画法画5 条等势面时, 各面的电势为u=linspace(1,3,5)+u0, 各等势面的半径为r=k6q. /u, 其中第一个球面的半径为rr=r(1)。
三维单位球面的坐标可由[X,Y,Z]=sphere 命令形成, 每维都是21×21 的网格矩阵, 由于外球会包围内球, 因此把球面的四分之一设为非数, 表示割去该部分Z(X<0&Y<0)=nan. 用曲面命令可画出第一个曲面surf(rr6X,rr6Y,rr6Z), 只要取不同的半径就能画出不同的等势面.为了使等势面好看, 可设置一个颜色浓淡连续变化的命令shading interp。
点电荷的立体电力线和等势面如图2, 旋转图片可从不同的角度观察。
x正电荷电场线等势面的三维图形yz图2源程序：%立体电力线的画法 q=1e-9;[X,Y,Z]=sphere(8); r0=; r1=; k=9e9; u0=k*q/r0; x=r1*X(:)'; y=r1*Y(:)'; z=r1*Z(:)';x=[x;zeros(size(x))]; y=[y;zeros(size(y))]; z=[z;zeros(size(z))]; plot3(x,y,z) hold on;%立体等势线之画法u=linspace(1,3,5)*u0; %画5 条等势面时, 各面的电势为u=linspace(1,3,5)+u0,r=k*q./u; %各等势面的半径为r=k6q. /u[X,Y,Z]=sphere;Z(X<0&Y<0)=nan;surf(r(1)*X,r(1)*Y,r(1)*Z); %第一到第五个球面surf(r(2)*X,r(2)*Y,r(2)*Z);surf(r(3)*X,r(3)*Y,r(3)*Z);surf(r(4)*X,r(4)*Y,r(4)*Z);surf(r(5)*X,r(5)*Y,r(5)*Z);shading interp %个颜色浓淡连续变化的命令shading interp。
xlabel('x','fontsize',16);ylabel('y','fontsize',16);zlabel('z','fontsize',16);title('正电荷电场线等势面的三维图形','fontsize',20);clear;2.一对点电荷平面等势线的画法仍然用MATLAB 的等高线命令画等势线。
对于正负两个点电荷, 电量不妨分别取q1=2e- 9,q2=- 1e- 9, 正电荷在x 轴正方, 负电荷在x 轴负方, 它们到原点的距离定为a=; 假设平面范围为xx0=,yy0=, 两个坐标向量分别x=linspace(- xx0,xx0,20)和y=linspace(- yy0,yy0,50)。
设置平面网格坐标为[X,Y]=meshgrid(x), 各点到两电荷的距离分别为r1=sqrt((X- a).^2+Y.^2)和r2=sqrt((X+a).^2+Y.^2)。
各点的电势为U=k6q1. /r1+k6q2. /r2, 取最高电势为u0=50, 最低电势取其负值。
在两者之间取11 个电势向量u=linspace(u0,- u0,11), 等高线命令contour(X,Y,U,u,'k- ' )用黑实线, 画出等势线如图2所示, 其中, 左边从里到外的第6 条包围负电荷的等势线为零势线。
平面电力线的画法利用MATLAB 的箭头命令, 可用各点的电场强度方向代替电力线。
根据梯度可求各点的场强的两个分量[Ex,Ey]=gradient(- U),合场强为E=sqrt(Ex.^2+Ey.^2)。
为了使箭头等长, 将场强Ex=Ex. /E,Ey=Ey. /E 归一化, 用箭头命令quiver(X,Y,Ex,Ey)可标出各网点的电场强度的方向,异号点电荷对的场点方向如图3 所示。
为了画出连续的电力线, 先确定电力线的起点。
电荷的半径可取为r=, 如图4 所示, 假设第一条电力线的起始角为30 度, 其弧度为q=30+pi /180, 起始点到第一个点电荷的坐标为x1=r0+cos(q),y=r0+sin(q), 到第二个点电荷的坐标只有横坐标x2=2+a+x1 不同。
用前面的方法可求出该点到两个电荷之间的距离r1 和r2, 从而计算场强的两个分量以及总场强Ex=q1+x1 /r1^3 +q2+x2 /r2^3, Ey=q1+y/r1^3+q2+y/r2^3, E=sqrt(Ex6Ex+Ey6Ey)。
下面只要用到场强分量与总场强的比值, 在计算场强分量时没有乘以静电力常量k。
由于电力线的方向与场强的切线方向相同, 取线段为s=,由此可求出终点的坐标为x1=x1+s#Ex/E,y=y+s+Ey/E, 从而计算x2。
以终点为新的起点就能计算其他终点。
当终点出界时或者到达另一点电荷时, 这个终点可作为最后终点. 这种计算电力线的方法称为切线法。
源程序：%一对电荷平面等势线和电场线图clear all;clf;%平面等势线的画法q1=2e-9;q2=-1e-9;a=;%到原点的距离xx0=;yy0=;k=9e9;x=linspace(-xx0,xx0,20);y=linspace(-yy0,yy0,50);[X,Y]=meshgrid(x);r11=sqrt((xx0/^2+(yy0/^2);r22=sqrt((xx0/+a)^2+(yy0/^2);r1=sqrt((X-a).^2+Y.^2); %各点到点电荷的距离r2=sqrt((X+a).^2+Y.^2);U=k*q1./r1+k*q2./r2; %各点的电势u0=k*q1/r11+k*q2/r22;u=linspace(u0,-u0,11); %取21个等势向量contour(X,Y,U,u,'k-');hold ongrid onplot(a,0,'o','MarkerSize',12);plot(-a,0,'o','MarkerSize',12);xlabel('x','fontsize',16);ylabel('y','fontsize',16);%平面电力线的画法[Ex,Ey]=gradient(-U);E=sqrt(Ex.^2+Ey.^2);Ex=Ex./E;Ey=Ey./E;hold on;quiver(X,Y,Ex,Ey);title('一对不相等的电荷的等势线图和电场线图','fontsize',20) clear;图3源程序：%一对电荷平面等势线和电场线图clear all;clf;%平面等势线的画法q1=1;q2=1;a=;xx0=;yy0=;k=9e9;x=linspace(-xx0,xx0,20);y=linspace(-yy0,yy0,50);[X,Y]=meshgrid(x);r11=sqrt((xx0/^2+(yy0/^2);r22=sqrt((xx0/+a)^2+(yy0/^2);r1=sqrt((X-a).^2+Y.^2);r2=sqrt((X+a).^2+Y.^2);U=k*q1./r1+k*q2./r2;u0=k*q1/r11+k*q2/r22;u=linspace(u0,-u0,11);contour(X,Y,U,u,'k-');hold on%平面电力线的画法[Ex,Ey]=gradient(-U);E=sqrt(Ex.^2+Ey.^2);Ex=Ex./E;Ey=Ey./E;dth1=20;th1=(dth1:dth1:180-dth1)*pi/180; r0=a/5;x1=r0*cos(th1)+a;y1=r0*sin(th1);streamline(X,Y,Ex,Ey,x1,y1); streamline(-X,-Y,-Ex,-Ey,x1,-y1); q=abs(q1/q2);dth2=dth1/q;th2=(180-dth2:-dth2:dth2)*pi/180;x2=r0*cos(th2)-a; y2=r0*sin(th2);streamline(X,Y,Ex,Ey,x2,y2); streamline(X,-Y,Ex,-Ey,x2,-y2); grid onplot(a,0,'o','MarkerSize',12); plot(-a,0,'o','MarkerSize',12); xlabel('x','fontsize',16); ylabel('y','fontsize',16); title('一对点电荷的电场分布图'); clear;xy一对点电荷的电场分布图-0.05-0.04-0.03-0.02-0.010.010.020.030.040.05-0.05-0.04-0.03-0.02-0.0100.010.020.030.040.05图4点电荷电场分布的3-D图图5四．实验心得本次电磁场实验是利用Matlab模拟点电荷电场的分布，刚收到实验指导书时，并不知道该怎么做，由于我们并没有正式学过Matlab，只是在部分课程如信号，自控等课上对该软件有所接触。