Python ニュートン法を実装する
ニュートン法とは
1変数では初期値 $x_0$ を設定し、以下の数列によって関数 $f$ の近似解を求める方法。
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変数
こちらに置きました。
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
例 2
$$e^x - 2 = 0$$
初期値
1.0
出力
0.693147180559945
計算機での $e^x - 2 = 0$ の計算結果
0.69314718055994530941723212145817656807
誤差
3.0941723212145817656806814692747176370 E-16
多変数
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
\left\{
\begin{array}{l}
x^2 - 2 = 0 \\
x^2 y^2 - 4 = 0 \\
x z ^2 - 1 = 0
\end{array}
\right.
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)
\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
参考記事
Author And Source
この問題について(Python ニュートン法を実装する), 我々は、より多くの情報をここで見つけました https://qiita.com/quryu/items/08035216f6ebabed9353著者帰属:元の著者の情報は、元のURLに含まれています。著作権は原作者に属する。
Content is automatically searched and collected through network algorithms . If there is a violation . Please contact us . We will adjust (correct author information ,or delete content ) as soon as possible .