きっかけ

電場(電界)の数値解析について研究しているが、授業など真面目に聞かず遊んでばかりいたため、基礎知識が不足していると感じた。
電場計算の基礎的な知識を得るために、『数値電界計算の基礎と応用』という本を買った。
普通、電磁気系の数値解析の本といえば、電場計算、磁場計算、電磁波計算などが合わせて一つの本となっているのだが、この本は電場計算のみに特化した素晴らしい本である。

ただ、買っても読まないと意味ないので今回自分用のアウトプットのためにも、学んだことをここに書き留めたいと思う。
また、この本には数式の記述や、理論的な話はあっても、サンプルプログラムなどはない。
ただ書いていることを自分なりに解釈して記事にするだけでは面白くないので、pythonで自分で簡単にサンプルプログラムを実装していく。

なお、上記の本は電界と称されているが、本記事では電場に統一する。

電磁気学の理論・公式

当然だが電磁気学の理論をある程度知っていなければ、電界計算をできるわけない。
といっても、今回主眼に置くのは数値解析の方のため、超簡単に、必要なものだけ話していく。

今更だが、数学的な知識レベルとしては大学一年生～二年生レベルで話していこうと思う。
といっても、$\nabla$の記号やテイラー展開、あとは行列をなんとなく知っていれば多分大丈夫だ。

物理的な話なんてどうでもいい！って人はポアソン方程式まで飛ばして欲しい。

クーロンの法則

電磁気学の教科書で最初に出てくるのはたいていクーロンの法則だ。
高校の教科書では次のように書かれていると思う。

$$
F=k \frac{q_1 q_2}{r^2} \tag{1}
$$
ご存知の通り、$F$は力、$q1$, $q2$は２つの粒子の電荷、そして$r$は距離である。

(高校でも習うと思うが)大学になると、この比例定数$k$は真空の誘電率$\varepsilon_0$を利用して$k=1/(4\pi \varepsilon_0)$と表現するのが普通になる。
さらに、ベクトル演算に慣れてくると、力$F$に方向の意味を持たせるためにも、次のように両辺ともベクトルで書くのが簡潔だ。

念の為行っておくが、ボールド体はベクトルを意味している。

\begin{align}
\boldsymbol{F} 
&= \frac{1}{ 4\pi \varepsilon_0 } \frac{q_1 q_2}{|\boldsymbol{r}|^2} \boldsymbol{\hat{r}} \\
&= \frac{1}{ 4\pi \varepsilon_0 } \frac{q_1 q_2}{|\boldsymbol{r}|^3} \boldsymbol{r} \tag{2}
\end{align}

$\boldsymbol{r}$は２つの粒子を結ぶベクトルで、$\boldsymbol{\hat{r}}$は$\boldsymbol{r}$と同じ向きで、長さ$1$のベクトルである。

一つ目は式(1)にただ２つの粒子を結ぶベクトルを加えただけで、最初はこっちのほうが見やすいと思う。
二つ目は$\boldsymbol{\hat{r}}$を、$\boldsymbol{\hat{r}}=\boldsymbol{r}/|\boldsymbol{r}|$と置き直しただけである。
個人的には、こっちのほうが文字の種類が少なくなって好きだ。プログラム的にもこっちのほうが実装しやすそう。
今後もこっちの書き方で書いていくと思う。

上記の話は現実世界と同じ三次元の話であり、二次元だと少し違う。
どうしてこうなるか、という話はせず、公式だけ上げておこう。

$$
\boldsymbol{F}
= \frac{1}{ 2\pi \varepsilon_0 } \frac{q_1 q_2}{|\boldsymbol{r}|^2} \boldsymbol{r}
$$

重要なのは、$2\pi$になっているところ、そして分母の乗数が一つ減っているところだ。
球の表面積が$4\pi$なのに対して円の周の長さが$2\pi$であること、三次元から二次元に次元が一つ減ったことからなんとなくわかってほしい。

電場と電束密度

クーロンの法則は２つの電荷間に働く力の関係を表す式である。
このことから、電荷を１つ置くことで、この世界に何らかの場ができて、もう１つの電荷を置いたときにその場によって力が働くと考えることができる。
……何言っているかわからないかもしれないので、式をを利用して説明する。

クーロンの法則を$\boldsymbol{F} = q_2 \boldsymbol{E}$のように書き直したとき、
$$
\boldsymbol{E}= \frac{1}{4\pi \varepsilon_0} \frac{\boldsymbol{r}}{|\boldsymbol{r}|^3}q_1
$$
が電荷$q_1$が作る電場$\boldsymbol{E}$である。

