Matplotlib で電場ベクトルの方向場を描く

正負の点電荷がつくる電場の向きを表す方向場を Python の Matplotlib で描く。

正負の点電荷がつくる電場

$\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}}$$

ライブラリの import

In [1]:
# NumPy 
import numpy as np
# Matplotlib 
import matplotlib.pyplot as plt

# グラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
# フォント設定
plt.rcParams['font.family'] = 'Noto Sans CJK JP'
plt.rcParams['font.size'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12

quiver() で電場の方向場を描く

電場ベクトルの定義

電場ベクトルを x, y の2変数関数として定義。

In [2]:
# 1個の点電荷がつくる電場
def Ex(x, y):
    return x/np.sqrt(x**2 + y**2)**3
def Ey(x, y):
    return y/np.sqrt(x**2 + y**2)**3

# 正負の2個の点電荷がつくる電場
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)

# 正負の2個の点電荷がつくる電場の大きさ
def E(x, y):
    return np.sqrt(E2x(x, y)**2 + E2y(x, y)**2)

# 規格化された電場ベクトル
def hE2x(x, y):
    return E2x(x, y)/E(x, y)
def hE2y(x, y):
    return E2y(x, y)/E(x, y)

ベクトル場を描く

点電荷が $y$ 軸上にあるので，リストを作成する際，ベクトルの始点が $y$ 軸に来ないようにしておく。（点電荷上では電場の分母がゼロになってしまうから。）
以下の例では，$x$ 座標は

$$-5.5, -4.5, \cdots, -0.5, 0.5, \cdots, 5.5$$

となり，$y$ 軸上 ($x = 0$) に始点がこない。

In [3]:
x = y = np.linspace(-5.5, 5.5, 12)
x, y = np.meshgrid(x, y)

plt.quiver() の書式は…

plt.quiver(始点のx座標, 始点のy座標, ベクトルのx成分, ベクトルのy成分)
In [4]:
# グラフのサイズ
fig = plt.figure(figsize=(6.4, 6.4))

# 表示範囲
plt.xlim(-6.5, 6.5)
plt.ylim(-6.5, 6.5)

# グリッドの表示
plt.grid()

# x の主目盛を 1 刻みに
plt.xticks(np.arange(-5, 6, 1))
# y の主目盛を 1 刻みに
plt.yticks(np.arange(-5, 6, 1))

# x, y を始点としたベクトル場の表示
plt.quiver(x, y, hE2x(x, y), hE2y(x, y))

# 念のために各ベクトルの始点を赤点で
plt.scatter(x, y, c='r', s=5);

pivot='mid' オプション

pivot='mid' オプションを使うと，ベクトルの真ん中を中心にしてベクトルの向きを変えて表示します。

In [5]:
# グラフのサイズ
fig = plt.figure(figsize=(6.4, 6.4))

# 表示範囲
plt.xlim(-6.5, 6.5)
plt.ylim(-6.5, 6.5)

# グリッドの表示
plt.grid()

# x の主目盛を 1 刻みに
plt.xticks(np.arange(-5, 6, 1))
# y の主目盛を 1 刻みに
plt.yticks(np.arange(-5, 6, 1))

# pivot = 'mid' としたベクトル場の表示
plt.quiver(x, y, hE2x(x, y), hE2y(x, y), pivot = 'mid')

# 念のために各ベクトルのpivotを赤点で
plt.scatter(x, y, c='r', s=5);

点電荷の位置にラベル

正電荷，負電荷の位置にラベルをつけてみる。

In [6]:
# グラフのサイズ
fig = plt.figure(figsize=(6.4, 6.4))

# 表示範囲
plt.xlim(-6.5, 6.5)
plt.ylim(-6.5, 6.5)

# グリッドの表示
plt.grid()

# x の主目盛を 1 刻みに
plt.xticks(np.arange(-5, 6, 1))
# y の主目盛を 1 刻みに
plt.yticks(np.arange(-5, 6, 1))

# pivot = 'mid' としたベクトル場の表示
plt.quiver(x, y, hE2x(x, y), hE2y(x, y), pivot = 'mid')

# 正電荷の位置に青丸と +
plt.scatter(0, 1, marker='o', 
            # size, color
            s = 200, c = 'white',
            edgecolors = 'blue', linewidths = 2)
plt.text(0, 1, "＋", color = "blue",
         horizontalalignment = "center", 
         verticalalignment = "center")

# 負電荷の位置に赤丸と -
plt.scatter(0, -1, marker='o', 
            # size, color
            s = 200, c = 'white',
            edgecolors = 'red', linewidths = 2)
plt.text(0, -1, "ー", color = "red",
         horizontalalignment = "center", 
         verticalalignment = "center");

Colormap (cmap)

plt.quiver() の書式で5番目の引数を設定すると，その大きさに応じたカラーマップが使える。

plt.quiver(x, y, vx, vy, C, cmap = ...)

以下の例では電場の大きさ E(x, y) に応じて矢印の色が変化する。
カラーマップ cmap = "Blues" を指定しているので，電場の強い（E(x, y) の値が大きい）ところは濃い青で，弱い（値が小さい）ところほど薄い青になる。

