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

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

自己紹介      Twitter: @dc1394 所属：放送大学教養学部教養学科4年 マテリアルサイエンス修士 量子力学の数値計算をやっています． 最も興味のある分野 第一原理計算  密度汎関数理論（Density Functional Theory, DFT）   第一原理計算やDFTについては，よろしければ以 下URLの拙作のスライドやQiitaの記事をご覧くだ さい． https://www.slideshare.net/dc1394/ss-26378208  https://qiita.com/dc1394/items/1bb597b27d7511c e31ba 

• https://www.slideshare.net/dc1394/ss-26378208
• https://qiita.com/dc1394/items/1bb597b27d7511c

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

Thomas-Fermi方程式とは   Thomas-Fermi方程式（以後T-F方程式と呼ぶ） とは，以下の非線形常微分方程式（導出は割愛） である． この方程式に解析解はない．   数値解を求める必要がある． 境界条件はy(0) = 1, y(∞) = 0 数値的に解く際には上記は無意味である．  y(xmin), y(xmax )の二点境界値問題にしてむりやり解く． 

常微分方程式の二点境界値問題  常微分方程式の二点境界値問題の数値的解法 の種類は以下である． 狙い撃ち法  有限差分法  有限体積法  有限要素法   狙い撃ち法と有限要素法の組合せでT-F方程式 を解く．

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

有限要素法の例    右図は有限要素法の 例． 赤線が厳密解，青線 が有限要素法で求め た数値解． 厳密解にそって区分 線形近似されている ことが分かる． https://qiita.com/Sunset_ Yuhi/items/4c4ddc25609 a7619cce0 より引用．

• https://qiita.com/Sunset_

原点に十分近い点での境界値y(xmin)  y(x)の展開式を用いて，原点に十分近い点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.

遠方での境界値y(xmax)    次に遠方xmaxでの境界値を決定する．論文[1]に よれば，遠方で漸進的に， とある．ここで，λ = 3.886, x0 = 5.2415であ る． [1] M. Desaix, D. Anderson, and M. Lisak, Eur. J. Phys. 25, 699 (2004).

T-F方程式の反復による解法  (1)式を直接，有限要素で離散化するのは難しい．  しかし(2)式なら，有限要素で離散化するのは容易である．   以下のような繰り返しによりT-F方程式を解く． 初期関数y0(x)はどうするか？  狙い撃ち法で求める．

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

• http://www.slis.tsukuba.ac.jp/

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

適合点への狙い撃ち法の問題   適合点の値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

要素方程式  最も単純な形状関数である1次要素によって，T-F 方程式を離散化すると，要素方程式は以下のよう になる．  ただし，Aeの要素Aij(i, j = 1, 2)は，  であり，beの要素bi(i = 1, 2)は，   である．ここで，x0e, x1eはe番目の要素をなす2節 点のx座標，l = x1e － x0eは線分要素の長さである． なおbeは，数値積分により求めなければならない．

要素方程式から全体行列と全体ベ クトルを生成する    全体行列と全体ベクトルを生成するには，全 体行列の対角要素（最初と最後を除く）が隣 接する要素からの寄与の合計であることに注 意して，要素ごとに上式を重ね合わせること で可能である． 最終的に，[K]{u} = {f}という連立一次方程式 に帰着する． ここで，[K]は全体行列，{f}は全体ベクトルで ある．

全体行列と全体ベクトル  全体行列[K]と，全体ベクトル{f}は，次式のよ うになる（簡単のために[K]が4次正方行列の 場合を考える） ． [K] {u} {f}

境界条件処理   連立一次方程式[K]{u} = {f}を数値的に解く前 に， 全体行列[K]と全体ベクトル{f}に境界条 件処理を施す必要がある． 境界条件にはDirichlet境界条件とNeumann境 界条件があるが，[K]と{f}に施す必要があるの は前者である．

Dirichlet境界条件の組み込み  Dirichlet境界条件を組み込む前の式は以下で ある（簡単のために[K]が4次正方行列の場合 を考える）． [K] {u} {f}

Dirichlet境界条件の組み込み  Dirichlet境界条件を組み込んだ後の式は以下 のようになる． [K] {u} {f}

[K] をメモリに格納するための工夫  1次要素だと，連立1次方程式[K]{u} = {F}の [K]が対称三重対角行列になるので，[K]は2次 元配列を使わなくても，単に配列を2つ用意す るだけで済む．   これはメモリの節約になる． Juliaでは，対称三重対角行列を格納するため に，SymTridiagonalという専用の型が存在す る．

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

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

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

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

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

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

Thomas-Fermi方程式の解のグラフ （y軸対数目盛）

実行画面

インプットファイル

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

並列化    このプログラムのボトルネックは，初期関数 y0(x)とbeを求める部分である． 前者は並列化できないが，後者は並列化でき る． 従って，この部分をプロセス並列化あるいは スレッド並列化できれば，プログラムの高速 化が見込める．

スレッド並列化       著者のコードではスレッド並列化で計算速度の向 上を図っているが，その効果を検証したのが次 ページの表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.2 - Cinnamon (64-bit) 物理メモリ:32GB (8GB×4, DDR4-3200)

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

課題  この方法で得たy(x)は，精度がそれほど良くな い．   Lobatto多項式のような，より複雑な形状関数を 用いる必要性 反復回数が多すぎる．  Broyden法やPulay法などの，より高度な混合法 を用いる必要性

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

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

• https://github.com/dc1394/thomasfermi_julia

参考サイト  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

• https://qiita.com/Sunset_Yuhi/items/4c4ddc25
• http://www.fem.gr.jp/index.html