点電荷による電気力線
点電荷による電界，電位
クーロンの法則は，$N$個の点電荷$Q_i$ [C]が位置ベクトル
$\boldsymbol{r}_i$ [m]（$i=1,2, \cdots, N$）にあるとき，観測点
$\boldsymbol{r}$ [m]の静電界$\boldsymbol{E}$ [N/C]は，
\begin{gather}
\boldsymbol{E} = \sum_{i=1}^N \boldsymbol{E}_i, \ \ \ \ \
\boldsymbol{E}_i = \frac{1}{4 \pi \epsilon_0} \frac{Q_i(\boldsymbol{r}-\boldsymbol{r}_i)}{|\boldsymbol{r}-\boldsymbol{r}_i|^3}
\end{gather}
ここで，$\epsilon_0$は真空中の誘電率を示し，
\begin{gather}
\frac{1}{4 \pi \epsilon_0} = 9 \times 10^9, \ \ \ \
\epsilon_0 = 8.842 \times 10^{-12}
\end{gather}
また，電位$V$は，
\begin{gather}
V = \frac{1}{4 \pi \epsilon_0} \sum_{i=1}^N \frac{Q_i}{R_i}, \ \ \ \ \
R_i = |\boldsymbol{r}-\boldsymbol{r}_i|
\end{gather}
モジュールのインポート
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import scienceplots
#Warning : As of version 2.0.0, you need to add import scienceplots before setting the style (plt.style.use('science')).
plt.style.use(['science', 'notebook'])
ユーザ関数
位置ベクトルより，単位ベクトル，2点間の距離を求める関数は，
def uvector(aa): # 多次元ベクトルの振幅，単位ベクトル
norm = np.linalg.norm(aa)
uv = aa / norm
return norm, uv
1つの点電荷による電界，電位は，
# r: 観測点での位置ベクトルの配列
# rd: 点電荷のある位置ベクトルの配列
# q: 点電荷の電荷量の配列
def p1charge_3d(r, rd, q): # 1つの点電荷による電界，電位
qk = q/4.0/np.pi
rr = r-rd
norm = np.linalg.norm(rr, axis=2, keepdims=True) # 観測点と点電荷との間の距離
v = qk / norm # 電位の配列
with np.errstate(divide="ignore", invalid="ignore"):
e = qk*rr/norm**3 # 電界の配列
return e, v, norm
複数の点電荷による電界，電位は，
def pncharge_3d(r, rdn, q): # 多数の点電荷による電界，電位
rd = rdn[0].reshape(1, 1, -1)
e, v, R = p1charge_3d(r, rd, q[0])
nq = len(q)
if nq>=2:
for i in range(nq-1):
rd = rdn[i+1].reshape(1, 1, -1)
ei, vi, Ri = p1charge_3d(r, rd, q[i+1])
e = e + ei
v = v + vi
R = np.append(R, Ri)
return e, v, R
テキストの図1.3(b)
テキストの図1.3(b)のように，同振幅・異符号の2つの点電荷のある場を考える．
# Fig.1.3(b) に対応
q = np.array([1.0, -1.0]) # 点電荷の電荷量
rd_p = np.array([[0.0, 1.0, 0.0]]) # 正の点電荷の位置ベクトル
rd_m = np.array([[0.0, -1.0, 0.0]]) # 負の点電荷の位置ベクトル
rd = np.concatenate([rd_p, rd_m])
rd
array([[ 0., 1., 0.],
[ 0., -1., 0.]])
点電荷の位置をプロットすると，
fig = plt.figure() # グラフ領域の作成
plt.axes().set_aspect('equal')
plt.xlim(-3.0, 3.0) # x軸範囲の設定
plt.ylim(-3.0, 3.0) # y軸範囲の設定
plt.xticks(np.linspace(-3.0, 3.0, 7))
plt.yticks(np.linspace(-3.0, 3.0, 7))
plt.xlabel("$x$") # 横軸のラベル設定
plt.ylabel("$y$") # 縦軸のラベル設定
plt.plot(rd_m[:,0], rd_m[:,1],'o')
plt.plot(rd_p[:,0], rd_p[:,1],'o')
fig.savefig('fig_1_3_b_none.eps')
電気力線の方程式
電界$\boldsymbol{E}$を，
\begin{eqnarray}
\boldsymbol{E}
&=& E \boldsymbol{a}
\nonumber \\
&=& E_x \boldsymbol{a}_x + E_y \boldsymbol{a}_y + E_z \boldsymbol{a}_z
\end{eqnarray}
より，電界に沿う単位ベクトル$\boldsymbol{a}$は，
\begin{gather}
\boldsymbol{a} = \frac{E_x}{E} \boldsymbol{a}_x
+ \frac{E_y}{E} \boldsymbol{a}_y + \frac{E_y}{E} \boldsymbol{a}_z
\end{gather}
ただし，
$\boldsymbol{a}_x$，$\boldsymbol{a}_y$，$\boldsymbol{a}_z$
は$x, y, z$方向の単位ベクトルを示す．
この電界に沿うベクトル線要素$d\boldsymbol{s}$を，
\begin{eqnarray}
d\boldsymbol{s}
&=& \boldsymbol{a} ds
\nonumber \\
&=& \boldsymbol{a}_x dx + \boldsymbol{a}_y dy + \boldsymbol{a}_z dz
\end{eqnarray}
とおくと，単位ベクトル$\boldsymbol{a}$の成分比較より次のような常微分方程式が得られる．
\begin{eqnarray}
f_1(s, \VEC{r}) = x' = \frac{dx}{ds} &=& \frac{E_x}{E} = \boldsymbol{a} \cdot \boldsymbol{a}_x
\\
f_2(s, \VEC{r}) = y' = \frac{dy}{ds} &=& \frac{E_y}{E} = \boldsymbol{a} \cdot \boldsymbol{a}_y
\\
f_3(s, \VEC{r}) = z' = \frac{dz}{ds} &=& \frac{E_z}{E} = \boldsymbol{a} \cdot \boldsymbol{a}_z
\end{eqnarray}
ここで，
\begin{gather}
\VECi{r} = x \VECi{a}_x + y \VECi{a}_y + z \VECi{a}_z
\end{gather}
微分方程式 $f_i(s,\VECi{r}) (i=1,2,3)$ は，
# q: 点電荷の電荷量のndarray
# rd: 点電荷の位置ベクトル(x,y)のndarray
def func_eline_3d(s, r, q, rd):
e, v, R = pncharge_3d(r, rd, q)
eamp, vs = uvector(e)
R_min = np.min(R)
if(R_min < 5.0e-4):
vs = vs*0.0
vs0 = vs[0,0,:]
return vs0
初期値$s=s_0$のとき，$x=x_0(s_0)$, $y=y_0(s_0)$，$z=z_0(s_0)$として，
積分形で表せば次のようになる．
\begin{eqnarray}
x(s) &=& x_0 + \int_{s_0}^s \frac{E_x(s')}{E(s')} ds'
\\
y(s) &=& y_0 + \int_{s_0}^s \frac{E_y(s')}{E(s')} ds'
\\
z(s) &=& z_0 + \int_{s_0}^s \frac{E_z(s')}{E(s')} ds'
\end{eqnarray}
変数$s$，および初期値の位置ベクトルを点電荷の周囲近傍とし，
s_ini, s_end, ns = 0.0, 5, 101
s = np.linspace(s_ini, s_end, ns)
nphi = 25
phi = np.linspace(0.0, 2.0*np.pi, nphi)
d_eps = 0.01
d_x, d_y = d_eps*np.cos(phi), d_eps*np.sin(phi)
Scipyの関数solve_ivpより，微分方程式を数値的に解くと，
def e_field_lines_2(q, rd, s, d_x, d_y, nphi):
nq = len(q)
ns = len(s)
s_ini = s[0]
x_end = s[-1]
ode = np.empty([ns,3,nphi*nq])
ij = 0
for j in range(nq):
y0 = rd[j,:]
for i in range(nphi):
yeps = np.array([d_x[i], d_y[i], 0.0])
y_ini = y0 + yeps
if q[j]>0: # プラスの電荷から計算した電気力線の軌跡
ode2 = solve_ivp(func_eline_3d, [s_ini, s_end], y_ini, \
dense_output=True, args=(q,rd,))
y = ode2.sol(s)
else: # マイナスの電荷から計算した電気力線の軌跡
ode2 = solve_ivp(func_eline_3d, [s_ini, -s_end], y_ini, \
dense_output=True, args=(q,rd,))
y = ode2.sol(-s)
yt = np.transpose(y)
plt.plot(yt[:,0], yt[:,1], linewidth=0.5, color="RebeccaPurple")
ode[:,:,ij] = yt
ij = ij+1
return ode
これより電気力線を描くと，
fig = plt.figure() # グラフ領域の作成
plt.axes().set_aspect('equal')
plt.xlim(-3.0, 3.0) # x軸範囲の設定
plt.ylim(-3.0, 3.0) # y軸範囲の設定
plt.xticks(np.linspace(-3.0, 3.0, 7))
plt.yticks(np.linspace(-3.0, 3.0, 7))
ode = e_field_lines_2(q, rd, s, d_x, d_y, nphi) # solve_ivp
plt.plot(rd_m[:,0], rd_m[:,1],'o')
plt.plot(rd_p[:,0], rd_p[:,1],'o')
fig.savefig('fig_1_3_b.pdf')
plt.show()
テキストの図1.4(a)
テキストの図1.4(a)のように，3つの正の点電荷と3つの負の点電荷による電界は，
q = np.array([1.0, 1.0, 1.0, -1.0, -1.0, -1.0])
# 正の点電荷の位置ベクトル
rd_p = np.array([[-1.0, 1.0, 0.0],[0.0, 1.0, 0.0],[1.0, 1.0, 0.0]])
# 負の点電荷の位置ベクトル
rd_m = np.array([[-1.0, -1.0, 0.0],[0.0, -1.0, 0.0],[1.0, -1.0, 0.0]])
rd = np.concatenate([rd_p, rd_m])
rd
array([[-1., 1., 0.],
[ 0., 1., 0.],
[ 1., 1., 0.],
[-1., -1., 0.],
[ 0., -1., 0.],
[ 1., -1., 0.]])
同様にして，微分方程式の数値計算より電気力線を描くと，
fig = plt.figure() # グラフ領域の作成
plt.axes().set_aspect('equal')
plt.xlim(-3.0, 3.0) # x軸範囲の設定
plt.ylim(-3.0, 3.0) # y軸範囲の設定
plt.xticks(np.linspace(-3.0, 3.0, 7))
plt.yticks(np.linspace(-3.0, 3.0, 7))
ode = e_field_lines_2(q, rd, s, d_x, d_y, nphi) # solve_ivp
plt.plot(rd_m[:,0], rd_m[:,1],'o')
plt.plot(rd_p[:,0], rd_p[:,1],'o')
fig.savefig('fig_1_4_a.pdf')
テキストの図1.4(b)
また，テキストの図1.4(b)のように，2つの正の点電荷と2つの負の点電荷による電界は，
q = np.array([1.0, 1.0, -1.0, -1.0])
rd_p = np.array([[2.0, 0.0, 0.0],[1.0, 0.0, 0.0]])
rd_m = np.array([[-1.0, 0.0, 0.0],[-2.0, 0.0, 0.0]])
rd = np.concatenate([rd_p, rd_m])
rd
array([[ 2., 0., 0.],
[ 1., 0., 0.],
[-1., 0., 0.],
[-2., 0., 0.]])
fig = plt.figure() # グラフ領域の作成
plt.axes().set_aspect('equal')
plt.xlim(-3.0, 3.0) # x軸範囲の設定
plt.ylim(-3.0, 3.0) # y軸範囲の設定
plt.xticks(np.linspace(-3.0, 3.0, 7))
plt.yticks(np.linspace(-3.0, 3.0, 7))
ode = e_field_lines_2(q, rd, s, d_x, d_y, nphi) # solve_ivp
plt.plot(rd_m[:,0], rd_m[:,1],'o')
plt.plot(rd_p[:,0], rd_p[:,1],'o')
fig.savefig('fig_1_4_b.pdf')
前のページに戻る
「目次」のページに戻る