解答：如何用matlab绘制三维空间球体的电场线

有朋友可以解答一下如何用matlab绘制三维空间球体的电场线吗，有偿

1 可视化方法

在MATLAB中，用于场的可视化的函数有streamline、contour、contourf、quiver。除此之外，mesh、surf、imagesc等很多函数也可以用于场的可视化，这里主要有下面四种。

函数名称 作用 调用方式
streamline 根据向量数据绘制流线图 streamline(X,Y,Z,U,V,W,startx,starty,startz)
contour 创建一个包含矩阵 Z 的等值线的等高线图 contour(X,Y,Z)
contourf 创建一个包含矩阵 Z 的等值线的填充等高线图 contourf(X,Y,Z)
quiver 将向量显示为箭头，用于向量场的可视化 quiver(x,y,u,v)
以下为官方文档中的这三种函数的示例，熟悉函数的使用。

streamline
[x,y] = meshgrid(0:0.1:1,0:0.1:1);
u = x;
v = -y;
startx = 0.1:0.1:1;
starty = ones(size(startx));
streamline(x,y,u,v,startx,starty)

Figure1-1 streamline绘图
contour
注：contourf用法和contour类似。
[X,Y,Z] = peaks;
contour(X,Y,Z,20)

Figure1-2 contour绘图
quiver
[x,y] = meshgrid(0:0.2:2,0:0.2:2);
u = cos(x).*y;
v = sin(x).*y;
quiver(x,y,u,v)

以下是代码：
clc
clear

a = 0.4; E0 = 2;
k1 = 0.05; k0 = 0.01;

[X,Y] = meshgrid(0:0.01:1);
[theta,r] = cart2pol(Y,X);
rout = r; rin = r;
rin(rin>a) = NaN; % 圆内部（设置外部为NaN）
rout(rout<a) = NaN; % 圆外部（设置内部为NaN）
Zout = -E0rout.cos(theta)+((k1-k0)/(k1+2k0))E0a^3cos(theta)./rout.^2;
Zin = -(3k0/(k1+2k0))E0rin.*cos(theta);

L1 = isnan(Zout); L2 = isnan(Zin);
Zout(L1==1) = 0; Zin(L2==1) = 0;
Z = Zout+Zin;
z = [-flip(Z);Z];
Z0 = [fliplr(z) z];

[PX,PY] = gradient(-Z0,0.05);
[x,y] = meshgrid(-1:0.01:1+0.01);
vy = 0; vx = 0.1:0.1:1;
[Vx,Vy] = meshgrid(vx,vy);

hold on
[C,h] = contourf(x,y,Z0,100);
set(h,'LineStyle','none')
M = 20;
quiver(x(M:M:end-M,M:M:end-M),y(M:M:end-M,M:M:end-M),PX(M:M:end-M,M:M:end-M),PY(M:M:end-M,M:M:end-M),1,'r');
contour(x,y,Z0,8,'b','ShowText','on');
rectangle('Position',[0-a,0-a,2a,2a],'Curvature',[1,1])
axis equal
colorbar
streamline(x,y,PX,PY,Vx,Vy);
streamline(-x,y,-PX,PY,-Vx,Vy);
streamline(x,-y,PX,-PY,Vx,-Vy);
streamline(-x,-y,-PX,-PY,-Vx,-Vy);

matlab绘制三维空间球体的电场线:

clc,clear,close all;
daotiqiudianchang
function daotiqiudianchang()
E0=1;%电场方向默认z轴
epsilon0=8.854187e-12;%真空中电常数
a=1;%导体球半径
Q=0;%导体带电量
syms r theta phi%采用球坐标
Uq=-(Q/(4*pi*a*epsilon0)+E0*r*cos(theta))+Q/(4*pi*epsilon0*r)+E0*a^3/r^2*cos(theta);%球外电势,设r=a处U=0
syms x y z%直角坐标
U=(E0*a^3*z)/(x^2 + y^2 + z^2)^(3/2) - E0*z + Q/(4*epsilon0*pi*(x^2 + y^2 + z^2)^(1/2)) - Q/(4*a*epsilon0*pi);
%以下是3D的电势分布
[xs,ys,zs]=sphere;
figure('Name','电势');
surf(xs,ys,zs,'FaceAlpha',.3,'EdgeColor','none');
hold on
title('外电场方向:z轴');
[x1,y1,z1]=meshgrid(-3*a:a/7:3*a,-3*a:a/30:3*a,-3*a:a/30:3*a);
U1=zeros(size(x1));
for k=1:numel(x1)
    if x1(k)^2+y1(k)^2+z1(k)^2>a^2
        %U1(k)=double(subs(U,[x,y,z],[x1(k),y1(k),z1(k)]));
        U1(k)=(E0*a^3*z1(k))/(x1(k)^2 + y1(k)^2 + z1(k)^2)^(3/2) - E0*z1(k) + Q/(4*epsilon0*pi*(x1(k)^2 + y1(k)^2 + z1(k)^2)^(1/2)) - Q/(4*a*epsilon0*pi);
    end
