gnuplot で電気力線を数値的に解いて描く

gnuplot だけで，電気力線を数値的に解き，電気力線に沿った電場ベクトルを描く。電場ベクトルが電気力線の接ベクトルであることから，微分方程式をルンゲ・クッタ法で解き，数値解をファイルに保存し，それを読み込んで電場ベクトルの始点の座標としてベクトル場を描く。これを全て gnuplot だけでやる。

参考：

gnuplot だけでプログラミングしたり，数値計算したりする際の参考として

2つの点電荷による電場

位置 $\boldsymbol{r}_1$ にある正の点電荷 $q>0$ と，位置 $\boldsymbol{r}_2$ にある負の点電荷 $\pm q$ がつくる電場ベクトル $\boldsymbol{E}$ は（適宜定数部分の規格化を行って）

\begin{eqnarray}
\boldsymbol{E}
&=& \frac{q}{4\pi \varepsilon_0} \frac{\boldsymbol{r} – \boldsymbol{r}_1}{|\boldsymbol{r} – \boldsymbol{r}_1|^3}
\pm \frac{q}{4\pi \varepsilon_0} \frac{\boldsymbol{r} – \boldsymbol{r}_2}{|\boldsymbol{r} – \boldsymbol{r}_2|^3} \\
&\Rightarrow& \frac{\boldsymbol{r} – \boldsymbol{r}_1}{|\boldsymbol{r} – \boldsymbol{r}_1|^3}
\pm
\frac{\boldsymbol{r} – \boldsymbol{r}_2}{|\boldsymbol{r} – \boldsymbol{r}_2|^3}
\end{eqnarray}

$z=0$ の $xy$ 平面で2次元ベクトル場を描く。

$$\boldsymbol{r} = (x, y, 0), \quad \boldsymbol{r}_1 = (x_1, y_1, 0), \quad \boldsymbol{r}_2 = (x_2, y_2, 0)$$\begin{eqnarray}
E_x(x,y) &=& \frac{x-x_1}{\left((x-x_1)^2 + (y-y_1)^2\right)^{3/2}}
\pm \frac{x-x_2}{\left((x-x_2)^2 + (y-y_2)^2\right)^{3/2}} \\
E_y(x,y) &=& \frac{y-y_1}{\left((x-x_1)^2 + (y-y_1)^2\right)^{3/2}}
\pm \frac{y-y_2}{\left((x-x_2)^2 + (y-y_2)^2\right)^{3/2}}
\end{eqnarray}

規格化した電場 $\hat{\boldsymbol{E}}$ は

\begin{eqnarray}
E(x,y) &=& \sqrt{\left\{E_x(x,y)\right\}^2 + \left\{E_y(x,y)\right\}^2} \\
\hat{E}_x(x,y) &=& \frac{E_x(x,y)}{E(x,y)}\\
\hat{E}_y(x,y) &=& \frac{E_y(x,y)}{E(x,y)}\\
\end{eqnarray}

正負の各1個の点電荷による電気力線

電場ベクトル $\boldsymbol{E}$ が電気力線の接ベクトルであるから，微分方程式

$$\frac{d\boldsymbol{r}}{ds} = \hat{\boldsymbol{E}}$$

を（ルンゲ・クッタ法で）数値的に解いて，電気力線 $\boldsymbol{r}(s)$ を求める。

\begin{eqnarray}
\frac{dx}{ds} &=& \hat{E}_x(x, y) \equiv F_1(x,y) \\
\frac{dy}{ds} &=& \hat{E}_y(x, y) \equiv F_2(x,y)
\end{eqnarray}

$\boldsymbol{r}_1 = (0, 1, 0), \ \boldsymbol{r}_1 = (0, -1, 0)$ の時は，以下の図のようになるから，数値計算の範囲としては，第1象限（$x > 0, \ y > 0$）だけを考えればよいであろう。

denba2

In [1]:
%gnuplot inline svg size 480,480 fixed enhanced font 'Noto Sans CJK JP,14'
In [2]:
# 電場 E
Ex(x, y) = x/sqrt(x**2 + y**2)**3
Ey(x, y) = y/sqrt(x**2 + y**2)**3

# 正電荷と負電荷がつくる電場
E2x(x, y) = Ex(x-x1, y-y1) - Ex(x-x2, y-y2)
E2y(x, y) = Ey(x-x1, y-y1) - Ey(x-x2, y-y2)

