点電荷による電気力線

点電荷による電界,電位

 クーロンの法則は,N個の点電荷Qi [C]が位置ベクトル ri [m](i=1,2,,N)にあるとき,観測点 r [m]の静電界E [N/C]は, (1)E=i=1NEi,     Ei=14πϵ0Qi(rri)|rri|3 ここで,ϵ0は真空中の誘電率を示し, (2)14πϵ0=9×109,    ϵ0=8.842×1012 また,電位Vは, (3)V=14πϵ0i=1NQiRi,     Ri=|rri|

モジュールのインポート

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を, E=Ea(4)=Exax+Eyay+Ezaz より,電界に沿う単位ベクトルaは, (5)a=ExEax+EyEay+EyEaz ただし, axayazx,y,z方向の単位ベクトルを示す. この電界に沿うベクトル線要素dsを, ds=ads(6)=axdx+aydy+azdz とおくと,単位ベクトルaの成分比較より次のような常微分方程式が得られる. (7)f1(s,r)=x=dxds=ExE=aax(8)f2(s,r)=y=dyds=EyE=aay(9)f3(s,r)=z=dzds=EzE=aaz ここで, (10)r=xax+yay+zaz 微分方程式 fi(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=s0のとき,x=x0(s0), y=y0(s0)z=z0(s0)として, 積分形で表せば次のようになる. (11)x(s)=x0+s0sEx(s)E(s)ds(12)y(s)=y0+s0sEy(s)E(s)ds(13)z(s)=z0+s0sEz(s)E(s)ds 変数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')