二次元のときは、クーロンの法則と同様に次のように表わされることも一応書いておこう。
$$
\boldsymbol{E}
= \frac{1}{ 2\pi \varepsilon_0 } \frac{\boldsymbol{r}}{|\boldsymbol{r}|^2} q_1
$$

ちょっと話が飛躍するが、微小領域での電荷$q$の値、電荷密度$\rho$を利用してうまいこと変形すると次式が得られる。
(この式の導出、当初はもう少し詳しく書こうと思ったのだが、いろいろな前提知識を説明しなければならないため省く。気になる人は自分で調べてみてね)
$$
\nabla \cdot \boldsymbol{E} = \frac{\rho}{\varepsilon_0}
$$

なお、二次元のときも上式は変わらない。

ところで、電束密度$\boldsymbol{D}$というものがある。
この電束密度$\boldsymbol{D}$と電場$\boldsymbol{E}$は、真空中では真空の誘電率$\varepsilon_0$を利用して、次のような比例関係が得られる。
$$
\boldsymbol{D}=\varepsilon_0 \boldsymbol{E}
$$
とても簡単な式だが、この$\boldsymbol{D}$と$\boldsymbol{E}$の定義の物理的な意味は全く違うことには注意しよう。
電束密度は電荷が作る電気力線の密度、電場は電荷に力を及ぼす場である。

ポアソン方程式

電位$\phi$は、次のように定義される。
$$-\nabla \phi = \boldsymbol{E} \tag{3}$$
この式は、電場$\boldsymbol{E}$という3次元ベクトルを1つのスカラー値に直すとても便利な式である。
このような値をスカラーポテンシャルと言う。

$\nabla \cdot \boldsymbol{E} = \rho/\varepsilon_0$に上式を代入すると、次式を求めることができる。

$$\nabla \cdot \nabla \phi = -\frac{\rho}{\varepsilon_0}$$
この式を解けば、電位を求められるし、それを微分することで電場も求めることができる。

右辺の$\varepsilon_0$などは定数なので、わざわざこれから書くのも面倒なので、もう少し一般的に書こう。
左辺を$u$に、右辺を$f$に置き直すと、次のように書き直すことができる。
$$\nabla \cdot \nabla u = f$$
このような式をポアソン方程式と呼ぶ。右辺のマイナスはとらないこともあるが、今回は書くのが面倒なので、とることにした。

長々と書いてきたが、結局、電場を求めるには、このポアソン方程式を解いて、電位を求めてから式(3)から電場を求めればOKだ。

境界条件(Boundary Condition, B.C.)

最後に、境界条件についてちょっとだけ話そう

無限遠点を電位の基準点とする(無限遠点での電位を0とする)と言う話はよく聞くと思う。
これは、境界条件を「無限遠点で0とする」として、上記のポアソン方程式を解く、と言う意味である。

数値解析では、ある計算領域内の値を求めるため、無限遠点を0と置くことは少ない(数値解析は無限遠点を境界条件とする問題は苦手であることが多い)。
そのため、境界条件の電位を基準値として解くことになる。

境界条件は主に以下の２つがある。

ディリクレ境界条件(第一境界条件)：計算領域の境界上において、$u=g$である。
ノイマン境界条件(第二境界条件)：計算領域の境界上において、$\displaystyle \frac{\partial u}{\partial n}=g$($n$は境界の法線方向の微分)

上の$g$は何らかの関数である。多分、定数のことが多いと思う。
他にも、周期境界条件やロビン条件(第三境界条件)などがあるが、今回説明は省く。

ポアソン方程式は、計算領域がディリクレあるいはノイマン境界条件で囲まれていれば解が存在することが証明されている。
注意すべき点として、ノイマン境界条件のみで囲まれている場合は基準値が指定されていないため、どこか任意の場所を基準値として指定する必要がある。

有限差分法(FDM)

ここまで来てついに数値解析の話に入っていく。長かった……。

上で述べたように、解くべき問題はポアソン方程式であり、計算領域はディリクレ／ノイマン境界条件に囲まれている。
これを、数式で表してまとめてみよう。

計算領域を$\Omega$、その境界を$\Gamma$としたとき、

