奥が深いニュートン法について解説してみた

ニュートン法とは

ニュートン法はある方程式の数値解を近似的に求める手法です。例えば下のような関数の数値解を求めたいとき、解は√2と求まりますが、小数での表現をしたい場合、数値解析の反復法を用いて計算することができます。

f(x) = x^2 - 2  

ニュートン法でつまづきやすい部分は初期値の設定です。これをミスると何回もループしても計算し切れないといった問題が発生します。今回取り扱う関数は2次関数で平易ですので、あまり深く考えずに実装しましょう。
具体的に下のような非線形方程式は初期値によって収束せずに解が振動するようです。
初期値(x0=3)とすると

g(x) = 3\arctan (x - 1) + \frac{x}{4}

このように接線が描ける。

求めたいg(x)=0の解からどんどん離れていくことが視覚的にわかると思います

このように解が振動していくため、最大のループ回数を設定する必要があリます。上の関数の場合は初期値を2.5に設定すれば収束するそうです！

ニュートン法のアルゴリズム

下の関数の数値解を近似的に求める。(√2を小数で表す) 初期値x0=3と設定すると

f(x) = x^2 - 2  

x0=3における接線は以下のように表せる。

y - f(x_0) = f'(x_0)(x-x_0)  

f'(x_0) = 2x_0  \quad f(x_0) = 7

y = 6(x-3)+7  \quad y = 6x - 11

ここで接線とy軸が交わる点を次のステップのx1とする。

x_1 = \frac{11}{6}

y - f(x_1) = f'(x_1)(x-x_1)  

y - f(\frac{11}{6}) = f'(\frac{11}{6})(x-\frac{11}{6})  

同様に接線とy軸が交わる点を次のステップのx2とする。
このような操作を繰り返してf(x)=0の解が求まります。
この操作を一般化し、漸化式で表すと

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

具体的な値を計算で求めたい方は以下のPythonコードを参考にしてください。最大のループ回数は100回、打ち切り誤差は1e-7に設定しました。このコードをいじって√3の値や2の1/10乗の値を求めてみてください。

ニュートン法

def f(x):
    return x**2.0-2

def newton(x0,eps=1e-7,error=1e-7,max_loop=100,cnt=0):
    while True:
        df_x=(f(x0+eps)-f(x0-eps))/(2*eps)
        if df_x <= error:
            print("slope is too small")
            exit()
        
        x1=x0-f(x0)/df_x
        cnt+=1
        if abs(x1-x0)<=error or cnt>=max_loop:
            break
        x0=x1
        print(x0)
    return x0
print(3.0)
print(newton(3.0))

実行結果

このような単純な2次関数だとループ回数は少なくて済みますね！