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

nano_exit

量子力学、固体物理、fortran、python、etc

逐次代入法を用いた超シンプルな例と収束因子

 x = 1-x

を考える。もちろん答えは x = 0.5である。
これを敢えて逐次代入法で解いてみよう。この方法は固体物理や第一原理計算で用いるSCF法と同じノリである。

最初に初期値を設定する。答えに近ければ近い程良いが、とりあえず x_0 = 0.7とおく。これを右辺に代入したものが左辺を与えるとして、

 x_1 = 1 - x_0 = 0.3
 x_2 = 1 - x_1 = 0.7

と次々に項が求まる。が、よく見ると全然0.5に近づいておらず、このままだとずっと振動して xが求まらない。

そこで収束因子 \alphaを導入する。 \alphaは新しい xをどれくらい古い xに混ぜるかというもので、例えば x_1は次のようになる。

 x'_0 = 1 - x_0
 x_1 = \alpha x'_0 + ( 1 - \alpha) x_0

 \alphaが小さい程、古い xの情報が残り、 xの変化が緩やかになる。
具体的に値を入れてその振舞を見てみると、

 \alpha = 0.2
 x_0 = 0.7
 x'_0 = 0.3
 x_1 = 0.2 \cdot 0.3 + 0.8 \cdot 0.7 = 0.06 + 0.56 = 0.62
 x_2 = 0.2 \cdot 0.38 + 0.8 \cdot 0.62 = 0.072 + 0.496 = 0.568

というように振動が抑えられ、0.5に近づいていくのが分かる。
では本当にこれを続けて0.5になるかどうかを解析的に解いてみよう。

 x_1 = \alpha ( 1 - x_0) + ( 1 - \alpha) x_0 = x_0 (1 - 2\alpha) + \alpha
 x_2 = x_1 (1 - 2\alpha) + \alpha = x_0 ( 1 - 2\alpha)^2 + \alpha (1 - 2\alpha) + \alpha
 \vdots
 x_n = x_0 ( 1 - 2\alpha)^n + \alpha ( 1 + (1 - 2\alpha) + (1 - 2\alpha)^2 + \cdots + (1 - 2\alpha)^{n-1})

よって、これを無限に繰り返すと、 |1 - 2\alpha| < 1 であれば、

 x_{\infty} = \alpha \frac{1}{1-(1-2\alpha)} = \frac{1}{2}

となり、初期値にも収束因子にも依存せず答えがちゃんと求まることが示せた。
 |1 - 2\alpha| < 1 、つまり 0 < \alpha < 1 という条件は、元々が新しい xと古い xの割合ということを思い出せば全く自然な範囲である。
ちなみに、最初に x_0=0.5を代入すると一発で求まるし、 \alpha=0.5でも一発で求まるので、問題によって適切な初期値と収束因子を選ぶことが重要であることがわかる。しかし、一般に収束因子に関しては物理的なものではないので類推は困難であり、経験的に適切だと思われる値を入れることが多いと個人的には思う。
最適な収束因子を求める問題とか数値計算の方でありそうな気がするが、不用意に飛び込むと全ての数値計算法に満足出来なくなりそうで怖いのでやめておこう。。。