Maxima で電場ベクトルの方向場を描く

正負の点電荷がつくる電場の向きを表す方向場を Maxima で描く。以下のページの抜粋。

正負の点電荷がつくる電場

$\boldsymbol{r}_1$ にある正電荷 $q\ (>0)$ と，$\boldsymbol{r}_2$ にある負電荷 $-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}
– \frac{q}{4\pi \varepsilon_0} \frac{\boldsymbol{r} – \boldsymbol{r}_2}{|\boldsymbol{r} – \boldsymbol{r}_2|^3}
\end{eqnarray}

$z=0$ の $xy$ 平面で2次元ベクトル場を描く。適宜無次元化を行い，
$\displaystyle \frac{q}{4\pi \varepsilon_0}\rightarrow 1$ とする。また，

$$\boldsymbol{r} = (x, y, 0), \quad \boldsymbol{r}_1 = (0, 1, 0),
\quad \boldsymbol{r}_2 = (0, -1, 0)$$

とすると，

$$E_x = \frac{x}{\left(x^2 + (y-1)^2\right)^{3/2}}-\frac{x}{\left(x^2 + (y+1)^2\right)^{3/2}}$$$$E_y = \frac{y-1}{\left(x^2 + (y-1)^2\right)^{3/2}}-
\frac{y+1}{\left(x^2 + (y+1)^2\right)^{3/2}}$$

$\boldsymbol{E}$ の向きを表す方向場（単位ベクトル場）

一般にベクトルは大きさと向きを表す。電場 $\boldsymbol{E}$ の大きさは電荷からの距離の2乗に反比例するので，]ベクトルの長さが大きく変化する。見やすくするために，以下のように $\boldsymbol{E}$ の向きを表す方向場，すなわち大きさが $1$ の単位ベクトル場を定義して，ベクトルの向きのみを正しく描画する。

$$\hat{\boldsymbol{E}} \equiv \frac{\boldsymbol{E}}{\sqrt{\boldsymbol{E}\cdot\boldsymbol{E}}}$$$$\hat{E}_x = \frac{E_x}{\sqrt{E_x^2 + E_y^2}}$$$$\hat{E}_y = \frac{E_y}{\sqrt{E_x^2 + E_y^2}}$$

draw2d() で電場の方向場を描く

グラフのサイズ等のデフォルト設定

draw2d() でその都度オプション設定することもできるが，set_draw_defaults() で draw2d() のデフォルト設定ができる。ただし，この set_draw_defaults() がちょっと曲者で，例え追加で dimensions を設定したくても，追加のオプションだけでなく，その時点で設定されている項目を全て列挙しなければならない。その割には，現在の draw のデフォルト設定がどうなっているかを知る関数はない。引数なしで set_draw_defaults(); を実行すると，全ての設定事項が削除されてしまう。

In [1]:
set_draw_defaults(
  /* 弘大 JupyterHub ではインライン SVG 表示 */
  file_name="~/.maxplot", terminal='svg,
  /* グラフのサイズ */
  dimensions=[640,640],
  /* フォント設定 */
  font="Noto Sans CJK JP", font_size=16,
  /* ベクトルの矢印のデフォルト設定 */
  head_length=0.1, head_angle=20)$

電場ベクトルの定義

電場ベクトルを x, y の2変数関数として定義。

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, y-1) - Ex(x, y+1)$
E2y(x, y):= Ey(x, y-1) - Ey(x, y+1)$

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

/* 規格化された電場 hat E */
hE2x(x, y):= scaling * E2x(x, y)/E(x, y)$
hE2y(x, y):= scaling * E2y(x, y)/E(x, y)$

ベクトル場を描く

Maxima の draw2d() でベクトルを描く基本は ここ を参照。

点電荷が $y$ 軸上にあるので，リストを作成する際，ベクトルの始点が $y$ 軸に来ないようにしておく。（点電荷上では電場の分母がゼロになってしまうから。）
以下の例では，coordxy の $x$ 座標は

$$-5.5, -4.5, \cdots, -0.5, 0.5, \cdots, 5.5$$

となり，$y$ 軸上 ($x = 0$) に始点がこない。

