MATLAB中使用streamline函数绘制正负点电荷及它们构成的电偶极子的电场线分布图

电场强度等于电势梯度乘以-1，而这些代码在调用gradient函数求偏导后并没有乘以-1。
先把上一篇博客中不严谨的代码拿出来：

% 在二维平面上绘制一个正点电荷的电场线图。

k = 8.9875e+9; % 比例系数：

e = 1.602e-19; % 指定点电荷电量为一个电子带电量绝对值

ke = k * e; % k by e

d = -2:0.05:2;

[x, y] = meshgrid(d, d);

% 计算电位

V = ke ./ sqrt(x.^2 + y.^2 + 0.01); % 加了一个0.01，防止分母为0。

% 求电势的梯度，即电场强度

[E_x, E_y] = gradient(V); % 警告：这里没有乘以-1

sx = [linspace(-1, 1, 10), linspace(1, 1, 10), linspace(-1, 1, 10), linspace(-1, -1, 10)];

sy = [linspace(1, 1, 10), linspace(-1, 1, 10), linspace(-1, -1, 10), linspace(-1, 1, 10)];

hold on;

streamline(x, y, E_x, E_y, sx, sy);

contour(x, y, V, linspace(min(V(:)), max(V(:)), 60));

hold off;

运行结果：



图片说明：其中， 由中心向四周辐射的射线就是描述正点电荷电场的电场线，等势线用色彩渐变的同心圆环表示（中央圆环颜色为暖色系，说明越靠近正电荷表面，电势越大；外围圆环颜色为冷色系，说明越远离正电荷表面，电势越小）。

代码说明：在以(-1,1) (1,1) (1,-1) (-1,-1)为顶点的正方形的四边上均匀取40个点作为起点，并绘制40条电场线（如果您数了一下不够40条的话，说明它们中有些线条发生了重叠）。

我们知道，正电荷的电场线由电荷发出，然后指向无穷远， 但在上面的代码中，由于我们求梯度的时候没有乘以-1，所以， 我们调用streamline函数绘制的电场线由sx、sy指定的起点向正电荷所在的原点汇聚，这造成一种我们成功绘

制了正点电荷的电场线分布图的假象。

现在，让我对症下药,

[E_x, E_y] = gradient(V);

改为：

[E_x, E_y] = gradient(-V);

运行结果：



这次，电场线从(sx,sy)出发，并指向无穷远。在靠近正电荷的四周便显得空荡荡的。要解决这个问题，我们需要把绘制流线的起点指定在正电荷表面，并在对电势求完梯度后乘以-1，这样绘制出来的电场线就从正电荷出发，指向无穷远。

下面是修改后的代码：

close all; clear; clc;

% 在二维平面上绘制一个正电荷的电场线图。

k = 8.9875e+9; % 比例系数

e = 1.602e-19;% 一个电子带电量绝对值

e_r = 2.8e-15; % 正电子的半径

ke = k .* e; % k by e

% 指定区间： d=<x,y<=d，并生成网格数据

d = -e_r*40:e_r:e_r*40;

[x, y] = meshgrid(d);

% 计算电位

V_max = ke / e_r; % 正电荷的表面电势最大，无穷远处电势为0

V = ke ./ sqrt(x.^2 + y.^2);

V(V==Inf) = V_max;

V(V>V_max) = V_max;

% 求电势的梯度，即电场强度

[E_x, E_y] = gradient(-V);

% 绘制电场线分布

t = 0:2*pi/20:2*pi;

sx = e_r * cos(t);

sy = e_r * sin(t);

hold on;

plot(0, 0, 'ro', 0,0, 'r+'); % 标出正电荷

streamline(x, y, E_x, E_y, sx, sy);

contour(x, y, V, linspace(min(V(:)), max(V(:)), 60));

axis([min(d), max(d), min(d), max(d)]);

hold off;

运行结果：



正点电荷的电场线分布图绘制出来了， 那么如何绘制负点电荷的电场线分布图呢？ 把电量e添加一个负号就行了吗？ 不行。

