import numpy as np 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'])
def p1charge(r, rd, q):
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
これより，複数の点電荷による電界，電位は，
def pncharge(r, rdn, q):
rd = rdn[0].reshape(1, 1, -1)
e, v = p1charge(r, rd, q[0])
if len(q)>=2:
for i in range(len(q)-1):
rd = rdn[i+1].reshape(1, 1, -1)
ei, vi = p1charge(r, rd, q[i+1])
e = e + ei
v = v + vi
return e,v
2次元分布の値をフォーマット付きで出力するため，
def table(sh,sv,h,v,f): # 数表の作成（フォーマットで整頓した出力）
print(f"{sh:>6s}",end=' ')
nh = len(h)
nv = len(v)
[print(f"{h[j]:6.3f}",end=' ') for j in range(nh)]
print('')
print(f"{sv:>6s}",end=' ')
print('')
for i in range(nv):
print(f'{v[i]:6.3f}',end=' ')
[print(f"{f[i][j]:6.3f}",end=' ') for j in range(nh)]
print('')
return
# Fig.1.3(b) に対応 q1 = np.array([1.0, -1.0]) # 点電荷の電荷量 rd1_p = np.array([[0.0, 1.0, 0.0]]) # 正の点電荷の位置ベクトル rd1_m = np.array([[0.0, -1.0, 0.0]]) # 負の点電荷の位置ベクトル rd1 = np.concatenate([rd1_p, rd1_m]) q1, rd1
(array([ 1., -1.]),
array([[ 0., 1., 0.],
[ 0., -1., 0.]]))
2次元分布を求めるため，メッシュグリッドのデータを求めておくと，
# 直角座標系(x,y)で表した観測点 x = np.linspace(-3, 3, 13) y = np.linspace(-3, 3, 13) xx, yy = np.meshgrid(x,y, indexing='xy')# xy座標の並び (デフォルト) zz = np.zeros(xx.shape) r = np.stack([xx,yy,zz], axis=-1)これより，電界および電位を求めると，
r00 = np.array([0.0, 0.0, 0.0]) e00, vv00 = pncharge(r00, rd1, q1) r0 = np.array([0.0, 0.9, 0.0]) e0, vv0 = pncharge(r0, rd1, q1) e1, vv1 = pncharge(r, rd1, q1) e_x = e1[:,:,0]/e00[0,0,1] e_y = e1[:,:,1]/e00[0,0,1] vvv = vv1[:,:,0]/vv0[0:,0:,0]ただし，電位は電荷の周辺のある点で規格化している．また，電界は原点の電界の振幅で値で規格化している． これより，電界の$x$成分$e_x$は，
table('x','y',x,y,e_x)
x -3.000 -2.500 -2.000 -1.500 -1.000 -0.500 0.000 0.500 1.000 1.500 2.000 2.500 3.000
y
-3.000 -0.020 -0.026 -0.033 -0.038 -0.038 -0.025 -0.000 0.025 0.038 0.038 0.033 0.026 0.020
-2.500 -0.024 -0.035 -0.049 -0.065 -0.075 -0.058 -0.000 0.058 0.075 0.065 0.049 0.035 0.024
-2.000 -0.028 -0.043 -0.068 -0.108 -0.161 -0.170 -0.000 0.170 0.161 0.108 0.068 0.043 0.028
-1.500 -0.028 -0.047 -0.084 -0.159 -0.332 -0.692 -0.000 0.692 0.332 0.159 0.084 0.047 0.028
-1.000 -0.024 -0.042 -0.081 -0.174 -0.455 -1.971 nan 1.971 0.455 0.174 0.081 0.042 0.024
-0.500 -0.014 -0.025 -0.050 -0.111 -0.272 -0.644 -0.000 0.644 0.272 0.111 0.050 0.025 0.014
0.000 -0.000 -0.000 -0.000 -0.000 -0.000 -0.000 -0.000 -0.000 -0.000 -0.000 -0.000 -0.000 -0.000
0.500 0.014 0.025 0.050 0.111 0.272 0.644 -0.000 -0.644 -0.272 -0.111 -0.050 -0.025 -0.014
1.000 0.024 0.042 0.081 0.174 0.455 1.971 nan -1.971 -0.455 -0.174 -0.081 -0.042 -0.024
1.500 0.028 0.047 0.084 0.159 0.332 0.692 -0.000 -0.692 -0.332 -0.159 -0.084 -0.047 -0.028
2.000 0.028 0.043 0.068 0.108 0.161 0.170 -0.000 -0.170 -0.161 -0.108 -0.068 -0.043 -0.028
2.500 0.024 0.035 0.049 0.065 0.075 0.058 -0.000 -0.058 -0.075 -0.065 -0.049 -0.035 -0.024
3.000 0.020 0.026 0.033 0.038 0.038 0.025 -0.000 -0.025 -0.038 -0.038 -0.033 -0.026 -0.020
同様にして，電界の$y$成分$e_y$は，
table('x','y',x,y,e_y)
x -3.000 -2.500 -2.000 -1.500 -1.000 -0.500 0.000 0.500 1.000 1.500 2.000 2.500 3.000
y
-3.000 -0.005 -0.011 -0.022 -0.038 -0.061 -0.084 -0.094 -0.084 -0.061 -0.038 -0.022 -0.011 -0.005
-2.500 -0.002 -0.008 -0.021 -0.047 -0.092 -0.150 -0.181 -0.150 -0.092 -0.047 -0.021 -0.008 -0.002
-2.000 0.004 -0.000 -0.013 -0.046 -0.129 -0.304 -0.444 -0.304 -0.129 -0.046 -0.013 -0.000 0.004
-1.500 0.012 0.013 0.010 -0.013 -0.115 -0.632 -1.920 -0.632 -0.115 -0.013 0.010 0.013 0.012
-1.000 0.021 0.030 0.044 0.064 0.089 0.114 nan 0.114 0.089 0.064 0.044 0.030 0.021
-0.500 0.029 0.045 0.077 0.142 0.307 0.897 2.222 0.897 0.307 0.142 0.077 0.045 0.029
0.000 0.032 0.051 0.089 0.171 0.354 0.716 1.000 0.716 0.354 0.171 0.089 0.051 0.032
0.500 0.029 0.045 0.077 0.142 0.307 0.897 2.222 0.897 0.307 0.142 0.077 0.045 0.029
1.000 0.021 0.030 0.044 0.064 0.089 0.114 nan 0.114 0.089 0.064 0.044 0.030 0.021
1.500 0.012 0.013 0.010 -0.013 -0.115 -0.632 -1.920 -0.632 -0.115 -0.013 0.010 0.013 0.012
2.000 0.004 -0.000 -0.013 -0.046 -0.129 -0.304 -0.444 -0.304 -0.129 -0.046 -0.013 -0.000 0.004
2.500 -0.002 -0.008 -0.021 -0.047 -0.092 -0.150 -0.181 -0.150 -0.092 -0.047 -0.021 -0.008 -0.002
3.000 -0.005 -0.011 -0.022 -0.038 -0.061 -0.084 -0.094 -0.084 -0.061 -0.038 -0.022 -0.011 -0.005
また，電位は，
table('x','y',x,y,vvv)
x -3.000 -2.500 -2.000 -1.500 -1.000 -0.500 0.000 0.500 1.000 1.500 2.000 2.500 3.000
y
-3.000 -0.008 -0.011 -0.014 -0.018 -0.022 -0.025 -0.026 -0.025 -0.022 -0.018 -0.014 -0.011 -0.008
-2.500 -0.009 -0.012 -0.016 -0.022 -0.030 -0.037 -0.040 -0.037 -0.030 -0.022 -0.016 -0.012 -0.009
-2.000 -0.008 -0.012 -0.018 -0.027 -0.041 -0.060 -0.070 -0.060 -0.041 -0.027 -0.018 -0.012 -0.008
-1.500 -0.008 -0.012 -0.018 -0.031 -0.055 -0.108 -0.169 -0.108 -0.055 -0.031 -0.018 -0.012 -0.008
-1.000 -0.006 -0.009 -0.015 -0.028 -0.058 -0.160 -inf -0.160 -0.058 -0.028 -0.015 -0.009 -0.006
-0.500 -0.003 -0.005 -0.009 -0.017 -0.036 -0.083 -0.141 -0.083 -0.036 -0.017 -0.009 -0.005 -0.003
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
0.500 0.003 0.005 0.009 0.017 0.036 0.083 0.141 0.083 0.036 0.017 0.009 0.005 0.003
1.000 0.006 0.009 0.015 0.028 0.058 0.160 inf 0.160 0.058 0.028 0.015 0.009 0.006
1.500 0.008 0.012 0.018 0.031 0.055 0.108 0.169 0.108 0.055 0.031 0.018 0.012 0.008
2.000 0.008 0.012 0.018 0.027 0.041 0.060 0.070 0.060 0.041 0.027 0.018 0.012 0.008
2.500 0.009 0.012 0.016 0.022 0.030 0.037 0.040 0.037 0.030 0.022 0.016 0.012 0.009
3.000 0.008 0.011 0.014 0.018 0.022 0.025 0.026 0.025 0.022 0.018 0.014 0.011 0.008
次に，電位の等高線図を描くために，サンプル点数を増やして，
# 直角座標系(x,y)で表した観測点 nxy = 201 x = np.linspace(-3, 3, nxy) y = np.linspace(-3, 3, nxy) xx, yy = np.meshgrid(x,y, indexing='xy')# xy座標の並び (デフォルト) zz = np.zeros(xx.shape) r = np.stack([xx,yy,zz], axis=-1)電位，電界を求めると，
r0 = np.array([0.0, 0.9, 0.0]) e0, vv0 = pncharge(r0, rd1, q1) e1, vv1 = pncharge(r, rd1, q1) e_x = e1[:,:,0] e_y = e1[:,:,1] vvv = vv1[:,:,0]/vv0[0:,0:,0]これより，プロットして，
fig = plt.figure(figsize=(8, 8)) # グラフ領域の作成
v = np.linspace(-1, 1, 21)
plt.streamplot(xx, yy, e_x, e_y, color='orange', density=1.0, linewidth=1.5)
ax = plt.contour(x, y, vvv, levels=v, linewidths=0.8)# 等高線表示
ax.clabel(fmt='%1.1f', inline=1, fontsize=7)# 等高線の値を表示
plt.contourf(x, y, vvv, levels=v, cmap="bwr")# twilight_shifted
plt.axis('scaled')
plt.colorbar(pad=0.1, shrink=0.77, orientation="vertical")
plt.plot(rd1_m[:,0], rd1_m[:,1],'o', color = 'green')
plt.plot(rd1_p[:,0], rd1_p[:,1],'o', color = 'red')
fig.savefig('e_v_fig_1_3_b.pdf')
plt.show()
fig = plt.figure() # グラフ領域の作成
plt.xlabel("$x, y$") # x軸のラベル設定
plt.ylabel("V") # y軸のラベル設定
plt.ylim(-1.5, 1.5) # y軸範囲の設定
plt.plot(r[:,nxy//2,1], vvv[:,nxy//2], label='$x=0$')
plt.plot(r[nxy//2,:,0], vvv[nxy//2,:], label='$y=0$')
plt.legend(ncol=1, loc='lower right', fancybox=False, frameon=True)
fig.savefig('e_v_fig_1_3_b_xy.pdf')
plt.show()
# Fig.1.4(a) に対応 q1 = np.array([1.0, 1.0, 1.0, -1.0, -1.0, -1.0]) # 正の点電荷の位置ベクトル rd1_p = np.array([[-1.0, 1.0, 0.0],[0.0, 1.0, 0.0],[1.0, 1.0, 0.0]]) # 負の点電荷の位置ベクトル rd1_m = np.array([[-1.0, -1.0, 0.0],[0.0, -1.0, 0.0],[1.0, -1.0, 0.0]]) rd1 = np.concatenate([rd1_p, rd1_m]) rd1
array([[-1., 1., 0.],
[ 0., 1., 0.],
[ 1., 1., 0.],
[-1., -1., 0.],
[ 0., -1., 0.],
[ 1., -1., 0.]])
電位，電界を求めると，
r0 = np.array([0.0, 0.9, 0.0]) e0, vv0 = pncharge(r0, rd1, q1) e1, vv1 = pncharge(r, rd1, q1) vvv = vv1[:,:,0]/vv0[0:,0:,0] e_x = e1[:,:,0] e_y = e1[:,:,1]これより，プロットして，
fig = plt.figure(figsize=(8, 8)) # グラフ領域の作成
v = np.linspace(-1.0, 1.0, 21)
plt.streamplot(xx, yy, e_x, e_y, color='orange', density=1.0, linewidth=1.5)
ax = plt.contour(x, y, vvv, levels=v, linewidths=0.8)# 等高線表示
ax.clabel(fmt='%1.1f', inline=1, fontsize=7)# 等高線の値を表示
plt.contourf(x, y, vvv, levels=v, cmap="bwr")# twilight_shifted
plt.axis('scaled')
plt.colorbar(pad=0.1, shrink=0.77, orientation="vertical")
plt.plot(rd1_m[:,0], rd1_m[:,1],'o', color = 'green')
plt.plot(rd1_p[:,0], rd1_p[:,1],'o', color = 'red')
fig.savefig('e_v_fig_1_4_a.pdf')
# Fig.1.4(b) に対応 q1 = np.array([1.0, 1.0, -1.0, -1.0]) rd1_p = np.array([[2.0, 0.0, 0.0],[1.0, 0.0, 0.0]]) rd1_m = np.array([[-1.0, 0.0, 0.0],[-2.0, 0.0, 0.0]]) rd1 = np.concatenate([rd1_p, rd1_m]) rd1
array([[ 2., 0., 0.],
[ 1., 0., 0.],
[-1., 0., 0.],
[-2., 0., 0.]])
電位，電界は，
r0 = np.array([-0.9, 0.0, 0.0]) e0, vv0 = pncharge(r0, rd1, q1) e1, vv1 = pncharge(r, rd1, q1) vvv = vv1[:,:,0]/vv0[0:,0:,0] e_x = e1[:,:,0] e_y = e1[:,:,1]プロットして，
fig = plt.figure(figsize=(8, 8)) # グラフ領域の作成
v = np.linspace(-1.0, 1.0, 21)
plt.streamplot(xx, yy, e_x, e_y, color='orange', density=1.0, linewidth=1.5)
ax = plt.contour(x, y, vvv, levels=v, linewidths=0.8)# 等高線表示
ax.clabel(fmt='%1.1f', inline=1, fontsize=7)# 等高線の値を表示
plt.contourf(x, y, vvv, levels=v, cmap="bwr")# twilight_shifted
plt.axis('scaled')
plt.colorbar(pad=0.1, shrink=0.77, orientation="vertical")
plt.plot(rd1_m[:,0], rd1_m[:,1],'o', color = 'green')
plt.plot(rd1_p[:,0], rd1_p[:,1],'o', color = 'red')
fig.savefig('e_v_fig_1_4_b.pdf')
【例題】正方形の各頂点に$Q$[C]の点電荷がおいてあるとき，その中心に何クーロンの点電荷をおけば各電荷 （各頂点の全電荷および中心の電荷）はつりあうか．
略解 中心に何も置かなければ，各頂点の電荷が同符号であるから反発力が生じ，つりあわないことは明らか． いま，中心におく電荷を $q$ [C]とし，中心を座標原点，各頂点を$x$軸および$y$軸上となるよう$xy$座標系を定義する． 各頂点の座標はA$(a,0)$，B$(0,a)$，C$(-a,0)$，D$(0,-a)$である．このとき，点Aでのクーロン力 $\VECi{F}_A$は，(C-A), (B-A), (D-A), (O-A)の順に，クーロンの法則より次のようになる． \begin{gather} \VECi{F}_A = \frac{1}{4\pi \epsilon _0} \left\{ \frac{Q^2}{(2a)^2} \VECi{i} + \frac{Q^2}{l^2} \frac{(\VECi{i} - \VECi{j})}{\sqrt{2}} + \frac{Q^2}{l^2} \frac{(\VECi{i} + \VECi{j})}{\sqrt{2}} + \frac{Q q}{a^2} \VECi{i} \right\} \end{gather} ここで，$l$は正方形の一辺の長さを示す． \begin{gather} l^2 = a^2 + a^2 = 2 a^2 \end{gather} $\VECi{F}_A$の$\VECi{j}$成分はキャンセルして，$\VECi{i}$成分だけが残り，次のようになる． \begin{eqnarray} \VECi{F}_A &=& \frac{Q}{4\pi \epsilon _0} \left\{ \frac{Q}{2 l^2} + \frac{Q}{l^2} \frac{1}{\sqrt{2}} + \frac{Q}{l^2} \frac{1}{\sqrt{2}} + \frac{2q}{l^2} \right\} \VECi{i} \nonumber \\ &=& \frac{Q}{4\pi \epsilon _0 l^2} \left\{ \left( \frac{1}{2} + \frac{\sqrt{2}}{2} + \frac{\sqrt{2}}{2} \right) Q + 2q \right\} \VECi{i} \nonumber \\ &=& \frac{Q}{4\pi \epsilon _0 l^2} \left( \frac{1+2\sqrt{2}}{2} Q + 2q \right) \VECi{i} \end{eqnarray} つりあうためには，$\VECi{F}_A = 0$とすればよいから， \begin{gather} \frac{1+2\sqrt{2}}{2} Q + 2q = 0 \end{gather} よって，電荷$q$ [C]は次のようになる． \begin{gather} q = -\frac{1+2\sqrt{2}}{4} Q \ \ \ \mbox{[C]} \end{gather} 原点のまわりに$90^\circ$ 回転させても同じ結果が得られ， 他の点B，C，D全て同じ結果が得られることは明らか． このとき，中心にある電荷に働く力は，対称性より電荷$q$の値に依らず常につりあう． したがって，全電荷がつりあうのは，中心に$-0.957 Q$ [C]の点電荷をおいたときである． なお，得られた結果は距離には依らない条件である． 規格化した電荷$q/Q$は， \begin{gather} \frac{q}{Q} = -\frac{1+2\sqrt{2}}{4}\simeq -0.957 \end{gather}q0 = -(1+2*np.sqrt(2))/4 q0
-0.9571067811865476
# 電荷に及ぼすクーロン力の釣り合い q1 = np.array([1.0, 1.0, 1.0, 1.0, q0]) # 正の点電荷の位置ベクトル rd1_p = np.array([[1.0, 0.0, 0.0],[0.0, 1.0, 0.0],[-1.0, 0.0, 0.0], [0.0, -1.0, 0.0]]) # 負の点電荷の位置ベクトル rd1_m = np.array([[0.0, 0.0, 0.0]]) rd1 = np.concatenate([rd1_p, rd1_m]) rd1
array([[ 1., 0., 0.],
[ 0., 1., 0.],
[-1., 0., 0.],
[ 0., -1., 0.],
[ 0., 0., 0.]])
r0 = np.array([0.9, 0.0, 0.0]) e0, vv0 = pncharge(r0, rd1, q1) e1, vv1 = pncharge(r, rd1, q1) vvv = vv1[:,:,0]/vv0[0:,0:,0] e_x = e1[:,:,0] e_y = e1[:,:,1]
fig = plt.figure(figsize=(8, 8)) # グラフ領域の作成
v = np.linspace(-1.0, 1.0, 21)
plt.streamplot(xx, yy, e_x, e_y, color='orange', density=1.0, linewidth=1.5)
ax = plt.contour(x, y, vvv, levels=v, linewidths=0.8)# 等高線表示
ax.clabel(fmt='%1.1f', inline=1, fontsize=7)# 等高線の値を表示
plt.contourf(x, y, vvv, levels=v, cmap="bwr")# twilight_shifted
plt.axis('scaled')
plt.colorbar(pad=0.1, shrink=0.77, orientation="vertical")
plt.plot(rd1_m[:,0], rd1_m[:,1],'o', color = 'green')
plt.plot(rd1_p[:,0], rd1_p[:,1],'o', color = 'red')
fig.savefig('e_v_fig_p4_m1.pdf')
fig = plt.figure() # グラフ領域の作成
plt.xlabel("$x, y$") # x軸のラベル設定
plt.ylabel("V") # y軸のラベル設定
plt.plot(r[:,nxy//2,1], vvv[:,nxy//2], label='$x=0$')
plt.plot(r[nxy//2,:,0], vvv[nxy//2,:], label='$y=0$')
plt.legend(ncol=1, loc='lower right', fancybox=False, frameon=True)
fig.savefig('e_v_fig_p4_m1_xy.pdf')