In [3]:
/* ベクトルの始点のリストを作成 */
coordxy: makelist()$
for i: 0 thru 11 do(
  for j: 0 thru 11 do(
    coordxy: append(coordxy, [[-5.5+i, -5.5+j]])
  )
)$

/* ベクトルは vector([始点のx, 始点のy], [x成分, y成分]) */
vechE(x, y):= vector([x, y], [hE2x(x, y), hE2y(x, y)])$
In [4]:
/* ベクトルの長さ調整用 */
scaling: 0.5$

draw2d(
  /* 表示範囲 */
  xrange = [-6.5, 6.5], yrange = [-6.5, 6.5], 
  xtics = 1, ytics = 1, grid=true,

  line_width = 2, 
  /* ベクトル場の表示 */
  makelist(vechE(k[1], k[2]), k, coordxy)
)$

点電荷の位置にラベル

正電荷，負電荷の位置にラベルをつけてみる。

In [5]:
draw2d(
  /* 表示範囲 */
  xrange = [-6.5, 6.5], yrange = [-6.5, 6.5], 
  xtics = 1, ytics = 1, grid=true,

  line_width = 2, 
  /* ベクトル場の表示 */
  makelist(vechE(k[1], k[2]), k, coordxy), 
  
  /* 正電荷の位置に青丸と + */
  color = blue, point_type = 6, point_size = 2,   
  points([[0, 1]]),
  label(["＋", 0, 1]), 
  
  /* 負電荷の位置に赤丸と - */
  /* 一旦白く塗りつぶす*/
  color = white, point_type = 7, point_size = 2,   
  points([[0, -1]]),
  /* あらためて赤丸を描く */
  color = red, point_type = 6, point_size = 2,   
  points([[0, -1]]),
  label(["ー", 0, -1])

)$

一応の完成作

Maxima ではベクトルの大きさに応じて色を palette にすることができないようなので，これで一応の完成作。サンプル数を増やし，少し短めのベクトルをたくさん描いてみた。

In [6]:
/* ベクトルの始点のリストを作成 */
/* -5 <= x <= 5, -5 <= y <= 5 で 0.5 刻みの格子点を始点に */
coordxy: makelist()$
for i: 0 thru 21 do(
  for j: 0 thru 21 do(
    coordxy: append(coordxy, [[-5.25+i/2, -5.25+j/2]])
  )
)$

/* 2次元ベクトルは vector([始点のx, 始点のy], [x成分, y成分]) */
vechE(x, y):= vector([x, y], [hE2x(x, y), hE2y(x, y)])$

/* ベクトルの長さ調整用 */
scaling: 0.3$

draw2d(
  /* 表示範囲 */
  xrange = [-6, 6], yrange = [-6, 6], 
  /* 目盛と外枠の非表示 */
  xtics = false, ytics = false,
  user_preamble = "unset border;",

  line_width = 4, 
  /* ベクトル場の表示 */
  makelist(vechE(k[1], k[2]), k, coordxy), 
  
  /* 正電荷の位置に青丸と + */
  color = blue, point_type = 6, point_size = 2,   
  points([[0, 1]]),
  label(["＋", 0, 1]), 
  
  /* 負電荷の位置に青丸と - */
  /* 一旦白く塗りつぶす*/
  color = white, point_type = 7, point_size = 2,   
  points([[0, -1]]),
  /* あらためて赤い輪を描く */
  color = red, point_type = 6, point_size = 2,   
  points([[0, -1]]),
  label(["ー", 0, -1])
)$

drawdf() で電場の方向場を描く

2次元ベクトルの方向場（大きさが一定で向きだけを表すベクトル場）に限れば，drawdf パッケージを使ってもっと簡単に描ける。

In [7]:
/* 最初に1回だけロードする */
load(drawdf)$
In [8]:
/* 電場 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, y-1) - Ex(x, y+1)$
E2y(x, y):= Ey(x, y-1) - Ey(x, y+1)$

方向場を描く

drawdf() はベクトル場の $x$ 成分，$y$ 成分と表示範囲を指定するだけで，ベクトルの長さを自動的に調節して描いてくれる。

In [9]:
drawdf(
  /* ベクトル場の x 成分, y 成分 */
  [E2x(x, y), E2y(x, y)], 
  /* 表示範囲 */
  [x, -5, 5], [y, -5, 5])$

