(初版:2020年3月、更新:2023年2月21日)
原点に配置した電荷 $Q$ による静電ポテンシャル $V$ は
$$\displaystyle{V = \frac{Q}{\sqrt{x^2 + y^2 + z^2}}}$$と書ける($4\pi\varepsilon_0 \to 1$とする)。電場 $\vec{E}$ は $\vec{E} = -\nabla V$ から
$$\begin{align*} \vec{E} &= \left( \frac{Q x}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{3}{2}}}, \frac{Q y}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{3}{2}}}, \frac{Q z}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{3}{2}}} \right)\\ &= \frac{Q}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{3}{2}}}\left(x, y, z \right) \end{align*}$$となる。
これをJuliaで書いてみる。
数式処理にSymPyを使う。
using Plots
using SymPy
using LinearAlgebra
@vars x y z Q;
静電ポテンシャルと電場の式:SymPy
Vsym(x,y,z,Q) = Q/sqrt(x^2+y^2+z^2)
Esym(x,y,z,Q) = -1 .* [diff(Vsym(x,y,z,Q),x), diff(Vsym(x,y,z,Q),y), diff(Vsym(x,y,z,Q),z)]
Esym (generic function with 1 method)
電荷の位置は原点、位置 $(x, y, z)$でのポテンシャルと電場が正しく表示されるか確認。
Vsym(x,y,z,Q)
Esym(x,y,z,Q)
電場の数値計算用の関数を作成
Ecal = lambdify(Esym(x,y,z,Q), [x,y,z,Q])
#122 (generic function with 1 method)
位置 $(-1, 0, 0)$ での電場を計算してみる( $Q = 1$ )。
Ecal(-1,0,0,1)
3-element Vector{Float64}: -1.0 0.0 0.0
x1 = 0; y1 = 0; z1 = 0; # 電荷の座標
Q1 = 1; # 電荷の大きさ
zc = 0; # 計算する位置のz座標
num = 5;
div = 0.3;
xs = -div*(num+0.5) : div : div*(num+0.5) # x座標
ys = -div*(num+0.5) : div : div*(num+0.5) # y座標
xxs = [x for x in xs for y in ys]
yys = [y for x in xs for y in ys]
println("z = ", z1)
EvecNormalize(x, y) = normalize(Ecal(x-x1, y-y1, zc-z1, Q1)) ./ 5
A = [norm(Ecal(x-x1, y-y1, zc-z1, Q1)) for x in xs, y in ys]'
contour(xs, ys, A, color=:Reds, fill=true, line=false, levs=100, aspect_ratio=:equal,
size=(400,300), xlabel="x", ylabel="y",xlims=[xs[begin]-div,xs[end]+div])
quiver!(xxs, yys, quiver=EvecNormalize)
z = 0
x1 = 1; y1 = 0; Q1 = 1; #電荷1の座標と大きさ
x2 = -1; y2 = 0; Q2 = 1; #電荷1の座標と大きさ
zc = 0; # 計算する位置のz座標
num = 10;
div = 0.25;
xs = -div*(num+0.5) : div : div*(num+0.5)
ys = -div*(num+0.5) : div : div*(num+0.5)
xxs = [x for x in xs for y in ys]
yys = [y for x in xs for y in ys]
println("z = ", zc)
EvecNormalize(x, y) = normalize(Ecal(x-x1, y-y1, zc-z1, Q1) + Ecal(x-x2, y-y2, zc-z1, Q2)) ./ 5
A = [norm(Ecal(x-x1, y-y1, zc-z1, 1) + Ecal(x-x2, y-y2, zc-z1, 1)) for x in xs, y in ys]'
contour(xs, ys, A, color=:Reds, fill=true, line=false, levs=100, aspect_ratio=:equal,
size=(400,300), xlabel="x", ylabel="y",xlims=[xs[begin]-div,xs[end]+div])
quiver!(xxs, yys, quiver=EvecNormalize)
z = 0
電場勾配は静電ポテンシャルの位置による2回微分で、以下のように書ける。
$$\begin{align*} V_{\alpha, \beta} = \frac{\partial^2 V}{\partial\alpha\partial\beta} &= \left[ \begin{array}{ccc}\frac{Q \left(2 x^{2} - y^{2} - z^{2}\right)}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{5}{2}}}&\frac{3 Q x y}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{5}{2}}}&\frac{3 Q x z}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{5}{2}}}\\\frac{3 Q x y}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{5}{2}}}&\frac{Q \left(- x^{2} + 2 y^{2} - z^{2}\right)}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{5}{2}}}&\frac{3 Q y z}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{5}{2}}}\\\frac{3 Q x z}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{5}{2}}}&\frac{3 Q y z}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{5}{2}}}&\frac{Q \left(- x^{2} - y^{2} + 2 z^{2}\right)}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{5}{2}}}\end{array}\right]\\ &= \frac{Q}{\left(x^{2} + y^{2} + z^{2}\right)^{\frac{5}{2}}} \left[ \begin{array}{ccc} 2x^{2}-y^{2}-z^{2}&3xy&3xz\\ 3xy&-x^{2}+2y^{2}-z^{2}&3yz\\ 3xz&3yz&-x^{2}-y^{2}+2z^{2} \end{array}\right] \end{align*}$$これをJuliaで書いてみる。
数式処理にSymPyを使う。
電場勾配(テンソル成分)の式:SymPy
V2difsym(a, b) = simplify(diff(Vsym(x,y,z,Q), a, b))
V2difsym (generic function with 1 method)
$\displaystyle{V_{x,x}}$ と $\displaystyle{V_{x,y}}$ を計算して、正しいか確認する。
Vxx = V2difsym(x, x)
Vxy = V2difsym(x, y)
電場勾配テンソル:SymPy
V2difmatsym = [V2difsym(x,x) V2difsym(y, x) V2difsym(z, x)
V2difsym(x,y) V2difsym(y, y) V2difsym(z, y)
V2difsym(x,z) V2difsym(y, z) V2difsym(z, z)
]
@vars a1
simplify(V2difmatsym.subs(x, 0-a1).subs(y, 0).subs(z, 0) + V2difmatsym.subs(x, 0+a1).subs(y, 0).subs(z, 0) +
V2difmatsym.subs(y, 0-a1).subs(z, 0).subs(x, 0) + V2difmatsym.subs(y, 0+a1).subs(z, 0).subs(x, 0) +
V2difmatsym.subs(z, 0-a1).subs(x, 0).subs(y, 0) + V2difmatsym.subs(z, 0+a1).subs(x, 0).subs(y, 0)
)
電場勾配はない。
@vars a2
simplify(V2difmatsym.subs(x, 0-a2).subs(y, 0).subs(z, 0) + V2difmatsym.subs(x, 0+a2).subs(y, 0).subs(z, 0) +
V2difmatsym.subs(y, 0-a2).subs(z, 0).subs(x, 0) + V2difmatsym.subs(y, 0+a2).subs(z, 0).subs(x, 0)
)
$\displaystyle{|V_{z,z}| > |V_{x,x}| = |V_{y,y}|}$ であることから、 $z$軸が電場勾配の最大主軸で、軸対称の異方性を有する(非対称パラメータ $\eta = 0$ )。
@vars a2
simplify(V2difmatsym.subs(x, 0-a2/2).subs(y, 0).subs(z, 0) + V2difmatsym.subs(x, 0+a2/2).subs(y, 0).subs(z, 0) +
V2difmatsym.subs(y, 0-a2).subs(z, 0).subs(x, 0) + V2difmatsym.subs(y, 0+a2).subs(z, 0).subs(x, 0)
)
$\displaystyle{|V_{x,x}| > |V_{z,z}| > |V_{y,y}|}$ であることから、 $x$軸が電場勾配の最大主軸で、非対称パラメータは $\displaystyle{\eta = \frac{1}{5}}$ となる。
@vars a3
simplify(V2difmatsym.subs(x, 0-a3/sympy.sqrt(3)).subs(y, 0).subs(z, 0) +
V2difmatsym.subs(x, 0+a3/(2sympy.sqrt(3))).subs(y, -a3/2).subs(z, 0) +
V2difmatsym.subs(x, 0+a3/(2sympy.sqrt(3))).subs(y, +a3/2).subs(z, 0)
)
$\displaystyle{|V_{z,z}| > |V_{x,x}| = |V_{y,y}|}$ であることから、 $z$軸が電場勾配の最大主軸で、軸対称の異方性を有する(非対称パラメータ $\eta = 0$ )。
数値計算用の関数を作成
V2difmatcal = lambdify(V2difmatsym, [x,y,z,Q])
V2difmatcal2(calpoint, Qpoint, Qval) = V2difmatcal(calpoint[1]-Qpoint[1], calpoint[2]-Qpoint[2], calpoint[3]-Qpoint[3], Qval)
V2difmatcal2 (generic function with 1 method)
calpoint = [0,0,0]
Qpoint1 = [ 1 0 0
0 1 0
-1 0 0
0 -1 0
0 0 1
0 0 -1];
Qval = 1
totalV2 = zeros(3,3)
for iii in eachindex(Qpoint1[:,1])
totalV2 += V2difmatcal2(calpoint, Qpoint1[iii,:], Qval)
end
display(totalV2)
plot([calpoint[1]], [calpoint[2]], [calpoint[3]], marker=:circle, label="calc. point")
plot!(Qpoint1[:,1], Qpoint1[:,2], Qpoint1[:,3],
xlims=[-1,1], ylims=[-1,1], zlims=[-1,1],
xlabel="x", ylabel="y", zlabel="z", label="charge point",
marker=:circle, line=:false, camera=(20,30), aspect_ratio=:equal, size=(350,300))
3×3 Matrix{Float64}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
calpoint = [0,0,0]
Qpoint1 = [ 1 0 0
0 1 0
-1 0 0
0 -1 0];
Qval = 1
totalV2 = zeros(3,3)
for iii in eachindex(Qpoint1[:,1])
totalV2 += V2difmatcal2(calpoint, Qpoint1[iii,:], Qval)
end
display(totalV2)
plot([calpoint[1]], [calpoint[2]], [calpoint[3]], marker=:circle, label="calc. point")
plot!(Qpoint1[:,1], Qpoint1[:,2], Qpoint1[:,3],
xlims=[-1,1], ylims=[-1,1], zlims=[-1,1],
xlabel="x", ylabel="y", zlabel="z", label="charge point",
marker=:circle, line=:false, camera=(20,30), aspect_ratio=:equal, size=(350,300))
3×3 Matrix{Float64}: 2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 -4.0
calpoint = [0,0,0]
Qpoint1 = [ 1/2 0 0
0 1 0
-1/2 0 0
0 -1 0];
Qval = 1
totalV2 = zeros(3,3)
for iii in eachindex(Qpoint1[:,1])
totalV2 += V2difmatcal2(calpoint, Qpoint1[iii,:], Qval)
end
display(totalV2)
plot([calpoint[1]], [calpoint[2]], [calpoint[3]], marker=:circle, label="calc. point")
plot!(Qpoint1[:,1], Qpoint1[:,2], Qpoint1[:,3],
xlims=[-1,1], ylims=[-1,1], zlims=[-1,1],
xlabel="x", ylabel="y", zlabel="z", label="charge point",
marker=:circle, line=:false, camera=(20,30), aspect_ratio=:equal, size=(350,300))
3×3 Matrix{Float64}: 30.0 0.0 0.0 0.0 -12.0 0.0 0.0 0.0 -18.0
calpoint = [0,0,0]
Qpoint1 = [ 2/sqrt(3) 0 0
-1/sqrt(3) 1 0
-1/sqrt(3) -1 0];
Qval = 1
totalV2 = zeros(3,3)
for iii in eachindex(Qpoint1[:,1])
totalV2 += V2difmatcal2(calpoint, Qpoint1[iii,:], Qval)
end
display(totalV2)
plot([calpoint[1]], [calpoint[2]], [calpoint[3]], marker=:circle, label="calc. point")
plot!(Qpoint1[:,1], Qpoint1[:,2], Qpoint1[:,3],
xlims=[-1.5,1.5], ylims=[-1.5,1.5], zlims=[-1,1],
xlabel="x", ylabel="y", zlabel="z", label="charge point",
marker=:circle, line=:false, camera=(20,30), aspect_ratio=:equal, size=(350,300))
3×3 Matrix{Float64}: 0.974279 0.0 0.0 0.0 0.974279 0.0 0.0 0.0 -1.94856