Thomas-Fermi方程式のFEMによる解法

スライド概要

Thomas-Fermi方程式を、有限要素法(Finite Element Method,FEM)でどう解くかを解説したスライドです。

profile-image

dc1394

@dc1394724

作者について:

スライド一覧
シェア
埋め込む

作成日

2021-08-01 22:52:22

各ページのテキスト

1. Thomas-Fermi方程式のFEMによる解 法 @dc1394

2. Thomas-Fermi方程式とは  Thomas-Fermi方程式とは,密度汎関数理論 における局所密度近似において,交換エネル ギー項を無視した際に得られる全エネルギー の式 運動エネルギー  (電子-核間 の)Coulomb エネルギー 電子-電子間のCoulomb エネルギー (Hartreeエネルギー) を,原子に適用した際に現れる非線形常微分 方程式であり,密度汎関数理論において重要 な方程式である.

3. Thomas-Fermi方程式とは    Thomas-Fermi方程式(以後T-F方程式と呼 ぶ)とは,以下の非線形常微分方程式である. T-F方程式は,解析的には解けないので,数値 的に解くしかない. ここでは,T-F方程式を,Julia言語を用いて, 有限要素法で数値的に解くことにする.

4. 有限要素法とは    有限要素法(Finite Element Method, FEM) とは,数値解析手法の一つであり,解析的に 解くことが難しい常微分・偏微分方程式の近 似解を数値的に得る方法の一つである. 方程式が定義された領域を小領域(要素)に 分割し,各小領域における方程式を比較的単 純で共通な補間関数で近似する. 構造力学分野で発達し,他の分野でも広く使 われている手法である.

5. 有限要素法の例    右図は有限要素法の 例. 赤線が厳密解,青線 が有限要素法で求め た数値解. 厳密解にそって区分 線形近似されている ことが分かる. Pythonで1次元有限要素法 (Poisson方程式)- Qiita より引用.

6. Thomas-Fermi方程式の数値解法上 の困難     T-F方程式を改めて書くと, であり,境界条件は,y(0) = 1, y(∞) = 0である. この境界条件は,この方程式を数値的に解く上で 何も言っていないのと同じである. この方程式は原点に特異点を持つので原点から解 けない.また当然ながらコンピュータは無限大を 扱えない.

7. 原点に十分近い点での境界値  仕方がないので,y(x)のホモトピー摂動法による 展開式[1]を用いて,原点に十分近い点xminでの境 界値を決定する. 論文[1]によれば,y(x)の展開式は,  であり,またy(x)の微分y’(x)は,     である.ここで,B = y’(0)であり,同論文によれ ば,B ≒ -1.588076779である. なお,このy(xmin)とy’(xmin)を使って,T-F方程式 を初期値問題として数値的に解こうとしても,遠 方で発散する. [1] M. A. Noor and S. T. Mohyud-Din. Homotopy Perturbation Method for Solving Thomas-Fermi Equation Using Pade Approximants. International Journal of Nonlinear Science, 8(1)(2009): 27-31.

8. 遠方での境界値      次に遠方xmaxでの境界値を決定する.文献[1]によ れば遠方で漸進的に, と書いてあるが,この近似は粗すぎる.ここで, 論文[2]によれば,やはり遠方で漸進的に, とある.ここで,λ = 3.886, x0 = 5.2415である. [1] R.G.パール, W.ヤング 『原子・分子の密度汎関数法』シュプリンガー・フェアラーク東京 (1996) [2] M. Desaix, D. Anderson, and M. Lisak, Eur. J. Phys. 25, 699 (2004).

9. T-F方程式の問題     これで,原点にごく近い点xminでの境界値と, 遠方xmaxでの境界値が分かったことになる. つまり,常微分方程式の2点境界値問題に帰着 したので,次に有限要素による離散化を考え てみる. T-F方程式は非線形常微分方程式なので,離散 化にも工夫が必要である. ここで,以下のような線形常微分方程式につ いて考えてみよう.

10. β(x)を用いた反復法    この方程式は線形常微分方程式であり,有限要素 による離散化によって,[K]{u} = {F}という形式 の連立1次方程式に帰着でき,解ける. ここでもし,β(x)を次のように置くならば, (1)式を解くことによって得たy(x)から,(2)式に よってβ(x)を得て,それからまた(1)式を解くこ とによってy(x)を得て…という反復法によって解 ける(T-F方程式を直接離散化するなどの別の方 法もあるだろうが,ここではこの方法を採用す る).

11. 初期関数y0(x)    この方法でT-F方程式を解くには,y(x)の初期 関数y0(x)が必要である. ここで,まずT-F方程式を適合点への狙い撃ち 法で解き(常微分方程式の解法には,予測子 修正子法の一種であるAdams-Moulton法を用 いる),その解をy0(x)とする. なお,Adams-Moulton法は,Juliaの DifferentialEquations.jlに実装されているので, これを使う.

12. 狙い撃ち法とは   狙い撃ち法とは,常微 分方程式の境界値問題 の数値解法の一つであ る. 常微分方程式の,初期 値問題と境界値問題の 違いは,後者は始点だ けでなく終点も指定さ れることである.  つまり,初期値問題と 同じ方法でfnまで求めて, 境界条件から外れてい たら,逐次修正すれば 良い. http://www.slis.tsukuba.ac.jp/ ~fujisawa.makoto.fu/lecture/m ic/text/09_derivative2.pdf より引用