グリッド等のオプションの設定

field_grid を指定して，grid = true でみると，drawdf() はグリッド線で囲まれた領域の中に自動的に長さを調節してベクトルを描いていることがわかる。

方向ベクトル場の色は field_color で指定するが，太さについては設定方法が見当たらない。
（line_width は効かない。）

In [10]:
drawdf(
  [E2x(x, y), E2y(x, y)], /* ベクトル場の x 成分, y 成分 */ 
  [x, -5, 5], [y, -5, 5], /* 表示範囲 */ 
  /* グリッドの行と列の数 */
  field_grid = [10, 10], 
  /* 目盛の設定とグリッドの表示 */
  xtics = 1, ytics = 1, grid = true, 
  /* ベクトル場の色 */
  field_color = "blue"
)$

点電荷の位置にラベル

正電荷，負電荷の位置にラベルをつけてみる。

In [11]:
drawdf(
  [E2x(x, y), E2y(x, y)], 
  [x, -5, 5], [y, -5, 5], 
  field_grid = [10, 10], 
  xtics = 1, ytics = 1, grid = true, 
  field_color = "blue", 

  /* 正電荷の位置を一旦白く塗りつぶし... */
  color = white, point_type = 7, point_size = 2,   
  points([[0, 1]]),
  /* あらためて青丸を描く */
  color = blue, point_type = 6, point_size = 2,   
  points([[0, 1]]),
  label(["＋", 0, 1]), 
  
  /* 負電荷の位置を一旦白く塗りつぶし... */
  color = white, point_type = 7, point_size = 2,   
  points([[0, -1]]),
  /* あらためて赤円を描く */
  color = red, point_type = 6, point_size = 2,   
  points([[0, -1]]),
  label(["ー", 0, -1])
)$

field_degree オプション

field_degree='solns と設定すると，（しばらく計算したあとに）方向場が 4 次の Runge Kutta で計算された小さな解曲線で示されるのがおもしろい。（弘大 JupyterHub では延々とスクロールしてはじめてグラフを確認できる。）電場の場合は，グリッド内の電気力線の一部が表示される。矢印の頭が円になるのもご愛嬌。（変更する方法が見当たらない。）

field_degree=2 の場合は2次スプラインで示す（ので，必ずしも正確な力線になっているとは限らない）。

In [12]:
drawdf(
  [E2x(x, y), E2y(x, y)], 
  [x, -5, 5], [y, -5, 5], 
  field_grid = [10, 10], 
  xtics = 1, ytics = 1, grid = true, 
  field_color = "blue", 

  color = white, point_type = 7, point_size = 2,   
  points([[0, 1]]),
  color = blue, point_type = 6, point_size = 2,   
  points([[0, 1]]),
  label(["＋", 0, 1]), 
  
  color = white, point_type = 7, point_size = 2,   
  points([[0, -1]]),
  color = red, point_type = 6, point_size = 2,   
  points([[0, -1]]),
  label(["ー", 0, -1]),
  
  field_degree = 'solns /* 1, 2, or 'solns */
)$

一応の完成作

In [14]:
drawdf(
  [E2x(x, y), E2y(x, y)], 
  [x, -5.5, 5.5], [y, -5.5, 5.5], 
  field_grid = [22, 22], 
  xtics = false, ytics = false, user_preamble="unset border;", 
  field_color = "blue", 

  color = white, point_type = 7, point_size = 2,   
  points([[0, 1]]),
  color = blue, point_type = 6, point_size = 2,   
  points([[0, 1]]),
  label(["＋", 0, 1]), 
  
  color = white, point_type = 7, point_size = 2,   
  points([[0, -1]]),
  color = red, point_type = 6, point_size = 2,   
  points([[0, -1]]),
  label(["ー", 0, -1]),
  
  field_degree = 'solns /* 1, 2, or 'solns */
)$