p3_tlt

点電荷による電気力線

点電荷による電界，電位

　クーロンの法則は，

N

個の点電荷

Q_{i}

[C]が位置ベクトル

r_{i}

[m]（

i = 1, 2, \dots, N

）にあるとき，観測点

r

[m]の静電界

E

[N/C]は，

\begin{matrix} (1) & E = \sum_{i = 1}^{N} E_{i}, E_{i} = \frac{1}{4 π ϵ_{0}} \frac{Q_{i} (r - r_{i})}{| r - r_{i} |^{3}} \end{matrix}

ここで，

ϵ_{0}

は真空中の誘電率を示し，

\begin{matrix} (2) & \frac{1}{4 π ϵ_{0}} = 9 \times 10^{9}, ϵ_{0} = 8.842 \times 10^{- 12} \end{matrix}

また，電位

V

は，

\begin{matrix} (3) & V = \frac{1}{4 π ϵ_{0}} \sum_{i = 1}^{N} \frac{Q_{i}}{R_{i}}, R_{i} = | r - r_{i} | \end{matrix}

モジュールのインポート

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')

電気力線の方程式

　電界

E

を，

\begin{array}{rcl} E & = & E a \\ (4) & = & E_{x} a_{x} + E_{y} a_{y} + E_{z} a_{z} \end{array}

より，電界に沿う単位ベクトル

a

は，

\begin{matrix} (5) & a = \frac{E_{x}}{E} a_{x} + \frac{E_{y}}{E} a_{y} + \frac{E_{y}}{E} a_{z} \end{matrix}

ただし，

a_{x}

，

a_{y}

，

a_{z}

は

x, y, z

方向の単位ベクトルを示す．この電界に沿うベクトル線要素

d s

を，

\begin{array}{rcl} d s & = & a d s \\ (6) & = & a_{x} d x + a_{y} d y + a_{z} d z \end{array}

とおくと，単位ベクトル

a

の成分比較より次のような常微分方程式が得られる．

\begin{array}{rcl} (7) & f_{1} (s, r) = x^{'} = \frac{d x}{d s} & = & \frac{E_{x}}{E} = a \cdot a_{x} \\ (8) & f_{2} (s, r) = y^{'} = \frac{d y}{d s} & = & \frac{E_{y}}{E} = a \cdot a_{y} \\ (9) & f_{3} (s, r) = z^{'} = \frac{d z}{d s} & = & \frac{E_{z}}{E} = a \cdot a_{z} \end{array}

ここで，

\begin{matrix} (10) & r = x a_{x} + y a_{y} + z a_{z} \end{matrix}

微分方程式

f_{i} (s, 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{array}{rcl} (11) & x (s) & = & x_{0} + \int_{s_{0}}^{s} \frac{E_{x} (s^{'})}{E (s^{'})} d s^{'} \\ (12) & y (s) & = & y_{0} + \int_{s_{0}}^{s} \frac{E_{y} (s^{'})}{E (s^{'})} d s^{'} \\ (13) & z (s) & = & z_{0} + \int_{s_{0}}^{s} \frac{E_{z} (s^{'})}{E (s^{'})} d s^{'} \end{array}

変数

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')