概要
有限要素法(FEM)や粒子法で非線形方程式を解く際、ニュートン法が使われます。しかし、ニュートン法は、適切な初期値を与えないと収束しません。
今回の記事では、初期値に依存しないような方法、ホモトピー法について説明したいと思います。
以下のスライドを参考にしました。
ニュートン法
最初に、ニュートン法について説明する。解きたい方程式を、
\begin{align}
f(x)=0
\end{align}
とする。ここで、既知の値を $x_n$ として、微小変化量 $\delta x$ を使い、
\begin{align}
f(x_n+\delta x)\approx0
\end{align}
となるような微小変化量 $\delta x$を求めることを考える。（ただし $f(x_n)\neq0$。$n$ は、暗に、繰り返し数を意味する。）
この式をテイラー展開することで、
\begin{align}
f(x_n+\delta x)&\approx f(x_n)+\frac{\partial f(x)}{\partial x}\Big|_{x=x_n} \delta x \\
& \approx0
\end{align}
\begin{align}
\therefore \delta x = -\left(\frac{\partial f(x)}{\partial x}\right)^{-1}\Big|_{x=x_n} f(x_n)
\end{align}
微小変化量 $\delta x$が求まる。
ニュートン法のアルゴリズムは、$f(x)=0$ を満たす $x$ を、
\begin{align}
x_{n+1} = x_n -\left(\frac{\partial f(x)}{\partial x}\right)^{-1}\Big|_{x=x_n} f(x_n)
\end{align}
として、繰り返し計算を行い求めていく。１次元の場合、プログラムはこんな感じになる（だろう）。
# 求めたい方程式
def f(x):
return 方程式
# 求めたい方程式の微分
def df_dx(x):
return 方程式の微分
iter = 1
old_x = 初期値
while(True):
print('--------- iter = ' + str(iter) +' ---------')
# δu
delta_x = -f(old_x)/df_dx(old_x)
new_x = old_x + delta_x
r = f(new_x)
iter = iter + 1
print(new_x,r)
# 閾値を決め、誤差が閾値より小さくなったとき終了
if np.abs(r) < 10**(-5) : break
old_x = new_x
ホモトピー法
$f(x)=0$ に対して、初期値 $x_0$ として以下の関数、
\begin{align}
g(x) &= f(x) -f(x_0) \\
h(t,x) &= tf(x)+(1-t)g(x)
\end{align}
を考える。
$t=0$ の場合、$h(0,x)=g(x)$ となり、$t=1$ の場合、$h(1,x)=f(x)$ となる。ホモトピー法は、$t$ を $0$ から $1$ に変化させて$h(t,x)=0$ を解いていき、$t=1$ のとき、求めたい方程式 $f(x)=0$ となる $x$ が求まる。
具体的な $h(t,x)=0$ 解き方として、常微分方程式に帰着して解く方法がある。$h(t,x)$ を全微分すると、
\begin{align}
\mathrm{d} h(t,x) &= \frac{\partial h}{\partial t} \mathrm{d}t + \frac{\partial h}{\partial x} \mathrm{d}x \\
&= 0 \\
\\
\Rightarrow \frac{\mathrm{d} h}{\mathrm{d} t} &= \frac{\partial h}{\partial t} + \frac{\partial h}{\partial x} \frac{\mathrm{d} x}{\mathrm{d} t} \\
&= 0 \\
\end{align}
となるので、
\begin{align}
\frac{\mathrm{d} x}{\mathrm{d} t}= -\left(\frac{\partial h}{\partial x}\right)^{-1}\frac{\partial h}{\partial t}
\end{align}
と求まる。
また
\begin{align}
&\frac{\partial h}{\partial t} = f(x) -g(x) = f(x_0) \\
&\frac{\partial h}{\partial x} = t\frac{\partial f}{\partial x} + (1-t)\frac{\partial f}{\partial x} = \frac{\partial f}{\partial x}
\end{align}
なので、微小量（刻み幅） $\delta$ として $dx/dt$ を単純に差分化すれば
\begin{align}
\frac{\mathrm{d} x}{\mathrm{d} t} &\approx \frac{x_{n+1}-x_n}{\delta} \\
&= -\left(\frac{\partial f}{\partial x}\right)^{-1} f(x_0)
\end{align}
と書ける。
ホモトピー法のアルゴリズムは、$f(x)=0$ を満たす$x$ を、
\begin{align}
x_{n+1} = x_n -\delta \left(\frac{\partial f(x)}{\partial x}\right)^{-1}\Big|_{x=x_n} f(x_0)
\end{align}
として繰り返し計算を行い、求めていく。アルゴリズムは、ニュートン法とほぼ変わらない。１次元の場合、プログラムはこんな感じになる（だろう）。
# 求めたい方程式
def f(x):
return 方程式
# 求めたい方程式の微分
def df_dx(x):
return 方程式の微分
iter = 1
delta = 刻み幅
t = 0.0
x0=初期値
old_x = x0
while(True):
print('--------- iter = ' + str(iter) +' ---------')
# δu
delta_x = -f(x0)/df_dx(old_x)
new_x = old_x + delta*delta_x
r = f(new_x)
t += delta
iter = iter + 1
print(new_x,r)
# t = 1 となったとき終了。小数点に気を付ける。
if t >= 1.0 : break
old_x = new_x
以上のような単純な差分化のホモトピー法では、刻み幅 $\delta$ を小さく取らなければ収束しない欠点がある。
まとめ
ニュートン法は、適切な初期値を与えなければ収束しない欠点がある。ホモトピー法は、そこまで初期値に依存はしないが、単純な差分化のホモトピー法では、刻み幅 $\delta$ を小さく取らなければ収束しない。
ホモトピー法を使う場合、もう少し高度な差分化が必要となる。