\begin{align}
\nabla \cdot \nabla u = f \ \ &(\text{ in } \Omega )\\
u = g_1(x,y) \ \ &(\text{ on } \Gamma_1) \\
\frac{\partial u}{\partial n} = g_2(x,y) \ \ &(\text{ on } \Gamma_2)
\end{align}

ただし、$\Gamma_1 \cup \Gamma_2 = \Gamma$, $\Gamma_1 \cap \Gamma_2 = \varnothing$である。

ちょっと抽象的な書き方なので、具体的な例で話す。
計算領域$\Omega$を長方形とする。そして、境界$\Gamma_1$は上下の辺、$\Gamma_2$は左右の辺とすると、次の図のようになる。

図１：計算領域

ところで、$\Gamma_1 \cup \Gamma_2 = \Gamma$, $\Gamma_1 \cap \Gamma_2 = \varnothing$の部分は流し読みせずちょっと立ち止まって考えて見てほしい。そんなに難しくないはずだ。
$\Gamma_1$と$\Gamma_2$の領域は被らず、その和集合が$\Gamma$となるということも認識しておこう。わからない人はこちら「数式を読めるエンジニアになるために」(宣伝)。

このように、境界条件を与えられた偏微分問題を解くことを、境界値問題という。

この境界値問題を数値解析するに当たって、大抵最も最初に取り扱われあるのがこの有限差分法(Finite Difference Method: FDM) である。

なお、これ以降、二次元ポアソン方程式について考えていくことにする。
というのも、三次元だと図示するのが面倒だったりするからである。まあ、二次元さえわかれば三次元への拡張は容易だ。
即ち、次の偏微分方程式を解くことを考える。
$$
\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=f
$$

テイラー展開

まず、テイラー展開について話しておこう。
FDMはこれが全てであり、これさえわかれば極めたも同然だ。

一次元のテイラー展開

$\Delta_x$を極めて小さい値として、一次元のテイラー展開は次のように表される。

\begin{align}
u(x+\Delta_x) &= u(x) + \Delta_x \frac{d u(x)}{d x} + \frac{\Delta_x^2}{2!} \frac{d^2 u(x)}{d x^2} +
\frac{\Delta_x^3}{3!} \frac{d^3 u(x)}{d x^3} + \cdots \\
&= \sum_{p=0}^\infty \frac{\Delta_x^p}{p!} \frac{d^p u(x)}{d x^p}
\end{align}

本当は、収束条件などがあるのが、今回は面倒なので無視しよう。

$\Delta_x$が小さい値なので、$\Delta_x^3$以降の項を無視してもまあ大体大丈夫だと考えよう。
$$
u(x+\Delta_x) \approx u(x) + \Delta_x \frac{d u(x)}{d x} + \frac{\Delta_x^2}{2!} \frac{d^2 u(x)}{d x^2}
$$
これは、$\Delta_x^3$程度の誤差があると言うことで、次のように表記する。
$$
u(x+\Delta_x)
= u(x) + \Delta_x \frac{d u(x)}{d x} + \frac{\Delta_x^2}{2!} \frac{d^2 u(x)}{d x^2} + \mathcal{O}(\Delta_x^3) \tag{4}
$$
これを、$\Delta^3$の誤差オーダーと言ったりもする。

また、同様に、次の式も得られる。
$$
u(x-\Delta_x)
= u(x) - \Delta_x \frac{d u(x)}{d x} + \frac{\Delta_x^2}{2!} \frac{d^2 u(x)}{d x^2} + \mathcal{O}(\Delta_x^3) \tag{5}
$$
式(4)と式(5)の二つを足すと、次のような式が得られる。
$$
u(x+\Delta_x) + u(x-\Delta_x)
= 2u(x) + \Delta_x^2 \frac{d^2 u(x)}{d x^2} + \mathcal{O}(\Delta_x^4)
$$
このとき、$\Delta_x$の項と同様に、$\Delta_x^3$のような奇数項全て打ち消してくれるため、誤差オーダーは$\mathcal{O}(\Delta_x^4)$となる。

移行して整理すると、次の式が得られる。
$$
\frac{d^2 u(x)}{d x^2}= \frac{u(x+\Delta_x) -2u(x) + u(x-\Delta_x)}{\Delta_x^2} + \mathcal{O}(\Delta_x^2) \tag{6}
$$
即ち、$u$の二階微分は、右辺を計算することで($\Delta_x^2$オーダーの誤差を持つ値を)求めることができる。
このように$\Delta_x^2$オーダーの誤差の精度を持つ(つまり、$\Delta_x$を半分にすると、その二乗分誤差が小さくなる)ことを二次精度と言う。