在调用streamline函数绘制电场线分布图时， 指定起点也相当重要。有了绘制正点电荷电场线分布图修改起点的发。 于是，聪明的你一定想到了：负点电荷的电场线从无穷远指向电荷本身，我们可以把绘制流线图的起点定在无穷远处（只要绘制流线时指定的起点离电荷的距离比电荷的半径大很多就算是无穷远处了）。

close all; clear; clc;

% 在二维平面上绘制一个负电荷的电场线图。

k = 8.9875e+9; % 比例系数

e = -1.602e-19; % 一个电子带电量

e_r = 2.8e-15; % 电子的半径

ke = k .* e; % k by e

% 指定区间： d=<x,y<=d，并生成网格数据

d = -e_r*40:e_r:e_r*40;

[x, y] = meshgrid(d);

V_min = ke / e_r; % 负电荷的表面电势最小，无穷远处电势为0

V = ke ./ sqrt(x.^2 + y.^2);

V(V==-Inf) = V_min;

V(V<V_min) = V_min;

% 求电势的梯度，即电场强度

[E_x, E_y] = gradient(-V);

% 绘制电场线分布

% 对于负电荷， 电场线从无穷远处指向电子，所以指定起点要指定为"无穷远处"

t = 0:2*pi/20:2*pi;

sx = max(d) * cos(t);

sy = max(d) * sin(t);

hold on;

plot(0, 0, 'ro', 0, 0, 'r-', 'MarkerSize', 10); % 标出负电荷

streamline(x, y, E_x, E_y, sx, sy);

contour(x, y, V, linspace(min(V(:)), max(V(:)), 60));

axis([min(d), max(d), min(d), max(d)]);

hold off;

运行结果：



最后，贴出来一对电偶极子的电场线分布图：

close all; clear; clc;

% 在二维平面上绘制一对电偶极子的电场线图。

k = 8.9875e+9; % 比例系数

e_p = 1.602e-19; % 正点电荷带电量

e_n = -e_p; % 负点电荷带电量

e_r = 2.8e-15; % 电荷的半径

% 指定区间： d=<x,y<=d，并生成网格数据

d = -e_r*40:e_r:e_r*40;

[x, y] = meshgrid(d);

dt = (max(d) - min(d)) / 10; % 设定两个电子间的距离

x_n = -dt / 2; y_n = 0; % 设定负电子的坐标值

x_p = dt / 2; y_p = 0; % 设定正电子的坐标值

% 分别计算正负电荷在周围一点的电势

V1_min = k * e_n / e_r; V2_max = k * e_p / e_r;

V1 = k * e_n ./ sqrt((x-x_n).^2 + (y-y_n).^2); % 负电荷

V2 = k * e_p ./ sqrt((x-x_p).^2 + (y-y_p).^2); % 正电荷

V1(V1==-Inf) = V1_min; V1(V1<V1_min) = V1_min;

V2(V2==Inf) = V2_max; V2(V2>V2_max) = V2_max;

% 利用叠加原理计算电势

V = V1 + V2;

[E_x, E_y] = gradient(-V);

hold on; grid on;

% 电偶极子一部分电场线从正点电荷出发， 并汇聚到负点电荷

% 绘制从正电荷发出的电场线， 这些电场线一部分汇聚到负点电荷

% 还有一部分射向无穷远

t = linspace(-pi, pi, 25);

sx = e_r * cos(t) + x_p; sy = e_r * sin(t) + y_p;

streamline(x, y, E_x, E_y, sx, sy);

% 为负电荷补充几条射向无穷远的电场线

sx = [min(d)/3*2, min(d), min(d), min(d), min(d)/3*2];

sy = [min(d), min(d)/3*1, 0, max(d)/3*1, max(d)];

streamline(x, y, E_x, E_y, sx, sy);

contour(x, y, V, linspace(min(V(:)), max(V(:)), 60)); % 绘制等势线

plot(x_n, y_n, 'ro', x_n, y_n, 'r-', 'MarkerSize', 8); % 标出负电荷

plot(x_p, y_p, 'ro', x_p, y_p, 'r+', 'MarkerSize', 8); % 标出正电荷

axis([min(d), max(d), min(d), max(d)]);

title('E-field of an electric dipole');

hold off;

运行结果：