In [7]:
# グラフのサイズ
fig = plt.figure(figsize=(6.4, 6.4))

# 表示範囲
plt.xlim(-6.5, 6.5)
plt.ylim(-6.5, 6.5)

# グリッドの表示
plt.grid()

# x の主目盛を 1 刻みに
plt.xticks(np.arange(-5, 6, 1))
# y の主目盛を 1 刻みに
plt.yticks(np.arange(-5, 6, 1))

# 電場の大きさに応じたカラーマップの設定
plt.quiver(x, y, hE2x(x, y), hE2y(x, y), 
           E(x, y), cmap = "Blues", pivot = 'mid');

cmap の色味，ベクトルの大きさ・太さの設定

電荷の位置から離れると，電場の大きさは急激に小さくなるので，このままでは薄くてほとんど見えないベクトルばかりになる。

これだと少し寂しいのでカラーマップの範囲を設定する clim を調整して，離れた位置のベクトルもそこそこの青色になるように調整する。

上記の例では clim = (0, 2) くらいであろうから，ゼロに近い大きさでも色がつくように手探りで数値をいれてみます。

また，ベクトルそのものの大きさ・太さも少し変更してみる。書式はまず units='xy' などと設定後，ベクトルの長さを伸ばしたいなら scale により小さい値を，太くしたいなら width により大きい値をいれます。これも手探りでいれてみました。

In [8]:
# グラフのサイズ
fig = plt.figure(figsize=(6.4, 6.4))

# 表示範囲
plt.xlim(-6.5, 6.5)
plt.ylim(-6.5, 6.5)

# グリッドの表示
plt.grid()

# x の主目盛を 1 刻みに
plt.xticks(np.arange(-5, 6, 1))
# y の主目盛を 1 刻みに
plt.yticks(np.arange(-5, 6, 1))

# 
plt.quiver(x, y, hE2x(x, y), hE2y(x, y), 
           E(x, y), cmap = "Blues", pivot='mid', 
           # 数値は手探りで
           clim=(-0.1, 0.8), 
           # units 設定後，長さは scale, 太さは width で
           units='xy', scale = 1, width=0.12);

clim を手探りで設定する以外の方法として，以下の例では電場の大きさそのものではなく， 自然対数 np.log(E(x, y)) に応じて矢印の色が変化するようにしてみた。大きさが急激に変わるような場合には，大きさそのものよりも，np.log() や np.sqrt() などと使って変化を緩やかにするという手もあるかもしれない。

In [9]:
# グラフのサイズ
fig = plt.figure(figsize=(6.4, 6.4))

# 表示範囲
plt.xlim(-6.5, 6.5)
plt.ylim(-6.5, 6.5)

# グリッドの表示
plt.grid()

# x の主目盛を 1 刻みに
plt.xticks(np.arange(-5, 6, 1))
# y の主目盛を 1 刻みに
plt.yticks(np.arange(-5, 6, 1))

# 
plt.quiver(x, y, hE2x(x, y), hE2y(x, y), 
           # E(x, y) ではなく log(E(x, y)) にして
           # 大きさの変化を緩やかに
           np.log(E(x, y)), cmap = "Blues", pivot='mid', 
           # units 設定後，長さは scale, 太さは width で
           units='xy', scale = 1, width=0.12);

一応の完成作

以上の設定をまとめてみたのが以下の例。少し短めのベクトルをたくさん描いてみた。

In [10]:
x = y = np.linspace(-5.25, 5.25, 22)
x, y = np.meshgrid(x, y)

# グラフのサイズ
fig = plt.figure(figsize=(6.4, 6.4))

# 枠と縦軸目盛を非表示に
ax = fig.add_subplot()
ax.axis('off')

# 表示範囲
plt.xlim(-5, 5)
plt.ylim(-5, 5)

plt.quiver(x, y, hE2x(x, y), hE2y(x, y), 
           # E(x, y) ではなく log(E(x, y)) にして
           # 大きさの変化を緩やかに
           np.log(E(x, y)), cmap = "Blues", pivot='mid', 
           # units 設定後，長さは scale, 太さは width で
           units='xy', scale = 1.8, width=0.08)

# 正電荷の位置に青丸と +
plt.scatter(0, 1, marker='o', 
            # size, color
            s = 200, c = 'white',
            edgecolors = 'blue', linewidths = 2)
plt.text(0, 1, "＋", color = "blue",
         horizontalalignment = "center", 
         verticalalignment = "center")
# 負電荷の位置に赤丸と -
plt.scatter(0, -1, marker='o', 
            # size, color
            s = 200, c = 'white',
            edgecolors = 'red', linewidths = 2)
plt.text(0, -1, "ー", color = "red",
         horizontalalignment = "center", 
         verticalalignment = "center");

streamplot() で電気力線を描く

streamplot() を使うと，電場ベクトル $\boldsymbol{E}$ を接ベクトルとする電気力線を描くことができる。