# 電場の大きさ
E(x, y) = sqrt(E2x(x, y)**2 + E2y(x, y)**2)

# 規格化された hat E 
F1(x, y) = E2x(x, y)/E(x, y)
F2(x, y) = E2y(x, y)/E(x, y)

# 点電荷の位置
x1 = 0
y1 = 1
x2 = 0
y2 = -1

# 計算範囲
xmin = 0
xmax = 3
ymin = 0
ymax = 3

# ルンゲ・クッタ法のパラメータ
h = 0.05
N = 100
array X[N]
array Y[N]

# データファイルの削除
# あとで append で開くので
! rm -f pmq.txt
! rm -f pmq-line.txt

# 初期条件
# 正の点電荷のまわりの半径 r0 の円上からはじめる
r0 = 0.2
# theta の分割数
Ntheta = 12

do for [j=0:Ntheta]{
  theta = -pi/2 + pi/Ntheta*j
  x0 = x1 + r0 * cos(theta)
  y0 = y1 + r0 * sin(theta)

  # 4次のルンゲ・クッタ法
  do for [i=1: N]{
    X[i] = x0
    Y[i] = y0
    k1 = h*F1(x0, y0)
    m1 = h*F2(x0, y0)
    k2 = h*F1(x0+0.5*k1, y0+0.5*m1)
    m2 = h*F2(x0+0.5*k1, y0+0.5*m1)
    k3 = h*F1(x0+0.5*k2, y0+0.5*m2)
    m3 = h*F2(x0+0.5*k2, y0+0.5*m2)
    k4 = h*F1(x0+k3, y0+m3)
    m4 = h*F2(x0+k3, y0+m3)
    X1 = x0 + (k1 + 2.*k2 + 2.*k3 + k4)/6
    Y1 = y0 + (m1 + 2.*m2 + 2.*m3 + m4)/6
    iend = i
    # 計算範囲を逸脱したら do loop から break
    if (X1 < xmin) {break}
    if (X1 > xmax) {break}
    if (Y1 < ymin) {break}
    if (Y1 > ymax) {break}
    x0 = X1
    y0 = Y1
  }

  # 電気力線の座標データファイル
  # with lines するので，適宜空行を挿入する
  set print "pmq-line.txt" append
  # 第1象限
  do for [i=3:iend]{
    print sprintf("%8.4f %8.4f", X[i], Y[i])
  }
  if (Y[iend] > 0.5) {print " "}
  # 第4象限
  do for [i=iend:3:-1]{
    print sprintf("%8.4f %8.4f", X[i], -Y[i])
  }
  print " "
  # 第2象限
  do for [i=3:iend]{
    print sprintf("%8.4f %8.4f", -X[i], Y[i])
  }
  if (Y[iend] > 0.5) {print " "}
  # 第3象限
  do for [i=iend:3:-1]{
    print sprintf("%8.4f %8.4f", -X[i], -Y[i])
  }
  print " "
  set print

  # 電場ベクトルの始点の座標データファイル
  set print "pmq.txt" append
  # 第1象限
  do for [i=11:iend-5:15]{
    print sprintf("%8.4f %8.4f", X[i], Y[i])
  }
  # 第4象限
  do for [i=11:iend-5:15]{
    print sprintf("%8.4f %8.4f", X[i], -Y[i])
  }
  # 第2象限
  do for [i=11:iend-5:15]{
    print sprintf("%8.4f %8.4f", -X[i], Y[i])
  }
  # 第3象限
  do for [i=11:iend-5:15]{
    print sprintf("%8.4f %8.4f", -X[i], -Y[i])
  }
  set print
}

reset
unset xtics
unset ytics
unset border

set xrange [-3.5:3.5]
set yrange [-3.5:3.5]
set size ratio 1
set zeroaxis

set title "正負の点電荷による電場"
unset colorbox
set palette cubehelix start 0 cycles 0 saturation 3 negative

# 正電荷の位置に青丸と +
set label 1 point pt 6 ps 2 lc "black" at x1, y1
set label 2 center at first x1, y1 "＋" tc "black" font ",24"

# 負電荷の位置に赤丸と -
set label 3 point pt 6 ps 2 lc "black" at x2, y2
set label 4 center at first x2, y2 "ー" tc "black" font ",24"

# 矢印の長さ調整用
scaling = 0.05
set key sample 1

