大学で物理学を専攻しています。今回（これが初投稿だけど・・・）はニュートン法と二分法を Python を使って解いてみようと思います。
ニュートン法とは
ニュートン法、またはニュートン・ラフソン法は、数値解析の分野いおいて、方程式を数値計算によって解くための反復法による球根アルゴリズムの 1 つです。
手順としては、関数$f(x)$について、
- ある$x_1$の数値を決定する
- $f(x_1)$を求め、その$x_1$の接戦の式より$y=0$と交わる$x$の値を$x_2$とする
- $f(x_2)$を求め、その$x_2$の接線の式より$y=0$と交わる$x$の値を$x_3$とする
- これを繰り返すことによって、$x_1, x_2$は$f(x)=0$となる$x=a$の値に近づいていく
といった流れになっています。
この一連の手順は、定式化されています。
接線の式は$f(x_n)=f'(x_n)(x_n - x_{n+1})$と表すことができるので、式変形すると、
$$
f(x_n)=f'(x_n)(x_n - x_{n+1})
\Leftrightarrow x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}
$$
といった形になります。
Pythonを用いて表現してみよう
具体例として、**$f(x)=e^x-4x$ に対して $f(x)=0$**を解いていきます。
微分については$f'(x)=e^x -4$であることを用いります。$f(x), f'(x)$をそれぞれ f(x), df(x)と定義し、初期値を$x=4$とします。
実際のコード
import numpy as np
import matplotlib.pylab as plt
def f(x):
return np.exp(x)-4*x
def df(x):
return np.exp(x)-4
#初期値
x=4
#epsilon1は要求精度。これよりも数値精度がよくなったら二分法の繰り返しを終了する。
epsilon1=0.00001
# n, c, f(c) を値を出力する。\t はタブという空白を適当にあける命令です。
print("n\t x\t f(x)")
n=1
# while文でf(x)=0になるまで繰り返す
while True:
x = x - f(x)/df(x)
print("{}\t{:.5f}\t{:.5f}".format(n, x, f(x)))
n += 1
if abs(f(x)) < epsilon1:
break
print("x= %f" % x)
すると結果は$x= 2.153292$という値になります。
確認してみよう
実際に Python を用いてグラフを作成して、確認してみたいと思います。
使用するコードは以下の通りです。
import numpy as np
import matplotlib.pylab as plt
def f(x):
return np.exp(x)-4*x
#0<x<4の間で描画
x=np.arange(0.0,4.0,0.05)
#zは、f(x)=0を描画
z=np.zeros(len(x))
plt.plot(x,f(x),x,ｚ)
これより、グラフは次のように描かれます。
このグラフから、実際に、$x= 2.153292$に答えがあるとわかります。
余談
ニュートン法は、scipy.optimize に入っているので
import scipy.optimize as spopt
def f(x):
return np.exp(x)-4*x
spopt.newton(f,4)
とすれば、簡単に求めることもできます。scipy.optimize についてはこちらよりご参照ください
ご高覧ありがとうございました。