#はじめに
数値計算関連の勉強をし直している中、知識の整理がてら、記事を投稿してみました。
少し古めの本ですが、加川幸雄先生らの『電気・電子境界要素法―基礎と応用』にあった
正則化境界要素法にチャレンジしています。
#正則化境界要素法について
境界要素法というのは、物理シミュレーション（偏微分方程式を数値的に解く作業）の一つです。
この手のシミュレーションは有限差分法や有限要素法といったものがメジャーな印象がありますが、
それとは少し毛色が違うものとして境界要素法が位置付けられているようです。
具体的には、解きたい偏微分方程式を変形するにあたり、
境界についての積分の形にするのが、境界要素法で、
領域についての積分の形にするのが、有限要素法になります。
（この二つの比較については別サイトや資料を確認してもらえると幸いです）
もちろん例外はありますが、通常、境界要素法を考える場合は境界そのものをいくつもの小さな"要素"に分けて解析をすることになります。
ただ、この少しマイナーな境界要素法には、一つ大きな弱点があります。
それが、要素同士が近すぎる場合に計算が不安定になってしまうということです（すっごい乱暴な言い方）。
そこで、先人たちはどうにかして、計算にかかわる要素間の距離を大きくして計算できないかと考えました。ここで出てきたのが正則化境界要素法とよばれる手法です。
雑に言ってしまえば、「境界」だけですべて完結させるのではなく、その境界で囲われた領域の外に「仮想電荷」を置いて、この不安定な計算を回避するというのが概要になります。
ここまでの話をまとめるとこんな感じになります。
#対象とする問題
理系学生おなじみ、定期テストで出てくる円筒ポテンシャルのラプラス問題をやってみました。
対称性とかを利用して、四分円の形で解析しています。
#計算の流れ
###関係式の記述
本来なら導出とかについても書くべきなんですが、
チートシートよろしく、答えだけゴリゴリ載せておきます。
次のように定式化された問題に対して、
\left\{
\begin{array}{lll}
∇^2 ϕ=0 & on \,\, Ω \\
ϕ=\hat{ϕ} & on \,\, Γ_1 \\
\frac{∂ϕ}{∂n} = \hat{q} &on \,\, Γ_2
\end{array}
\right.
正則化境界要素法においては、境界上（に設置した要素）における解が次の関係を満たします。
$$ [K]\bigr\{\tilde{ϕ}\bigr\}-[G]\bigr\{\tilde{q}\bigr\}=\bigr\{0 \bigr\} $$
[K]≡[F][R] \\
[R]≡[[G]^{-1} ]^T [L]
ただし、$[ ]$はマトリックス、$ \bigr\{ \bigr\} $はベクトルを意味しています(原著の表記)。
なお、各マトリックスは次のように定義されています。
[G]_{jm}≡\int_{Γ_m}^{} ϕ_{jm}^* dΓ\\
[F]_{jk}≡\frac{1}{2} \sum_{m=1}^{N} \int_{Γ_m}^{} ϕ_{jm}^* \frac{∂ϕ_{km}^*}{∂n}+ ϕ_{km}^* \frac{∂ϕ_{jm}^*}{∂n} dΓ \\
[L]_{jk}≡ δ_{mn} S_m
ここで、$n, m = (1, 2,..., N)$は要素に関する添え字、
対して、$j, k = (1, 2,..., M)$は仮想電荷に関する添え字を意味しています。
また、$δ_{ij}$はディラックのデルタ関数を表しています。
\left\{
\begin{array}{lll}
δ_{ij}=1 & (i=j) \\
δ_{ij}=0 ̂& (i \neq j) \\
\end{array}
\right.
ここでは、2次元のラプラス問題を解いているため、
グリーン関数の具体的な形として次に示す基本解を代入します。
ϕ_{jm}^*=-\frac{1}{2π}ln\,{r_{jm}} \\
\frac{∂ϕ_{jm}^*}{∂n} = \frac{-n_{mx}X_{jm}-n_{my}Y_{jm}}{2πr_{jm}^2} \\
r_{jm}≡\sqrt{X_{jm}^2+Y_{jm}^2 } \\
X_{jm}≡(x_j-x_m ), \,\,\,\, Y_{jm}≡(y_j-y_m )\\
なお、$n_{mx}$および$n_{my}$は$n$番目の(線)要素に垂直なベクトルのx成分, y成分を表しています。
###境界条件の反映と連立一次方程式の記述
境界要素ごとの境界条件を確認すると、ポテンシャルそのものが境界条件として与えられている要素、
もしくはフラックス（ポテンシャルの勾配）が既に与えられている要素の二つがあることが確認できます。
この場合、例えば、ポテンシャルが固定されていたら、フラックスは未知になります（逆もしかり）。
そこで、ポテンシャルのベクトル$\bigr\{\tilde{ϕ}\bigr\}$を既知のもの$\bigr\{\hat{ϕ}\bigr\}$、未知のもの$\bigr\{\check{ϕ}\bigr\}$にわけ、
フラックス$\bigr\{\tilde{q}\bigr\}$についても同様に分けていきました。
したがってメイン関係式はまず、次のようになります。
[K]\bigr\{\tilde{ϕ} \bigr\}-[G]\bigr\{\tilde{q}\bigr\}=\bigr\{0 \bigr\} \\
[K]\begin{bmatrix} \hat{ϕ} \\ \check{ϕ} \end{bmatrix}
-[G]\begin{bmatrix} \check{q} \\ \hat{q}\end{bmatrix}=\bigr\{0 \bigr\} \\
すると、マトリックスについても分割できることがわかります。すなわち、
[\hat{K}|\check{K}]\begin{bmatrix} \hat{ϕ} \\ \check{ϕ} \end{bmatrix}-[\check{G}|\hat{G}]\begin{bmatrix} \check{q} \\ \hat{q}\end{bmatrix}=\bigr\{0 \bigr\} \\
なお、この式は次のように書き換えることができます。
[\hat{K}|-\hat{G}]\begin{bmatrix} \hat{ϕ} \\ \hat{q} \end{bmatrix}-[\check{G}|-\check{K}]\begin{bmatrix} \check{q} \\ \check{ϕ}\end{bmatrix}=\bigr\{0 \bigr\} \\
したがって解くべき連立方程式は次のようにまとめられます。
$$ A\boldsymbol{x} = \boldsymbol{b} $$
\left\{
\begin{array}{lll}
A=[\check{G}|-\check{K}] \\
\boldsymbol{x} =\begin{bmatrix} \check{q} \\ \check{ϕ}\end{bmatrix} \\
\boldsymbol{b} =[\hat{K}|-\hat{G}]\begin{bmatrix} \hat{ϕ} \\ \hat{q} \end{bmatrix}
\end{array}
\right.
これを解くことで、境界(要素)上の値を求めることができます。
###領域内部の値について
最後に、領域内部の求めたい位置におけるポテンシャルを$ϕ_{i}$は、
受信点からj番目の仮想電荷までの距離を$r_{ji}$
連立方程式を解くことで導かれた境界上のポテンシャルを$\bigr\{\tilde{ϕ}\bigr\}$として、
次のように求められます。
ϕ_i=\bigr\{{ϕ_i^*} \bigr\}^T[R]\bigr\{\tilde{ϕ}\bigr\} \\
\bigr\{{ϕ_i^*} \bigr\}_j=-\frac{1}{2π}ln\,{r_{ji}} \\
#実装
####動作環境
- OS: Ubuntu 20.04.2 LTS
- 言語: C++
- コンパイラ: gcc version 9.3.0
- その他ライブラリ: eigen 3.3.9 (行列計算に利用)
- コード
####動かし方
- Makefileを編集(必要に応じて)
- コンソールにてmakeと入力
- 作成したアプリを呼び出せば行けるは
- 要素(端)と仮想電荷, 領域内の点におけるポテンシャルがcsvで出力される
Makefileはこんな感じ(きたない)
CC = g++
TARGETS = clean
FILES := geometry_common geometry_sample integralGL general_func core_matrix solve main
Appname = app
LIBS = -I/usr/include/eigen3
MAIN: $(Appname) clean
@printf "#\n# the App has been built!!! \n#\n"
$(Appname): CPP
$(CC) -o $@ $(addsuffix .o,$(FILES))
CPP:
$(CC) -c $(addsuffix .cpp,$(FILES)) -std=c++17 $(LIBS)
clean: clean_msg
rm -f *.o
clean_msg:
@printf "#\n# cleaning obj files: \n#\n"
環境によって適宜変更してみてください。
#結果
####テスト1
ここで、先ほどの図を再度載せておきます。
こんな感じで要素や仮想電荷を設置して($a=2,, b=5,, ϕ_a = 4,, ϕ_b=15,, r=2,, c=2$)
計算した結果はこうなります(計算したのは左、厳密解が右)。
なお、厳密解として以下のものを右にプロットしています。
これだけだと違いが判らないので、誤差を表示するとこうなります。
内側になるにつれて、誤差が増えているようです。
####テスト2
次に、仮想電荷を外側ギリギリまで近づけて($a=2,, b=5,, ϕ_a = 4,, ϕ_b=15,, r=1,, c=2$)、
計算するとこうなります。
若干嫌な予感がするので、誤差を出力してみると、
このように外側の解析に支障が出ていることがわかります。
#注意点
####文献の誤植？
中古で本を買ったせいもあると思いますが、
行列の表現が怪しい部分がありました。
文献『電気・電子境界要素法―基礎と応用』で次のように書かれているところを
$$[K]≡[R]^T[F][R] $$
下のように直すと、正しい解が得られます。
$$[K]≡[F][R] $$
本を読んでる最中におかしくね？とは思いましたが、
改めて実装してみたら、やっぱり間違っていたようでした。
####仮想電荷の設置について
仮想電荷についてはピシッと四分円の形に整列させないで設置するのがセオリーのようです。
整列させてしまうと、連立方程式の係数行列が正則でなくなり、計算ができなくなるとのこと。
その代わり、1点だけ（本当に）ビミョーにずらすことで、計算できました。
↑こんな感じで仮想電荷の位置を少しずらしています。
####マトリックスの導出に用いる境界積分について
ガウス・ルジャンドル積分(3点)で近似しています。
#追記(21/02/24)
発散を回避する手段について、アドバイスをいただきました。
マトリックスに単位行列の微小定数倍を足して正則行列にすることで達成されるようです。
氏の指摘では、連立方程式の係数行列にこれを施せばよいとのことでしたが、
それ以前に逆行列を計算する部分があったため、
そちらに修正をかけました。
具体的には、$[K]$を計算する前の、$[G]$を次のようにしています。（関係式の記述を参照）
$$[G]\leftarrow[G]+ε[I]$$
(なお、$ε=10^{-10}$としています)
この方法を用いれば、このように規則正しく外部電荷を並べても、
こんなふうに計算できていることがわかります。
真摯にアドバイスをくださったNaokiHori氏には多大なる感謝の意を表します
#雑記・コメント
- 境界要素法に関する情報が少なすぎて焦った（大学の講義でもあまり耳にしないレベルでマイナーな気がする、面白かったけど）
- 境界要素法についてはPythonのオープンソースソフトウェアがあるので、本当に境界要素法が気になる方はそっちをいじった方がいいかも
- コードについて改善点とかあったら、些細な事でも指摘してほしいです
- 初投稿だったので、この記事やこの記事には大変お世話になりました。
#参考資料
- 加川幸雄, 榎園正人, 武田毅. 電気・電子境界要素法 : 基礎と応用 / 加川幸雄, 榎園正人, 武田毅共著. 東京, 森北出版, 2001, 気になった方はこちら