nano_exit

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

逐次代入法を用いた最もシンプルな例と収束因子:図のテスト

前回、簡単な例で逐次代入法を考察した。koideforest.hatenadiary.com

図の使い方のテストもかねて、その内容を図にしてみた。

収束因子を使わない場合、 x = 1-xは以下のようになる。
f:id:koideforest:20161228135911p:plain
ただの長方形をグルグルまわるだけで一向に収束しないのが分かる。

収束因子 \alpha = 0.7を使うと、
f:id:koideforest:20161228135918p:plain
のように y = x y = 1-xの交点にドンドン近づいていくのが分かる。

ちなみに、これらはgnuplotで出力しているが、直線の連続はデータ点を読み込ませていて、 y = x y = 1-xgnuplotをそのまま使っている。つまり、

pl 'iteration.dat' u 1:2 w l
rep x
rep 1-x

と入力している。
各ファイルの中身は、

iteration.dat:
-0.50000 -0.50000
-0.50000 1.50000
1.50000 1.50000
1.50000 -0.50000
-0.50000 -0.50000
-0.50000 1.50000
1.50000 1.50000
1.50000 -0.50000
-0.50000 -0.50000
-0.50000 1.50000

iteration_broyden.dat:
-0.50000 -0.50000
-0.50000 1.50000
0.90000 0.90000
0.90000 0.10000
0.34000 0.34000
0.34000 0.66000
0.56400 0.56400
0.56400 0.43600
0.47440 0.47440
0.47440 0.52560

である。
これらのデータはpythonで計算させて出力させていて、収束因子を使った方は、

#!/usr/bin/env python

def fn( x ):
    return 1.0 - x

f = open( 'iteration_broyden.dat', 'w' )

x     = -0.5
alpha = 0.7
for i in map( float, range(5) ):
    f.write( '{:10.5f} {:10.5f}\n'.format( x, x       ) )
    f.write( '{:10.5f} {:10.5f}\n'.format( x, fn( x ) ) )
    x = ( 1.0 - alpha ) * x + alpha * fn( x )

f.close()

というような感じになっている。