正負の点電荷がつくる電場の向きを表す方向場を Python の SymPy Plotting Backends (SPB) で描く。
正負の点電荷がつくる電場
$\boldsymbol{r}_1$ にある正電荷 $q\ (>0)$ と，$\boldsymbol{r}_2$ にある負電荷 $-q$ がつくる電場ベクトル $\boldsymbol{E}$ は
\begin{eqnarray}
\boldsymbol{E}
&=& \frac{q}{4\pi \varepsilon_0} \frac{\boldsymbol{r} – \boldsymbol{r}_1}{|\boldsymbol{r} – \boldsymbol{r}_1|^3}
– \frac{q}{4\pi \varepsilon_0} \frac{\boldsymbol{r} – \boldsymbol{r}_2}{|\boldsymbol{r} – \boldsymbol{r}_2|^3}
\end{eqnarray}
$z=0$ の $xy$ 平面で2次元ベクトル場を描く。適宜無次元化を行い，
$\displaystyle \frac{q}{4\pi \varepsilon_0}\rightarrow 1$ とする。また，
$$\boldsymbol{r} = (x, y, 0), \quad \boldsymbol{r}_1 = (0, 1, 0),
\quad \boldsymbol{r}_2 = (0, -1, 0)$$
とすると，
$$E_x = \frac{x}{\left(x^2 + (y-1)^2\right)^{3/2}}-\frac{x}{\left(x^2 + (y+1)^2\right)^{3/2}}$$$$E_y = \frac{y-1}{\left(x^2 + (y-1)^2\right)^{3/2}}-
\frac{y+1}{\left(x^2 + (y+1)^2\right)^{3/2}}$$
$\boldsymbol{E}$ の向きを表す方向場（単位ベクトル場）
一般にベクトルは大きさと向きを表す。電場 $\boldsymbol{E}$ の大きさは電荷からの距離の2乗に反比例するので，ベクトルの長さが大きく変化する。見やすくするために，以下のように $\boldsymbol{E}$ の向きを表す方向場，すなわち大きさが $1$ の単位ベクトル場を定義して，ベクトルの向きのみを正しく描画することも必要になる場合があるかもしれない。
$$\hat{\boldsymbol{E}} \equiv \frac{\boldsymbol{E}}{\sqrt{\boldsymbol{E}\cdot\boldsymbol{E}}}$$$$\hat{E}_x = \frac{E_x}{\sqrt{E_x^2 + E_y^2}}$$$$\hat{E}_y = \frac{E_y}{\sqrt{E_x^2 + E_y^2}}$$
SPB の場合は別途定義しなくても，normalize = True とすれば適宜規格化された方向ベクトル場を描くことができる。
ライブラリの import
from sympy import *
# 1文字変数の Symbol の定義が省略できる
from sympy.abc import *
# SymPy Plotting Backends (SPB)
from spb import *
# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']
SPB で電場の方向場を描く
電場ベクトルの定義
電場ベクトルを x, y の2変数関数として定義。
def Ex(x, y):
return x/sqrt(x**2 + y**2)**3
def Ey(x, y):
return y/sqrt(x**2 + y**2)**3
def E2x(x, y):
return Ex(x, y-1) - Ex(x, y+1)
def E2y(x, y):
return Ey(x, y-1) - Ey(x, y+1)
def E(x, y):
return sqrt(E2x(x, y)**2 + E2y(x, y)**2)
ベクトル場を描く
SPB において2次元ベクトル場を描く際の書式は…
plot_vector([x成分, y成分], (x, xmin, xmax), (y, ymin, ymax), ...);
デフォルトでは，ベクトルの大きさの contour plot や colormap も描くので，まずはそれらを False にし，normalize = True として方向ベクトル場のみを描いてみる。このオプションによって，ベクトルの長さは同じ長さに自動調整する。
plot_vector(
[E2x(x, y), E2y(x, y)],
(x, -5.5, 5.5), (y, -5.5, 5.5),
# contour plot なし，colormap なし，長さ自動調整して
scalar = False, use_cm = False, normalize = True,
size = (6.4, 6.4)
);
上記の例では，ベクトルの本数（サンプル数）を特に指定せず自動で描かせている。点電荷が $y$ 軸上にあるので，ベクトルの始点が点電荷上にならないように，サンプル数を調整してみる。
以下の例では，$x, y$ 座標ともに
$$-5.5, -4.5, \cdots, -0.5, 0.5, \cdots, 5.5$$
となり，$y$ 軸上 ($x = 0$) に始点がこない。ベクトル場の始点がこのサンプル位置にあることがよくわかるように，グリッド間隔を 1 にしてみる。
{"pivot":"mid"} オプションを使うと，ベクトルの真ん中を中心にしてベクトルの向きを変えて表示します。
p = plot_vector(
[E2x(x, y), E2y(x, y)],
# -5.5 から 5.5 までを...
(x, -5.5, 5.5), (y, -5.5, 5.5),
# 等分して 12 点で評価する
n1 = 12, n2 = 12,
# 始点ではなく真ん中を固定して向きを変える
quiver_kw={"pivot":"mid"},
scalar = False, use_cm = False, normalize = True,
xlim = (-6.5, 6.5), ylim = (-6.5, 6.5),
size = (6.4, 6.4), show = False
)
ax = p.ax
# -6 から (7未満の) 6 まで x の目盛を 1 刻みに
ax.set_xticks([i for i in range(-6,7)])
# -6 から (7未満の) 6 まで y の目盛を 1 刻みに
ax.set_yticks([i for i in range(-6,7)]);
点電荷の位置にラベル
正電荷，負電荷の位置にラベルをつけてみる。
p = plot_vector(
[E2x(x, y), E2y(x, y)],
(x, -5.5, 5.5), (y, -5.5, 5.5),
n1 = 12, n2 = 12,
quiver_kw={"pivot":"mid"},
scalar = False, use_cm = False, normalize = True,
xlim = (-6.5, 6.5), ylim = (-6.5, 6.5),
size = (6.4, 6.4), show = False
)
ax = p.ax
ax.set_xticks([i for i in range(-6,7)])
ax.set_yticks([i for i in range(-6,7)])
# 正電荷の位置に青丸と +
ax.scatter(0, 1, marker='o',
# size, color
s = 200, c = 'white',
edgecolors = 'blue', linewidths = 2)
ax.text(0, 1, "＋", color = "blue",
ha = "center", va = "center", size=12)
# 負電荷の位置に赤丸と -
ax.scatter(0, -1, marker='o',
# size, color
s = 200, c = 'white',
edgecolors = 'red', linewidths = 2)
ax.text(0, -1, "ー", color = "red",
ha = "center", va = "center", size=12);
カラーマップ (use_cm = True)
デフォルトでは use_cm = True なので，強弱（ベクトルの大きさ $E(x,y)$）に応じてベクトルに色がつく。
p = plot_vector(
[E2x(x, y), E2y(x, y)],
(x, -5.5, 5.5), (y, -5.5, 5.5),
n1 = 12, n2 = 12,
quiver_kw={"pivot":"mid"},
scalar = False, normalize=True,
xlim = (-6.5, 6.5), ylim = (-6.5, 6.5),
size = (6.4, 6.4), show = False
)
ax = p.ax
ax.set_xticks([i for i in range(-6,7)])
ax.set_yticks([i for i in range(-6,7)]);
カラーマップ の色味，ベクトルの大きさ・太さの設定
電荷の位置から離れると，電場の大きさは急激に小さくなるので，このままではほとんど同じ色のベクトルばかりになる。
これだと少し寂しいのでカラーマップの範囲を設定する clim を調整して，離れた位置のベクトルもそこそこの青色になるように調整する。
上記の例では clim = (0, 2) くらいであろうから，ゼロに近い大きさでも色がつくように手探りで数値をいれてみます。
また，ベクトルそのものの大きさ・太さも少し変更してみる。書式はまず units='xy' などと設定後，ベクトルの長さを伸ばしたいなら scale により小さい値を，太くしたいなら width により大きい値をいれます。これも手探りでいれてみました。
右側の colorbar を非表示にするには colorbar = False です。
quiver_kw に設定可能なオプションは matplotlib.pyplot.quiver() と同じ。以下を参照。
p = plot_vector(
[E2x(x, y), E2y(x, y)],
(x, -5.5, 5.5), (y, -5.5, 5.5),
n1 = 12, n2 = 12,
scalar = False, use_cm = True, normalize=True,
xlim = (-6.5, 6.5), ylim = (-6.5, 6.5),
size = (6.4, 6.4), show = False,
# colorbar 非表示
colorbar=False,
quiver_kw = {
"cmap":"Blues",
"clim":(-0.1, 0.8),
"pivot":"mid",
"units":"xy", "scale":1, "width":0.12
}
)
ax = p.ax
ax.set_xticks([i for i in range(-6,7)])
ax.set_yticks([i for i in range(-6,7)]);
一応の完成作
以上の設定をまとめてみたのが以下の例。少し短めのベクトルをたくさん描いてみた。
p = plot_vector(
[E2x(x, y), E2y(x, y)],
(x, -5.25, 5.25), (y, -5.25, 5.25),
n1 = 22, n2 = 22,
scalar = False, use_cm = True, normalize=True,
colorbar = False,
xlim = (-5, 5), ylim = (-5, 5),
size = (6.4, 6.4), show = False,
legend = False,
quiver_kw = {
"cmap":"Blues",
"clim":(-0.1, 0.8),
"pivot":"mid",
"units":"xy", "scale":2, "width":0.06
}
)
ax = p.ax
# 正電荷の位置に青丸と +
ax.scatter(0, 1, marker='o',
# size, color
s = 200, c = 'white',
edgecolors = 'blue', linewidths = 2)
ax.text(0, 1, "＋", color = "blue",
ha = "center", va = "center", size=12)
# 負電荷の位置に赤丸と -
ax.scatter(0, -1, marker='o',
# size, color
s = 200, c = 'white',
edgecolors = 'red', linewidths = 2)
ax.text(0, -1, "ー", color = "red",
ha = "center", va = "center", size=12)
# 枠と縦軸目盛を非表示に
ax.axis('off');
電気力線を描く
streamlines = True とすると，電場ベクトル $\boldsymbol{E}$ を接ベクトルとする電気力線を描くことができる。
plot_vector(
[E2x(x, y), E2y(x, y)],
(x, -5, 5), (y, -5, 5),
# 電場の流線すなわち電気力線を描く
streamlines = True,
scalar = False, grid = False, normalize=True,
use_cm = False,
size = (6.4, 6.4)
);
オプションを指定して描く例。電気力線は切れたりしないので stream_kw に "broken_streamlines":False とし，あまり混雑しないように "density" の値を 1 より小さくしてみる。また，電気力線を「先に」描いたあとに点電荷の位置に○を描くようにするため，"zorder":0 としてみる。
stream_kw に設定可能なオプションは，matplotlib.axes.Axes.streamplot() のオプションと同じ。以下を参照。
p = plot_vector(
[E2x(x, y), E2y(x, y)],
(x, -5, 5), (y, -5, 5),
streamlines = True,
scalar = False, grid = False, normalize=True,
use_cm = False, show = False,
size = (6.4, 6.4),
stream_kw = {
"linewidth":1.5, "zorder":0,
"density":0.3,
"color":"tab:blue",
"broken_streamlines":False
}
)
ax = p.ax
# 正電荷の位置に青丸と +
ax.scatter(0, 1, marker='o',
# size, color
s = 200, c = 'white',
edgecolors = 'blue', linewidths = 2)
ax.text(0, 1, "＋", color = "blue",
ha = "center", va = "center", size=12)
# 負電荷の位置に赤丸と -
ax.scatter(0, -1, marker='o',
# size, color
s = 200, c = 'white',
edgecolors = 'red', linewidths = 2)
ax.text(0, -1, "ー", color = "red",
ha = "center", va = "center", size=12)
# 枠と縦軸目盛を非表示に
ax.axis('off');
等電位面（線）を描く
正負の点電荷がつくる静電ポテンシャルは（適宜無次元化を行って）
$$\phi(x, y) = -\frac{1}{\sqrt{x^2 + (y-1)^2}}+\frac{1}{\sqrt{x^2 + (y+1)^2}}$$
電場ベクトル $\boldsymbol{E}$ は
$$\boldsymbol{E} = – \nabla \phi$$
であるから，$\phi = \mbox{const.}$ の等電位面（線）に対して電場ベクトル $\boldsymbol{E}$ は直交する。
def phi(x, y):
return -1/sqrt(x**2 + (y-1)**2) + 1/sqrt(x**2 + (y+1)**2)
これまでは scalar = False としていたのでスカラー場の contour は描いていなかったが，scalar = True とすればベクトル場の大きさをスカラー場として描く。
以下の例のように scalar = phi(x, y) と指定すればスカラー関数 phi(x, y) の contour を描く。上で定義された $\phi(x, y)$ は $-\infty$ から $+\infty$ までとりうるので，等電位面（線）を描くレベルを quiver_kw の levels に指定してみる。
静電ポテンシャル $\phi(x,y)$ の等電位面（線）に対して電場ベクトルが直交していることがわかる。
p = plot_vector(
[E2x(x, y), E2y(x, y)],
(x, -5.5, 5.5), (y, -5.5, 5.5),
n = 12, grid = False,
scalar = phi(x, y), use_cm = False, normalize=True,
xlim = (-6.5, 6.5), ylim = (-6.5, 6.5),
size = (6.4, 6.4), show = False,
# ベクトル場の色・大きさ・太さ
quiver_kw = {
"color":"tab:blue", "pivot":"mid",
"units":"xy", "scale":1.5, "width":0.07
},
# 等電位線
# 塗りつぶし無し，contour にラベル無し
is_filled=False, clabels = False,
contour_kw = {
"levels":[-1, -0.3, -0.1, -0.05,
0.05, 0.1, 0.3, 1],
# 等電位線は全て同色に
"cmap":None,
"colors":['lightgreen']*8
}
)
ax = p.ax
# 正電荷の位置に青丸と +
ax.scatter(0, 1, marker='o',
# size, color
s = 200, c = 'white',
edgecolors = 'blue', linewidths = 2)
ax.text(0, 1, "＋", color = "blue",
ha = "center", va = "center", size=12)
# 負電荷の位置に赤丸と -
ax.scatter(0, -1, marker='o',
# size, color
s = 200, c = 'white',
edgecolors = 'red', linewidths = 2)
ax.text(0, -1, "ー", color = "red",
ha = "center", va = "center", size=12)
# 枠と縦軸目盛を非表示に
ax.axis('off');
電気力線と等電位面（線）を両方描く
電気力線と等電位面（線）を両方描いてみると，確かに電気力線は常に等電位面（線）に直交していることが，よりはっきりとわかる。
p = plot_vector(
[E2x(x, y)/E(x,y), E2y(x, y)/E(x,y)],
(x, -5, 5), (y, -5, 5),
# 電場の流線すなわち電気力線
streamlines = True,
# scalar ここでは静電ポテンシャル
scalar = phi(x, y),
grid = False,
use_cm = False, show = False,
size = (6.4, 6.4),
is_filled=False, clabels = False,
contour_kw = {
"levels":[-1, -0.3, -0.1, -0.05,
0.05, 0.1, 0.3, 1],
"cmap":None,
"colors":['lightgreen']*8
},
stream_kw = {
"linewidth":1.5, "zorder":0,
"density":0.3,
"color":"tab:blue",
"broken_streamlines":False
}
)
ax = p.ax
# 正電荷の位置に青丸と +
ax.scatter(0, 1, marker='o',
# size, color
s = 200, c = 'white',
edgecolors = 'blue', linewidths = 2)
ax.text(0, 1, "＋", color = "blue",
ha = "center", va = "center", size=12)
# 負電荷の位置に赤丸と -
ax.scatter(0, -1, marker='o',
# size, color
s = 200, c = 'white',
edgecolors = 'red', linewidths = 2)
ax.text(0, -1, "ー", color = "red",
ha = "center", va = "center", size=12)
# 枠と縦軸目盛を非表示に
ax.axis('off');