end
xslice=-3*a:a/7:3*a;
contourslice(x1,y1,z1,U1,xslice,[],[],9);
axis equal
hold off
%以下是3D的电场
figure('Name','电场');
title('外电场方向:z轴');
hold on
surf(xs,ys,zs,'FaceAlpha',.3,'EdgeColor','none');
%电场线的疏密程度和电场场强有关，电场线密集的地方电场场强大
%算导体表面的电场，电场线的起始点的密集程度应与导体表面的电场成正相关
Ebiaomian=subs(diff(Uq,r),r,a);
frho=abs(Ebiaomian*sin(theta));
thetaqu=0:pi/1000:pi;
thetaquan=double(subs(frho,theta,thetaqu));
randtheta=rand(size(thetaqu));
theta2=thetaqu(randtheta*30<thetaquan);%电场线起始点的theta
phi2=0:pi/numel(theta2):2*pi-pi/numel(theta2);%电场线起始点的phi
diffr=diff(Uq,r);
difftheta=diff(Uq,theta)/r;
diffrtheta=sqrt(diffr^2+difftheta^2);
diffr1=diffr/diffrtheta;
difftheta1=difftheta/diffrtheta;
for k=theta2(1:end)
    theta3=k;
    r3=a;
    if subs(diffr,[r,theta],[r3,theta3])>0%流进
        liu=1;
    else%流出
        liu=-1;
    end
    while abs(r3(end)*cos(theta3(end)))<2*a && abs(r3(end)*sin(theta3(end)))<4*a
        r3=[r3 r3(end)+liu/20*double(subs(diffr1,[r,theta],[r3(end),theta3(end)]))];
        theta3=[theta3 theta3(end)+liu/20*double(subs(difftheta1,[r,theta],[r3(end),theta3(end)]))];
    end
    for phi3=phi2(1:end)
        x3=r3.*sin(theta3)*cos(phi3);
        y3=r3.*sin(theta3)*sin(phi3);
        z3=r3.*cos(theta3);
        plot3(x3,y3,z3);
    end
end
axis equal
view(3);
hold off
end

