自己無撞着方程式
数値計算で方程式を解く簡単な方法の１つとして逐次代入法を紹介します。
x=f(x)
という方程式を解くことを考えます。
逐次代入法
ニュートン法など色々な方法があると思いますが、最も簡単な方法の一つが逐次代入法です。適当な初期値$x^0$から初めて
\begin{align}
x^{i+1} &= f(x^i) \\
\text{STOP if }& |x^{i+1}-x^i|<\epsilon
\end{align}
と更新していく方法です(もちろん$f$の種類によっては収束しない場合もあります)。
この方法は簡単ですが一般に収束が遅いことでも知られています。そこでいくつか収束を早める方法があります。
割線法
抵当な定数$\beta \in (0,1] $を用いて
x^{i+1} = (1-\beta)x^i + \beta f(x^i)
と更新する方法です。これによって変化分が内分されて収束点周りの振動を抑えることができます。
アンダーソン法(アンダーソンミキシング)
各ステップ$i$に対して定数$\beta^i = \in (0,1] $, と整数$m^i \in [1,i]$を用いて
\begin{align}
x^1 &= f(x^0) \\
x^{i+1} &= (1-\beta^i)\sum_{j~0}^{m^i} \alpha_j^i u^{i-m^i+j}+ \beta^i \sum_{j~0}^{m_i} \alpha_j^i f( u^{i-m^i+j} ) \\
,where& \\
\min_{ \sum_{j=0}^{m^i}\alpha_j^i = 1 }& \left|\sum_{j=0}^{m^i} \alpha_j^i F(x^{k-m^i+j}) \right|^2 \\
F(x) &= f(x) - x
\end{align}
と更新します。これによって更新後の内分点が最小になるように自動調整する方法です。
例
例えば、$f$が実関数で$\beta^i=1,m^i=1$であるとします。
すると更新の式は
u^{i+1} = \alpha^i_0 f(x^{i-1}) + \alpha^i_1 f(x^i)
各ステップで$\alpha$を求めるためにはラグランジュの未定乗数法を使って
g(\alpha_0^i,\alpha_1^i,\lambda) = [ \alpha^i_0 F(x^{i-1}) + \alpha^i_1 F(x^i) ]^2 -\lambda( \alpha^i_0 + \alpha^i_1 -1 )
を最小化します。そのための停留条件は
\frac{\partial g}{\partial \alpha_0^i} = 0 \\
\frac{\partial g}{\partial \alpha_1^i} = 0 \\
\frac{\partial g}{\partial \lambda} = 0
の3つであり、この連立方程式を解くと
\alpha_0^i = \frac{F(x^i)}{ F(x^i) - F(x^{i-1}) } \\
\alpha_1^i = \frac{F(x^{i-1})}{ F(x^{i-1}) - F(x^i) }
となります。ちなみに $f:\mathbb{C}^N \rightarrow \mathbb{C}^M$のときは
\begin{align}
g(\alpha_0^i,\alpha_1^i,\lambda) &= ( \alpha^i_0 F^*(x^{i-1}) + \alpha^i_1 F^*(x^i) )^T \cdot ( \alpha^i_0 F(x^{i-1}) + \alpha^i_1 F(x^i) ) \\
&-\lambda( \alpha^i_0 + \alpha^i_1 -1 )
\end{align}
を最小化して
\begin{align}
\alpha_0^i &= Re\left[\frac{ F^*(x^{i})\cdot\{ F(x^i)-F(x^{i-1}) \} }{ |F(x^i)-F(x^{i-1})|^2 }\right] \\
\alpha_1^i &= Re\left[\frac{ F^*(x^{i-1})\cdot\{ F(x^{i-1})-F(x^{i}) \} }{ |F(x^i)-F(x^{i-1})|^2 }\right]
\end{align}
となります。
実装
ここでは試しに
x = f(x) = 0.1 x^5 - 0.08x^4+0.5x^3-10x^2+3x
という方程式の解の一つを求めてみます。$y=f(x)$は次図の青線のような曲線です。オレンジの線は$y=x$です。
初期値は$x_0=0.2$としました。アルゴリズムの性質上、$y=x$より上に出ている部分からスタートすれば右隣、下なら左隣の交点の$x$の値が解として求まります(もし、そのような交点がなければ計算は発散します)。
A：単純な逐次代入法、B：割線法、C：アンダーソン法でそれぞれ計算してみます。
import numpy as np
from matplotlib import pyplot as plt
def f(x):
return 0.01*x*x*x*x*x - 0.08*x*x*x*x + 0.5*x*x*x-10*x*x + 3*x
def updataA(x0):
xin=x0
log = []
for i in range(1000):
x=f(xin)
err = abs(x-xin)
if(err<1e-8): break
log += [np.log10(err)]
xin=x
return x, log
def updataB(x0):
b=0.3
xin=x0
log = []
for i in range(1000):
x=b*f(xin) + (1.-b)*xin
err = abs(x-xin)
if(err<1e-8): break
log += [np.log10(err)]
xin=x
return x, log
def updataC(x0):
b=0.3
xin=f(x0)
xin2=x0
log = []
for i in range(1000):
xn = f(xin)
d1 = xin-xin2
d2 = xn-xin
x= (d1/(d1-d2))*xn + (d2/(d2-d1))*xin
err = abs(x-xin)
if(err<1e-8): break
log += [np.log10(err)]
xin2=xin
xin=x
return x, log
次のように実行結果を出力すると
x0 = 0.2
x1,log1 = updataA(x0)
x2,log2 = updataB(x0)
x3,log3 = updataC(x0)
plt.plot(log1,label="simple")
plt.plot(log2,label="divid")
plt.plot(log3,label="Anderson")
plt.xlim([0,20])
plt.legend()
最終的な解は
\begin{align}
A:& 0.20197545803632116 \\
B:& 0.2019754464821037 \\
C:& 0.20197545459296495
\end{align}
かかったステップ数は
\begin{align}
A:& 666 \\
B:& 13 \\
C:& 5
\end{align}
です。誤差の変化をプロットすると次のようになります($\log_{10}$スケール)。
確かにアンダーソン>割線>何もしないの順に収束が早くなっています。