plot "pmq-line.txt" u 1:2:(log(E($1,$2))) w l lw 3 lc palette notitle, \
  "pmq.txt"  \
  u 1:2:(scaling*F1($1,$2)):(scaling*F2($1,$2)):(log(E($1,$2))) \
  w vec lc palette lw 6 filled head title "電場"

In [3]:
set output "./vec-pmq.svg"
replot
set output

2つの正電荷による電気力線

In [4]:
# 電場 E
Ex(x, y) = x/sqrt(x**2 + y**2)**3
Ey(x, y) = y/sqrt(x**2 + y**2)**3

# 正電荷と負電荷がつくる電場
E2x(x, y) = Ex(x-x1, y-y1) + Ex(x-x2, y-y2)
E2y(x, y) = Ey(x-x1, y-y1) + Ey(x-x2, y-y2)

# 電場の大きさ
E(x, y) = sqrt(E2x(x, y)**2 + E2y(x, y)**2)

# 規格化された hat E 
F1(x, y) = E2x(x, y)/E(x, y)
F2(x, y) = E2y(x, y)/E(x, y)

# 点電荷の位置
x1 = 0
y1 = 1
x2 = 0
y2 = -1

# 計算範囲
xmin = 0
xmax = 3
ymin = 0
ymax = 3

# ルンゲ・クッタ法のパラメータ
h = 0.05
N = 100
array X[N]
array Y[N]

# データファイルの削除
# あとで append で開くので
! rm -f ppq.txt
! rm -f ppq-line.txt

# 初期条件
# 正の点電荷のまわりの半径 r0 の円上からはじめる
r0 = 0.2
# theta の分割数
Ntheta = 12

do for [j=2:Ntheta]{
  theta = -pi/2 + pi/Ntheta*j
  x0 = x1 + r0 * cos(theta)
  y0 = y1 + r0 * sin(theta)

  # 4次のルンゲ・クッタ法
  do for [i=1: N]{
    X[i] = x0
    Y[i] = y0
    k1 = h*F1(x0, y0)
    m1 = h*F2(x0, y0)
    k2 = h*F1(x0+0.5*k1, y0+0.5*m1)
    m2 = h*F2(x0+0.5*k1, y0+0.5*m1)
    k3 = h*F1(x0+0.5*k2, y0+0.5*m2)
    m3 = h*F2(x0+0.5*k2, y0+0.5*m2)
    k4 = h*F1(x0+k3, y0+m3)
    m4 = h*F2(x0+k3, y0+m3)
    X1 = x0 + (k1 + 2.*k2 + 2.*k3 + k4)/6
    Y1 = y0 + (m1 + 2.*m2 + 2.*m3 + m4)/6
    iend = i
    # 計算範囲を逸脱したら do loop から break
    if (X1 < xmin) {break}
    if (X1 > xmax) {break}
    if (Y1 < ymin) {break}
    if (Y1 > ymax) {break}
    x0 = X1
    y0 = Y1
  }

  # 電気力線の座標データファイル
  # with lines するので，適宜空行を挿入する
  set print "ppq-line.txt" append
  # 第1象限
  do for [i=3:iend]{
    print sprintf("%8.4f %8.4f", X[i], Y[i])
  }
  print " "
  # 第4象限
  do for [i=iend:3:-1]{
    print sprintf("%8.4f %8.4f", X[i], -Y[i])
  }
  print " "
  # 第2象限
  do for [i=3:iend]{
    print sprintf("%8.4f %8.4f", -X[i], Y[i])
  }
  print " "
  # 第3象限
  do for [i=iend:3:-1]{
    print sprintf("%8.4f %8.4f", -X[i], -Y[i])
  }
  print " "
  set print

  # 電場ベクトルの始点の座標データファイル
  set print "ppq.txt" append
  # 第1象限
  do for [i=11:iend-5:15]{
    print sprintf("%8.4f %8.4f", X[i], Y[i])
  }
  # 第4象限
  do for [i=11:iend-5:15]{
    print sprintf("%8.4f %8.4f", X[i], -Y[i])
  }
  # 第2象限
  do for [i=11:iend-5:15]{
    print sprintf("%8.4f %8.4f", -X[i], Y[i])
  }
  # 第3象限
  do for [i=11:iend-5:15]{
    print sprintf("%8.4f %8.4f", -X[i], -Y[i])
  }
  set print
}

reset
unset xtics
unset ytics
unset border

