-Python- ニュートン法

方程式

 x^3 - 5x + 1 = 0

をNewton法で解いてみるスクリプトを考えます.


# Definition of Function
def newton1dim(f, df, x0, eps = 1e-10, max_iter = 1000): 
    x = x0
    iter = 0
    while True: 
        x_new = x - f(x)/df(x)
        if abs(x - x_new) < eps: 
            break
        x = x_new
        iter += 1
        if iter == max_iter
            break
    return x_new
 
def f(x): 
    return x ** 3 - 5 * x + 1
 
def df(x): 
    return 3 * x**2 - 5
 
# Display Results
print(newton1dim(f, df, 2))
print(newton1dim(f, df, 0))

print(newton1dim(f, df, -3))
 
このスクリプトでは,newton1dimで解を計算しています.引数のfは解を求めたい関数で,dfはその導関数です.x0は初期値で,epsは収束条件に用いるεです.デフォルト値付き引数として max_iter を指定していますが,これは初期値によっては無限ループとなる可能性があるので,繰り返しの最大回数を指定しています.最大繰り返し回数を超えた場合には,そのままループから抜け出します.

実行結果は以下のようになります.

2.1284190638445777
0.20163967572340463
-2.330058739567982
 
ニュートン法は多次元の場合に拡張することができます.

 

# Import Module

import numpy as np

from numpy import linalg

 

# Definition of Class

class Newton: 

    def __init__(self, f, df, eps = 1e-10, max_iter = 1000): 

        self.f = f

        self.df = df

        self.eps = eps

        self.max_iter = max_iter

    

    def solve(self, x0): 

        x = x0

        iter = 0

        self.path_ = x0.reshape(1, -1)

        while True: 

            x_new = x - np.dot(linalg.inv(self.df(x)), self.f(x))

            self.path_ = np.r_[self.path_, x_new.reshape(1, -1)]

            if ((x - x_new) ** 2).sum() < self.eps*self.eps: 

                break

                x = x_new

                iter += 1

                if iter == self.max_iter: 

                    break

            return x_new

 

ここではクラスとして実装し,xの値の軌跡を保存するようにしています.

18行目では,更新式の計算を行なっています.逆行列の計算には関数 np.linalg.invを使っています.


# Import Module

import numpy as np

import matplotlib.pyplot as plt

import newton

 

# Definition of Function

def f1(x, y): 

    return x**3 - 2 * y

 

def f2(x, y): 

    return x**2 + y**2 - 1

 

def f(xx): 

    x = xx[0]

    y = xx[1]

    return np.array([f1(x, y), f2(x, y)])

 

def df(xx): 

    x = xx[0]

    y = xx[1]

    return np.array([[3 * x**2, -2], [2 * x, 2 * y]])

 

# Drawing Graph

xmin, xmax, ymin, ymax = -3, 3, -3, 3

plt.xlim(xmin, xmax)

plt.ylim(ymin, ymax)

x = np.linspace(xmin, xmax, 200)

y = np.linspace(ymin, ymax, 200)

xmesh, ymesh = np.meshgrid(x, y)

z1 = f1(xmesh, ymesh)

z2 = f2(xmesh, ymesh)

plt.contour(xmesh, ymesh, z1, colors = "r", levels = [0])

plt.contour(xmesh, ymesh, z2, colors = "k", levels = [0])

 

# Solving problem

solver = newton.Newton(f, df)

 

initials = [np.array([1, 1]), 

            np.array([-1, -1]), 

            np.array([1, -1])]

markers = ["+", "*", "x"]

 

for x0, m in zip(initials, markers): 

    sol = solver.solve(x0)

    plt.scatter(solver.path_[:, 0], 

                solver.path_[:, 1], color = "k", marker = m)

    print(sol)

 

plt.show()

この例では初期値(1, -1),(-1, -1),(1, -1)を使って方程式を解いています.その時のxkの動きを可視化しています.f1(x, y) = 0とf2(x, y) = 0を表す曲線は等高線の仕組みを使って描画しています.


実行結果は以下のようになります.

[1.  0.5]
[-1.  -0.5]
[-1.  -2.5]
 

f:id:HidehikoMURAO:20181106200617p:plain

実行結果のグラフを見ると,(1, 1),(-1, -1)と比べて初期値(1, -1)から始めた場合(図中のX印)は,少し遠回りしてから収束に向かっていることがわかります.