二次元テイラー展開

一次元のテイラー展開は有名だが、多次元のテイラー展開は知らない人も多いかもしれない。
これは、次式で表される。

$$
u(x+\Delta_x, y+\Delta_y)=\sum_{p=0}^\infty \frac{1}{p!}\left( \Delta_x \frac{\partial}{\partial x} + \Delta_y \frac{\partial}{\partial y} \right)^p u(x,y)
$$
ただし、このとき、
$$
\left(\frac{\partial}{\partial x} \right)^2 u= \frac{\partial^2 u}{\partial x^2}
$$
のように定義することにする。

例えば$\Delta_x=\Delta_y=\Delta$とすると、次のように展開できる。

\begin{multline}
u(x+\Delta, y+\Delta)=u(x,y)
+ \Delta \left(\frac{\partial u(x,y)}{\partial x}+ \frac{\partial u(x,y)}{\partial y} \right)\\
+ \frac{\Delta^2}{2} \left( \frac{\partial^2 u(x,y)}{\partial x^2} + 2 \frac{\partial^2 u(x,y)}{\partial x\partial y}+ \frac{\partial^2 u(x,y)}{\partial y^2} \right)
+ \mathcal{O}(\Delta^3) \tag{7}
\end{multline}

$x$方向または$y$方向にが一定の場合は一次元の式と同じになる。

\begin{align}
u(x+\Delta_x, y)=\sum_{p=0}^\infty \frac{\Delta_x^p}{p!}\frac{\partial^p u(x,y)}{\partial x^p} \\
u(x, y+\Delta_y)=\sum_{p=0}^\infty \frac{\Delta_y^p}{p!}\frac{\partial^p u(x,y)}{\partial y^p} \tag{8}
\end{align}

よし、準備はできた。
二次元ポアソン方程式のFDMを定式化していこう。

二次精度FDM

離散化

一応、もう一度ポアソン方程式をここに記す。

$$
\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f
$$

式(8)を利用すれば、一次元の二階微分を求めたときのように、次のような式が得られる。

\begin{align}
\frac{\partial^2 u(x+\Delta_x,y)}{\partial x^2}= \frac{u(x+\Delta_x,y) -2u(x,y) + u(x-\Delta_x,y)}{\Delta_x^2} + \mathcal{O}(\Delta_x^2) \\
\frac{\partial^2 u(x,y+\Delta _y)}{\partial y^2}= \frac{u(x,y+\Delta_y) -2u(x,y) + u(x,y-\Delta_y)}{\Delta_y^2} + \mathcal{O}(\Delta_y^2)
\end{align}

これをポアソン方程式に代入すると、次のような離散式が得られる。

\begin{multline}
 \frac{u(x+\Delta_x,y) -2u(x,y) + u(x-\Delta_x,y)}{\Delta_x^2} + \frac{u(x,y+\Delta_y) -2u(x,y) + u(x,y-\Delta_y)}{\Delta_y^2}\\ = f + \mathcal{O}(\Delta_{x,y}^2)
\end{multline}

ここで、$\mathcal{O}(\Delta_{x,y}^2)=\mathcal{O}(\Delta_x^2)+\mathcal{O}(\Delta_y^2)$である。

ここまではなんとなく理解できただろうか。ちょっと式が多いが、言っていることはそんな難しくはないと思うので、ちょっと混乱したらもう一度読み直しておこう。

次に、次の図のように計算領域を格子分割する。

図２：格子分割

各格子点における$u$の値を、$u_{m,n}=u(m\Delta_x, n\Delta_y)$というように定義して、上の離散式に代入して整理すると、次のような式が得られる。
$$
\frac{u_{m-1,n} -2u_{m,n} + u_{m+1,n}}{\Delta_x^2} + \frac{u_{m,n-1} -2u_{m,n} + u_{m,n+1}}{\Delta_y^2} = f_{m,n} + \mathcal{O}(\Delta_{x,y}^2)
$$
ここで$\Delta_x =\Delta_y = \Delta$とすれば、次のように簡単になる。
$$
\frac{u_{m-1,n} + u_{m+1,n} + u_{m,n-1} + u_{m,n+1} - 4u_{m,n}}{\Delta^2} = f_{m,n} + \mathcal{O}(\Delta^2)
$$

この$m$, $n$を整数とすれば、$u_{m,n}$を未知数とした連立方程式を立てることができることに気づける。
この連立方程式を解くことで、$u$の値を求めることができる。

