Juliaで解くThomas-Fermi方程式
放送大学教養学部教養学科
川井弘之

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より引用．

y(xmin)の推定

論文[1]によれば，y(x)の展開式は，である．
 ここで，B = y’(0)であり，同論文によれば，B≒ -1.588076779である．
 上式により，原点に十分近い点xminでの境界値y(xmin)が求められる．


[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)の推定

論文[1]によれば，y(x)は遠方で漸進的に，である．ここで，λ = 3.886, x0 = 5.2415である．
 上式により，遠方の点xmaxにおけるy(xmax)が求められる．


[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/mic/text/09_derivative2.pdf
より引用

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

最も単純な形状関数である1次要素によって，T-F方程式を離散化すると，要素方程式は以下のようになる．

ただし，Aeの要素Aij(i, j = 1, 2)は，であり，beの要素bi(i = 1, 2)は，である．ここで，x0e, x1eはe番目の要素をなす2節点のx座標，l = x1e － x0eは線分要素の長さである．
e なおb は，数値積分により求めなければならない．
要素方程式から全体行列と全体ベクトルを生成する全体行列と全体ベクトルを生成するには，全体行列の対角要素（最初と最後を除く）が隣接する要素からの寄与の合計であることに注意して，要素ごとに上式を重ね合わせることで可能である．
 最終的に，[K]{u} = {f}という連立一次方程式に帰着する．
 ここで，[K]は全体行列，{f}は全体ベクトルである．


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

[K]

{u}

{f}

境界条件処理

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

Dirichlet境界条件の組み込み


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

[K]

{u}

{f}

Dirichlet境界条件の組み込み


Dirichlet境界条件を組み込んだ後の式は以下のようになる．

[K]

{u}

{f}

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

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

課題


この方法で得たy(x)は，精度がそれほど良くない Lobatto多項式のような，より複雑な形状関数を用いる必要性

反復回数が多すぎる Broyden法やPulay法などの，より高度な混合法を用いる必要性

まとめThomas-Fermi方程式を，狙い撃ち法および有限要素（1次要素）による離散化と反復法によって，数値的に解いた．