Python ニュートン法を実装する

# ニュートン法とは

1変数では初期値 $x_0$ を設定し、以下の数列によって関数 $f$ の近似解を求める方法。

```math
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} 
```

多変数では初期値 $\boldsymbol{x_0} \in \mathbb{R}^n$ を設定し 、関数 $f: \mathbb{R}^n \rightarrow \mathbb{R}^n$ の近似解を求める方法。

$$\boldsymbol{x_{n + 1}} = \boldsymbol{x_n} - H^{-1} f(\boldsymbol{x_n})$$

ここで $H$ はヘッセ行列。

# 1変数

こちらに置きました。

https://github.com/quryu-1/newton

```python
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt

def newton_method(f, x, init_value, num):
    df = sp.diff(f, x)
    next_value = 0
    for i in range(num):
        next_value = init_value - f.subs(x, init_value) / df.subs(x, init_value)
        init_value = next_value
    return next_value
    
if __name__ == '__main__':
    x = sp.Symbol('x')
    # ここで関数を設定する
    f = sp.sin(x) - 0.5
    # ここで初期値を設定する
    init_value = 0.1
    # ここで繰り返し回数を設定する
    num = 10
    solution = newton_method(f, x, init_value, num)
    print(solution)
```

## 例 1

$${\rm sin}(x) - 0.5 = 0$$

初期値

```
0.5
```

出力

```
0.523598775598299
```

計算機での $\pi / 6$ の計算結果

```
0.52359877559829887307710723054658381403
```

誤差

```
1.2692289276945341618597101189910806429 E-16
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/ddb5231f-a7fb-b300-9d37-ebd1ce4da1fb.png)

## 例 2

$$e^x - 2 = 0$$

初期値

```
1.0
```

出力

```
0.693147180559945
```

計算機での $e^x - 2 = 0$ の計算結果

```
0.69314718055994530941723212145817656807
```

誤差

```
3.0941723212145817656806814692747176370 E-16
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/679deec8-0085-a960-ec83-1438f4466794.png)

# 多変数

```python
import sympy as sp
import numpy as np

def newton_method(f, variables, init_value, num):
    for i in range(num):
        hesse = get_hesse(f, variables, init_value)
        subs_vector = np.array(
            [[f[i].subs(init_value)] for i in range(len(f))],
            dtype=float
        )
        init_vector = np.array(
            [[v] for v in list(init_value.values())],
            dtype=float
        )
        
        next_value = init_vector - np.linalg.inv(hesse) @ subs_vector
        init_value = {key: next_value[n][0] for n, key in enumerate(init_value)}
    return init_value

def get_hesse(f, variables, init_value):
    hesse = []
    for i in range(len(variables)):
        l = []
        row = []
        l.append(variables[i])
        for j in range(len(variables)):
            if len(l) == 1:
                l.append(variables[j])
            else:
                l[1] = variables[j]
            row.append(
                sp.diff(f[i], *l).subs(
                    init_value
                )
            )
        hesse.append(row)
    return np.array(hesse, dtype=float)

if __name__ == '__main__':
    x = sp.Symbol('x')
    y = sp.Symbol('y')
    z = sp.Symbol('z')
    variables = [x, y, z]
    # ここで関数を設定する
    f = [x ** 2 - 2, (x ** 2) * (y ** 2) - 4, x * z ** 2 - 1]
    # ここで初期値を設定する
    init_value = {x: 1, y: 2, z: 0.3}
    # ここで繰り返し回数を設定する
    solution = newton_method(f, variables, init_value, 50)
    print(solution)
```

## 例 1

```math
\left\{
\begin{array}{l}
x^2 - 2 = 0 \\
x^2 y^2 - 4 = 0 \\
x z ^2 - 1 = 0
\end{array}
\right.
```

初期値

```
{x: 1, y: 2, z: 0.3}
```

出力

```
{x: 1.414213562373095, y: 1.4142135623730954, z: 0.8408964152537146}
```

解（初期値に近い値）

$$x = \sqrt{2}, y = \sqrt{2}, z = \sqrt{1/\sqrt{2}}$$

上記の解の計算機での計算結果

```
x = 1.4142135623730950488016887242096980786
y = 1.4142135623730950488016887242096980786
z = 0.84089641525371454303112547623321489504
```

誤差

```
x: 4.8801688724209698078596369014359646933 E-17
y: 3.5119831127579030192140409244737913280 E-16
z: 5.6968874523766785104960243772098740359 E-17
```

# 参考記事

https://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95

https://science-log.com/%E6%95%B0%E5%AD%A6/%E3%80%90%E6%9C%80%E9%81%A9%E5%8C%96%E5%95%8F%E9%A1%8C%E3%81%AE%E5%9F%BA%E7%A4%8E%E3%80%91%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95%E3%81%A8%E3%83%98%E3%83%83%E3%82%BB%E8%A1%8C%E5%88%97/