ホモトピー法

# 概要
　有限要素法(FEM)や粒子法で非線形方程式を解く際、ニュートン法が使われます。しかし、ニュートン法は、適切な初期値を与えないと収束しません。

　今回の記事では、初期値に依存しないような方法、ホモトピー法について説明したいと思います。

以下のスライドを参考にしました。

https://www.slideshare.net/HideoHirose/homotopy-67988638

# ニュートン法

　最初に、ニュートン法について説明する。解きたい方程式を、

```math
\begin{align}
f(x)=0
\end{align}
```
とする。ここで、既知の値を $x_n$ として、微小変化量 $\delta x$ を使い、

```math
\begin{align}
f(x_n+\delta x)\approx0
\end{align}
```
となるような微小変化量 $\delta x$を求めることを考える。（ただし $f(x_n)\neq0$。$n$ は、暗に、繰り返し数を意味する。）

　この式をテイラー展開することで、

```math
\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}
```

```math
\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$ を、

```math
\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}
```
として、繰り返し計算を行い求めていく。１次元の場合、プログラムはこんな感じになる（だろう）。


```php
#　求めたい方程式
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$ として以下の関数、

```math
\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)$ を全微分すると、 

```math
\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}
```
となるので、

```math
\begin{align}
\frac{\mathrm{d} x}{\mathrm{d} t}= -\left(\frac{\partial h}{\partial x}\right)^{-1}\frac{\partial h}{\partial t} 
\end{align}
```
と求まる。
　また

```math
\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$ を単純に差分化すれば

```math
\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$ を、

```math
\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}
```
として繰り返し計算を行い、求めていく。アルゴリズムは、ニュートン法とほぼ変わらない。１次元の場合、プログラムはこんな感じになる（だろう）。

```php
#　求めたい方程式
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$ を小さく取らなければ収束しない。
　
　ホモトピー法を使う場合、もう少し高度な差分化が必要となる。