1 浅析 Matlab 模拟静电场三维图 与恒定电流场模拟静电场 摘要：本文利用恒定电流场模拟静电场，通过寻找等势点的方法描绘出点 电荷对的电力线和电势面；而后利用 Matlab 来实现点电荷对在三维空间里的电 力线以及等势面的描绘，并且给出其详尽的计算程序以及注释，使 Matlab 初学 者能够轻松的看懂程序；同时对 Matlab 模拟静电场和恒定电流模拟静电场两种 方法描绘点电荷对之间的电力线图以及等势线图进行分析对比。 关键词：Matlab ；电力线；等势面；三维 1 引言 对于静电场的描绘有很多方法以及改进。代伟等人对传统的恒定电流法模 拟静电场的实验做出了导电介质、等位点观测以及等位点记录等方面做了改进， 使实验结果更加精确 [1] 。而对于 Matlab 描绘静电场中，王明美利用 streamline 命令描绘出了一对点电荷的二维电力线和等势线 [2] 。王静将两点电荷的电荷量 改为比值，对 Matlab 描绘静电场实验进行了优化 [3] 。周胜利用循环和 ode45 解 微分方程的方法描绘出点电荷的电场 [4] 。张雅男等人对恒定电流模拟静电场和 matlab 模拟静电场二维情况下绘制出的图形进行比较，并且通过分析得出两种 方法所得的结果相似却并不完全一致 [5] 。 本文通过比较 matlab 来模拟描绘电荷对之间的静电场的方法与恒定电流法 描绘静电场的方法，对两种实验的原理、过程以及结果进行比较，进而了解两 种方法之间的区别、联系以及优缺点。 2 利用恒定电流场模拟静电场 2.1 简介恒定电流场模拟静电场实验原理 带电体在周围空间产生的电场可以用电场强度 E 或者电势 U 来描述。由于 静电场中不会有电流，不能够用直流电表直接测量。而静电式仪表要用到金属 制的探头，当探头伸入静电场中时，静电场会发生显著变化。不能够直接在静 电场中绘制等势线。而从静电场和电流场都引入电势 U ，都遵守高斯定理等相 似的地方，所以可以利用恒定电流场来对静电场进行模拟 [6] 。 2.2 恒定电流场模拟静电场实验 当绘制点电荷对电场时，通过两个电极接到导电介质上，再在电极上加上 恒定直流电压，就可以得到了恒定电流场。 导电介质可以选取导电纸、水、导电玻璃等，本文选用的导电介质是导电 纸。 实验结果可以利用等臂记录法、复写纸法、放大尺法等方法来记录。本文2 利用了补偿法电路 [6] 和复写纸法来寻找等势点并减小误差。并且绘制出了等量 异号点电荷对形成的等势线以及电力线，并且取点在 excel 中拟合出图形，如图 1。图1 等量异种点电荷的等势线和电力线 Fig.1 The power line and potential of a pair of diffient class equivalent point charges 图 1 显示：等量异种点电荷等势线越靠近电荷越密集。电力线起于正电荷 终于负电荷。 3 利用 Matlab 模拟静电场 3.1 简介 Matlab 部分编程命令 Plot3 是画三维曲线的命令，可以描绘出空间中立体电力线。 Surf 是将三维网格连成曲面的命令，可以形成三维空间下的电势面。 Contour 是等高线命令，可以画出平面等势线。 Gradient 是求梯度的命令。由于电场强度是电势的负梯度 [7] 公式： ，利用命令[Ex,Ey]=gradient(-U),求出电场在空间各点的 x 分量和 y      n e n V E 分量。 Ode45 是 matlab 中一个常用的解微分方程的命令 [8] 。 3.2 实现 Matlab 模拟静电场编程 3.2.1 点电荷对电力线画法 常用的点电荷对电力线画法有两种：第一种叫做切线法，第二种是解微分 方程 [3] ，本文应用第二种方法。 设电荷量为 q1、q2 的两点电荷在(-1,0,0)处和(1,0,0)处，空间任意一点 p(x,y) 。由于电场里面任意一点电场线的切线方向就是该点的场强方向，可以3 得到： ，引入参变量 t： ，利用库伦定理和场强叠加原理， dy dx Ey Ex  t Ey dy Ex dx   则可以求出两点电荷在 p 点的场强分别为：2 3 2 2 1 1 ] ) 1 [( ] ) 1 [( y x yj i x q k E      2 3 2 2 2 2 ] ) 1 [( ] ) 1 [( y x yj i x q k E      计算其和场强为： j E i E j y x y q y x y q k i y x x q y x x q k E E E y x                   } ] ) 1 [( ] ) 1 [( { } ] ) 1 [( ) 1 ( ] ) 1 [( ) 1 ( { 2 3 2 2 2 2 3 2 2 1 2 3 2 2 2 2 3 2 2 1 2 1 由此我们可以得到电力线的微分方程： 2 3 2 2 2 2 3 2 2 1 ] ) 1 [( ) 1 ( ] ) 1 [( ) 1 ( y x x kq y x x kq E dt dx x          2 3 2 2 2 2 3 2 2 1 ] ) 1 [( ] ) 1 [( y x y kq y x y kq E dt dy y        在计算公式中静电力常量 ，由于我们运用 matlab 模拟绘 2 2 9 10 0 . 9      C m N k 图，可以将 k 值取为 1，所得出的静电场图形不变 [3] 。 将此微分方程编成函数文件： function dxdy=fun1(t,p,flag,q1,q2); dxdy=[q1p(1)./(sqrt((p(2)+1).^2+p(1).^2).^3)+q2p(1)./(sqrt((p(2)- 1).^2+p(1).^2).^3); q1*(p(2)+1)./(sqrt((p(2)+1).^2+p(1).^2).^3)+q2*(p(2)-1)./(sqrt((p(2)- 1).^2+p(1).^2).^3)]; 命名为 fun1.m 。 接下来利用上面编辑好的微分方程函数来绘出等量同种点电荷对的电力线。 首先可以将电荷量设为 e 的倍数，我们在输入电荷量的时候就可以简化为输入 实数来描绘静电场了。 clear,clc,close all % 清除命令 q1=2;q2=2; %确定两点电荷的电荷量 a=1; %设定两点电荷到原点的距离 a0=0.1; %设定点电荷的半径 figure (1); %建立图形窗口 1 box on;