set xrange [-3.5:3.5]
set yrange [-3.5:3.5]
set size ratio 1
set zeroaxis

set title "2つの正電荷による電場"
unset colorbox
set palette cubehelix start 0 cycles 0 saturation 3 negative

# 正電荷の位置に青丸と +
set label 1 point pt 6 ps 2 lc "black" at x1, y1
set label 2 center at first x1, y1 "＋" tc "black" font ",24"
set label 3 point pt 6 ps 2 lc "black" at x2, y2
set label 4 center at first x2, y2 "＋" tc "black" font ",24"

# 矢印の長さ調整用
scaling = 0.1
set key sample 1

plot "ppq-line.txt" u 1:2:(log(E($1,$2))) w l lw 3 lc palette notitle, \
  "ppq.txt"  \
  u 1:2:(scaling*F1($1,$2)):(scaling*F2($1,$2)):(log(E($1,$2))) \
  w vec lc palette lw 6 filled head title "電場"

In [5]:
set output "./vec-ppq.svg"
replot
set output

電気双極子による電場

原点にある電気双極子 $\boldsymbol{p}$ がつくる電場ベクトル $\boldsymbol{E}$ は（適宜定数部分の規格化を行って）

$$\boldsymbol{E} \Rightarrow 3 \frac{\boldsymbol{r}\cdot\boldsymbol{p}}{r^5} \boldsymbol{r}
– \frac{\boldsymbol{p}}{r^3}$$

$\boldsymbol{p} = (0, 1, 0), \ z = 0$ として $xy$ 平面で2次元ベクトル場を描く。

\begin{eqnarray}
E_x(x,y) &=& \frac{3}{\left( x^2 + y^2 \right)^{5/2}} x y \\
E_y(x,y) &=& \frac{3}{\left( x^2 + y^2 \right)^{5/2}} y^2 – \frac{1}{\left( x^2 + y^2 \right)^{3/2}} \\
\end{eqnarray}

規格化した電場 $\hat{\boldsymbol{E}}$ は

\begin{eqnarray}
E(x,y) &=& \sqrt{\left\{E_x(x,y)\right\}^2 + \left\{E_y(x,y)\right\}^2} \\
\hat{E}_x(x,y) &=& \frac{E_x(x,y)}{E(x,y)} \equiv F_1(x,y)\\
\hat{E}_y(x,y) &=& \frac{E_y(x,y)}{E(x,y)} \equiv F_2(x,y)\\
\end{eqnarray}

電気双極子による電気力線

電場ベクトル $\boldsymbol{E}$ が電気力線の接ベクトルであるから，微分方程式

$$\frac{d\boldsymbol{r}}{ds} = \hat{\boldsymbol{E}}$$

を（ルンゲ・クッタ法で）数値的に解いて，電気力線 $\boldsymbol{r}(s)$ を求める。

\begin{eqnarray}
\frac{dx}{ds} &=& \hat{E}_x(x, y) \equiv F_1(x,y) \\
\frac{dy}{ds} &=& \hat{E}_y(x, y) \equiv F_2(x,y)
\end{eqnarray}

電気双極子が原点にある場合は，以下の図のようになるから，数値計算の範囲としては，第1象限（$x > 0, \ y > 0$）だけを考えればよいであろう。

dipole

In [6]:
# 電場 E
Ex(x, y) = 3.* x*y/sqrt(x**2 + y**2)**5 
Ey(x, y) = 3.* y*y/sqrt(x**2 + y**2)**5 - 1./sqrt(x**2 + y**2)**3

# 規格化された hat E
E(x, y) = sqrt(Ex(x, y)**2 + Ey(x, y)**2)
F1(x,y) = Ex(x, y)/E(x, y)
F2(x,y) = Ey(x, y)/E(x, y)

# 計算範囲
xmin = 0
xmax = 3
ymin = 0
ymax = 3

# ルンゲ・クッタ法のパラメータ
h = 0.05
N = 100
array X[N]
array Y[N]

# データファイルの削除
# あとで append で開くので
! rm -f dipole.txt
! rm -f dipole-line.txt