説明のために、gif動画を作った。
以下のように、それぞれの点に注目して式を立てていく。

図３：定式化(思ったより大きくなってしまった……)

ここで、上式の未知数が12個($u_{1,0}, u_{2,0}, u_{0,1}, u_{1,1}, u_{2,1},u_{3,1}, u_{0,2},u_{1,2},u_{1,3},u_{2,1},u_{2,2}$)なのに対して、式はたった４つだ。これでは解けるわけがない。
ちょっと考えてみればわかるが、上のような2×2だけでなく、どんだけ大きくしても、式の数よりも未知数のほうが多くなってしまう。

そこで、境界条件の出番だ。

境界条件

図2のように計算領域をみると、境界上に格子点があることがわかる。
ここもいろいろなやり方があるが、このように格子点を置くのが最も簡単なのでこれで説明しよう。

ディリクレ境界条件

ディリクレ境界条件は次のように表されるのであった。
$$
u = g_1(x,y) \text{ (on }\Gamma_1 \text{ )}
$$
これはそのまま代入すればいい。

ノイマン境界条件

ノイマン境界条件は、一階偏微分だった。
$$
\frac{\partial u}{\partial n} = g \text{ (on }\Gamma_1 \text{ )}
$$
これも、(長方形領域の場合)とても簡単で、テイラー展開で一階微分の近似を行えばいい。
以下のような感じだ。

\begin{align}
\frac{u_{m+1,n} - u_{m,n}}{\Delta_x} = g_2(x,y) \\
\frac{u_{m,n+1} - u_{m,n}}{\Delta_y} = g_2(x,y)
\end{align}

なお、長方形領域のように、$x$座標もしくは$y$座標と平行なら上式で簡単に求められるが、実際の計算領域は斜めになることが多い。この場合は、次回の実践編で説明する。

境界条件の右辺を$g_1(x,y)=0$, $g_2(x,y)=0$として、簡単な図を作ってみた。

図４：境界条件

領域の上下の辺をディリクレ境界条件、左右をノイマン境界条件にしたものである。

これによって、連立方程式の式の数=未知数の数となり、あとはその連立方程式を頑張って解くだけだ。

連立方程式

連立方程式の解き方は……、これはめちゃめちゃいろいろな手法があるのでここで説明はしない。
一応軽く触れておくと、連立方程式は行列を利用して次のように表すことができる。
$$
\boldsymbol{A}\boldsymbol{x} = \boldsymbol{b}
$$
この$\boldsymbol{A}$の逆行列$\boldsymbol{A}^{-1}$を求めて、$\boldsymbol{b}$にかける直接法と、ガウス＝ザイデル法など何度も代入などを繰り返して解に近づけていく反復法がある。
直接法は時間がかかるが安定しており、反復法は発散する可能性があるが高速と、それぞれ一長一短がある。
ちなみにnumpyのnumpy.linalg.solveは直接法のLU分解という方法で解いている。

他の物理量への変換

最初に説明したように、ポアソン方程式から求まる解$u$は電位$\phi$に値する。
この電位$\phi$から電場$E$を求めるにはどうすればいいだろうか？

これも上のようにテイラー展開を利用しよう。
電位と電場は次のような関係があったのだった。
$$
\nabla \phi = \boldsymbol{E}
$$
$x$成分は次のようになる。
$$
\frac{\partial \phi}{\partial x} = E_x
$$
そして、テイラー展開による近似を行う。
$$
\frac{\phi_{m+1,n} - \phi_{m,n}}{\Delta_x} = E_x(m\Delta_x,n\Delta_y)
$$
これは一次精度だが、求める電場を２つの格子点の中心とすれば、二次精度になる。
$$
\frac{\phi_{m+1,n} - \phi_{m,n}}{\Delta_x} = E_x(\left(m+\frac{1}{2}\right)\Delta_x,n\Delta_y)
$$

図５：電場

上図のように、電場を求める際は格子点の間の値を求めるのが一般的である。

まとめ

電界解析はポアソン方程式さえ解ければいいということを説明した。

その解析方法として、今回は有限差分法(FDM)について説明した。
この有限差分法はテイラー展開に基づいており、これさえ理解できれば、あとは連立方程式を立てて、それを解けば解を求めることができる。

次回は、実践編と称して、この差分法を利用して実際にポアソン方程式を解くプログラムを作ってみよう。