読者です 読者をやめる 読者になる 読者になる

定常ラプラス方程式

定常のラプラス方程式は、任意の関数{\phi}について、ラプラシアンをかけて


{\frac{\partial^2 \phi}{\partial x^2} = 0}

なだけ。femの最初の練習には一番ちょうどいい課題ですね。重み付け残差法さえ理解すれば終わりです。このへんはちゃんと教科書に書いてあるからいい加減に書きます。
重み関数をかけて積分形にしてまとめた後、重積分・部分積分ガウスの発散定理がわかる人なら、式展開してると気が付けば二階の微分項が一階の微分項に変わるという。弱形式化はこれで終了。

離散化の際の領域分割ですが、めんどいので二次元の場合は三角形一次要素でしばらく押し通します。愛想パラメトリック要素とか二次要素とかあるけど、めんどいのでしばらく無視します。
形状関数も教科書通りに三元連立方程式を解いて終わり。

さっさと行列方程式の形でまとめて離散化終了。
拡散方程式くらいだとほとんどの教科書にちゃんと書いてあるから楽でした。


連立方程式のソルバーは結構悩む問題かもしれない。中学校の解法みたいにガウスの消去法で解くと時間のかかり方が半端ないし、ガウス・ザイデル法やSOR法のような反復法だってそんなに速くない。共役勾配(CG)法ベースの解法が一番だろうけど、この理論を一から理解するのは辛そうだし…。
妥協して、下の本のアルゴリズムをそのままコードに直す。

反復法の数理 (応用数値計算ライブラリ)

反復法の数理 (応用数値計算ライブラリ)

後々のことを考えると、非対称行列の連立方程式が解ける手法がいいし…ということで、BiCGStab法を実装。本のアルゴリズムの通りにコードを書いて、適当な式解いて合ってるか確認する。この作業、正味30分。何この本超便利。

ここまで来たらやっとコードが書ける。数値計算に慣れてたら、1から勉強を始めてコードの実装まで1、2日で終わるでしょう。むちゃくちゃお手軽です。