# 初期条件
# 原点のまわりの半径 r0 の円上からはじめる
r0 = 0.22
# theta の分割数
Ntheta = 8
th0 = pi/3
thend = pi/2
do for [j=1:Ntheta]{
  theta = th0 + (thend-th0)/Ntheta*j
  x0 = r0 * cos(theta)
  y0 = r0 * sin(theta)

  # 4次のルンゲ・クッタ法
  do for [i=1: N]{
    X[i] = x0
    Y[i] = y0
    k1 = h*F1(x0, y0)
    m1 = h*F2(x0, y0)
    k2 = h*F1(x0+0.5*k1, y0+0.5*m1)
    m2 = h*F2(x0+0.5*k1, y0+0.5*m1)
    k3 = h*F1(x0+0.5*k2, y0+0.5*m2)
    m3 = h*F2(x0+0.5*k2, y0+0.5*m2)
    k4 = h*F1(x0+k3, y0+m3)
    m4 = h*F2(x0+k3, y0+m3)
    X1 = x0 + (k1 + 2.*k2 + 2.*k3 + k4)/6
    Y1 = y0 + (m1 + 2.*m2 + 2.*m3 + m4)/6
    iend = i
    # 計算範囲を逸脱したら do loop から break
    if (X1 < xmin) {break}
    if (X1 > xmax) {break}
    if (Y1 < ymin) {break}
    if (Y1 > ymax) {break}
    x0 = X1
    y0 = Y1
  }

  # 電気力線の座標データファイル
  # with lines するので，適宜空行を挿入する
  set print "dipole-line.txt" append
  # 第1象限
  do for [i=1:iend]{
    print sprintf("%8.4f %8.4f", X[i], Y[i])
  }
  if (Y[iend] > 0.5) {print " "}
  # 第4象限
  do for [i=iend:1:-1]{
    print sprintf("%8.4f %8.4f", X[i], -Y[i])
  }
  print " "
  # 第2象限
  do for [i=1:iend]{
    print sprintf("%8.4f %8.4f", -X[i], Y[i])
  }
  if (Y[iend] > 0.5) {print " "}
  # 第3象限
  do for [i=iend:1:-1]{
    print sprintf("%8.4f %8.4f", -X[i], -Y[i])
  }
  print " "
  set print  

  # 電場ベクトルの始点の座標データファイル
  set print "dipole.txt" append
  # 第1象限
  do for [i=25:iend-5:15]{
    print sprintf("%8.4f %8.4f", X[i], Y[i])
  }
  # 第4象限
  do for [i=25:iend-5:15]{
    print sprintf("%8.4f %8.4f", X[i], -Y[i])
  }
  # 第2象限
  do for [i=25:iend-5:15]{
    print sprintf("%8.4f %8.4f", -X[i], Y[i])
  }
  # 第3象限
  do for [i=25:iend-5:15]{
    print sprintf("%8.4f %8.4f", -X[i], -Y[i])
  }
  set print  
}

reset
unset xtics
unset ytics
unset border

set xrange [-3:3]
set yrange [-3:3]
set size ratio 1
set zeroaxis

set title "電気双極子による電場"
unset colorbox
set palette cubehelix start 0 cycles 0 saturation 3 negative

# 電気双極子 p
set arrow 1 from 0, -0.2 to 0, 0.05 filled head lw 6 lc "black" lt 1

# 矢印の長さ調整用
scaling = 0.1
set key sample 1
plot "dipole-line.txt" u 1:2:(log(E($1,$2))) w l lw 3 lc palette notitle, \
  "dipole.txt"   \
  u 1:2:(scaling*F1($1,$2)):(scaling*F2($1,$2)):(log(E($1,$2))) \
  w vec lc palette lw 6 filled head title "電場"

In [7]:
set output "./vec-ele-dipole.svg"
replot
set output

磁気双極子による磁力線

電気双極子による電気力線と同じ。色だけ変えて描いてみる。

In [8]:
set title "磁気双極子による磁場"
unset colorbox
set palette cubehelix start 2 cycles 0 saturation 3 negative

# 電気双極子 p
set arrow 1 from 0, -0.2 to 0, 0.05 filled head lw 6 lc "red" lt 1

# 矢印の長さ調整用
scaling = 0.1
set key sample 1

plot "dipole-line.txt" u 1:2:(log(E($1,$2))) w l lw 3 lc palette notitle, \
  "dipole.txt"   \
  u 1:2:(scaling*F1($1,$2)):(scaling*F2($1,$2)):(log(E($1,$2))) \
  w vec lc palette lw 6 filled head title "磁場"

In [9]:
set output "./vec-mag-dipole.svg"
replot
set output