以下の例では，単位ベクトル $\displaystyle \frac{\boldsymbol{E}}{\sqrt{\boldsymbol{E}\cdot\boldsymbol{E}}}$ を接ベクトルとしている。まずはオプションなしで…

In [11]:
# グラフのサイズ
fig = plt.figure(figsize=(6.4, 6.4))

# 枠と縦軸目盛を非表示に
ax = fig.add_subplot()
ax.axis('off')

plt.streamplot(x, y, hE2x(x, y), hE2y(x, y));

オプションを指定して描く例。電気力線は切れたりしないので broken_streamlines=False とし，あまり混雑しないように density の値を 1 より小さくしてみる。また，電気力線を「先に」描いたあとに点電荷の位置に○を描くようにするため，zorder = 0 としてみる。

In [12]:
# グラフのサイズ
fig = plt.figure(figsize=(6.4, 6.4))

# 枠と縦軸目盛を非表示に
ax = fig.add_subplot()
ax.axis('off')

plt.streamplot(x, y, hE2x(x, y), hE2y(x, y),
               # 各オプションの設定
               broken_streamlines = False, 
               density=0.3, 
               zorder = 0, 
               color = 'tab:blue');

# 正電荷の位置に青丸と +
plt.scatter(0, 1, marker='o', 
            # size, color
            s = 200, c = 'white',
            edgecolors = 'blue', linewidths = 2)
plt.text(0, 1, "＋", color = "blue",
         horizontalalignment = "center", 
         verticalalignment = "center")

# 負電荷の位置に赤丸と -
plt.scatter(0, -1, marker='o', 
            # size, color
            s = 200, c = 'white',
            edgecolors = 'red', linewidths = 2)
plt.text(0, -1, "ー", color = "red",
         horizontalalignment = "center", 
         verticalalignment = "center");

contour() で等電位面（線）を描く

正負の点電荷がつくる静電ポテンシャルは（適宜無次元化を行って）

$$\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}$ は直交する。

In [13]:
# 静電ポテンシャル
def phi(x, y):
    return -1/np.sqrt(x**2 + (y-1)**2) + 1/np.sqrt(x**2 + (y+1)**2)
In [14]:
# グラフのサイズ
fig = plt.figure(figsize=(6.4, 6.4))

# 枠と縦軸目盛を非表示に
ax = fig.add_subplot()
ax.axis('off')

# 電場ベクトルの始点
x = y = np.linspace(-5.5, 5.5, 12)
x, y = np.meshgrid(x, y)
# 電場ベクトルを描く
plt.quiver(x, y, hE2x(x, y), hE2y(x, y), 
           pivot = 'mid', color='tab:blue')

# contour の meshgrid
delta = 0.1
xrange = np.arange(-6, 6, delta)
yrange = np.arange(-6, 6, delta)
x, y = np.meshgrid(xrange, yrange)
# 等電位面（線）を描く
plt.contour(x, y, phi(x, y), 
            # 等電位線を描く phi(x, y) の値のリスト
            [-1,  -0.3, -0.1, -0.05, 
             0.05, 0.1, 0.3, 1], 
            # 等電位線の色，今回は全て同色に
            colors=['lightgreen']*8);

電気力線と等電位面（線）を両方描く

電気力線と等電位面（線）を両方描いてみると，確かに電気力線は常に等電位面（線）に直交していることが，よりはっきりとわかる。

In [15]:
# グラフのサイズ
fig = plt.figure(figsize=(6.4, 6.4))

# 枠と縦軸目盛を非表示に
ax = fig.add_subplot()
ax.axis('off')

# 電場ベクトルの始点
x = y = np.linspace(-5.5, 5.5, 12)
x, y = np.meshgrid(x, y)
# 電気力線を描く
plt.streamplot(x, y, hE2x(x, y), hE2y(x, y),
               broken_streamlines = False, 
               density=0.3, 
               zorder = 0, 
               color = 'tab:blue');

# 正電荷の位置に青丸と +
plt.scatter(0, 1, marker='o', 
            # size, color
            s = 200, c = 'white',
            edgecolors = 'blue', linewidths = 2)
plt.text(0, 1, "＋", color = "blue",
         horizontalalignment = "center", 
         verticalalignment = "center")

# 負電荷の位置に赤丸と -
plt.scatter(0, -1, marker='o', 
            # size, color
            s = 200, c = 'white',
            edgecolors = 'red', linewidths = 2)
plt.text(0, -1, "ー", color = "red",
         horizontalalignment = "center", 
         verticalalignment = "center");

# contour の meshgrid
delta = 0.1
xrange = np.arange(-6, 6, delta)
yrange = np.arange(-6, 6, delta)
x, y = np.meshgrid(xrange, yrange)
# 等電位面を描く
plt.contour(x, y, phi(x, y), 
            # 等電位線を描く phi(x, y) の値のリスト
            [-1,  -0.3, -0.1, -0.05, 
             0.05, 0.1, 0.3, 1], 
            # 等電位線の色，今回は全て同色に
            colors=['lightgreen']*8);