13. 適合点への狙い撃ち法の問題   適合点の値xmatchは任意であるが,概ねxmatch = 10.0~15.0程度でないと収束しない. なお,適合点での狙い撃ち法は,適合点で値 が「異なる」ので,精度の高い解が得られな いが,初期関数y0(x)とするには十分である. x xminから解いたy(x) 9.998 0.0243354 9.999 0.0243308 10.000 (xmatch) 0.0243262 xmaxから解いたy(x) 0.0244253 10.001 0.0244206 10.002 0.0244160

14. 1次混合    入力と出力が一致する解を得るために,最も 簡単な1次混合法を使う. すなわち,i + 1段階での改善された入力関数 yi+1in は,i段階でのyiinとyioutを用いて,次式 で与えられる. ここでαは任意のパラメータであるが,α = 0.1程度としないと収束しない.

15. 反復法の収束の判定     yiは(数値計算上は)ベクトルと見なせる. ここで,yioutとyiinの差のベクトルの大きさを NormRDと適当に名付け, NormRD = |yiout – yiin|と定義する. このNormRDがある閾値(criterion)未満に なった場合(NormRD < criterion),収束し たと判定する. ここで,criterionは任意のパラメータである が,10-13程度が限界のようである.

16. 有限要素による離散化   最も単純な1次要素によって,T-F方程式を離 散化する(2次要素は,試してみたが収束しな かった). 1次要素だと,連立1次方程式[K]{u} = {F}の [K]が対称正定値三重対角行列になるので, [K]は2次元配列を使わなくても,単に配列を2 つ用意するだけで済む.  これはメモリの節約になる.

17. フローチャート 初期関数y0(x), [K]の生成 β(x), {f}の生成 [K], {f}に境界条件処理 連立方程式を解きy(x)を求める 収束したか? YES 終了 NO

18. T-F方程式の数値解とエネルギー      以上の解法で得たT-F方程式の数値解は,文献 [1]に載っている,T-F方程式の解の数表と「ほ ぼ」一致している. T-F方程式の解から得られる中性原子のエネル ギーは,Zを原子番号として, である. [1] E. U. Condon, Halis Odabasi. Atomic Structure, Cambridge University Press, Cambridge, 1980 ちなみにこの本はGoogle ブックスで(全部ではないが)読める.数表のペー ジのURLはこちら: http://bit.ly/1fRQ71T

19. T-F方程式の数値解とエネルギー    ここで,y’(0)は,有限要素法と補外によって得た 値を用いる. 著者のコードでは,y’(0) = -1.5459(エネルギー の値は,Zを原子番号として,E = -0.7483Z7/3 (Hartree))であり,この値は文献[1]の値とまあま あ一致している. [1] E. U. Condon, Halis Odabasi. Atomic Structure, Cambridge University Press, Cambridge, 1980

20. Thomas-Fermi方程式の解のグラフ

21. Thomas-Fermi方程式の解のグラフ (y軸対数目盛)

22. 実行画面

23. インプットファイル

24. 並列化による高速化    昨今のCPUには,少なくとも2つ以上のCPUコ アが内蔵されており,この複数のCPUコアを 同時に使うことで,処理を高速化できる. ただし,複数のCPUコアを同時に使うために は,プログラムのコードを適宜書き換える必 要がある(自動ではやってくれない!). なお,複数のCPUコアを同時に使う(並列化 する)ためには,マルチスレッドによる並列 化とマルチプロセスによる並列化がある.

25. スレッド並列化の環境        著者のプログラムでは,マルチスレッドによる並 列化(スレッド並列化)を行い,計算速度の向上 を図っている. その効果を検証したのが次ページの表1である. 計測環境: CPU: Intel Core i9-10980XE (Cascade Lake-X, Hyper Threading ON (物理18コア論理36コア), Enhanced Intel SpeedStep Technology OFF, Turbo Mode OFF) Juliaのバージョン: v1.7.0-beta3 OS: Linux Mint 20.1 - Cinnamon (64-bit) 物理メモリ:32GB (8GB×4, DDR4-3200)

26. スレッド並列化のパフォーマンス 表1 並列化の有無の効果の比較(三回の平均) 計算時間(秒) 並列化効率 収束時のNormRD スレッド並列化なし 1001.6 1 9.3×10-13 スレッド並列化あり 56.1 17.9 9.3×10-13 ※計算時間とは,総経過時間から,結果出力処理の時間を差し引いた時間を 差す   並列化により,17.9倍のパフォーマンス向上 に成功していることが分かる. 並列化による誤差は確認できない.

27. ソースコードへのリンク    このプログラムのソースコードは,GitHub上 で公開しています. https://github.com/dc1394/thomasfermi_juli a ライセンスは2条項BSDライセンスとします.

28. まとめ      Thomas-Fermi方程式を,有限要素(1次要 素)による離散化と反復法によって,数値的 に解いた. 得られた結果は,文献の結果とほぼ一致して いた. 高速化をめざして,マルチスレッドによる並 列化(スレッド並列化)を試みた. スレッド並列化は非常に効果的であった. 並列化による誤差は確認できなかった.

29. 参考サイト  Pythonで1次元有限要素法(Poisson方程式) - Qiita  https://qiita.com/Sunset_Yuhi/items/4c4ddc25 609a7619cce0  Internet-College of Finite Element Method  http://www.fem.gr.jp/index.html