Pythonでニュートン法を実装しました。

この本のニュートン法を実装していきます。
www.amazon.co.jp

1変数の場合

第3章のアルゴリズム3.3(1変数のニュートン法)を、Pythonで実装してみました。

import numpy as py

def f(x):
	return - x ** 3
def f_p(x):
	return - 3 * x ** 2

def f_pp(x):
	return - 6 * x

def newton(x,fp,fpp):
	iteration = 100
	delta = 0.00001
	for i in range(iteration):
		x_new = x- ( fp(x)/fpp(x) )
		if abs(x_new -x) < delta:
			break
		x = x_new
	return(x)

print(newton(-5.,f_p,f_pp))

多変数の場合

次に、多変数の場合で実装しました。下記は2変数の場合に対応しています。教科書ではアルゴリズム3.4(多変数のニュートン法)です。関数は例題3.2に合わせました。

import numpy as np

def f(x):
	return x[0] ** 3 + x[1] ** 3 - 9 * x[0] * x[1] + 27

def df(x):
	return np.array([3 * x[0] ** 2 - 9 * x[1], 3 * x[1] ** 2 - 9 * x[0]])

def H(x):
	return np.array([[6*x[0],-9],[-9,6*x[1]]])

def newton_multi(x,H):
	iteration = 10000
	delta = 0.00001
	for i in range(iteration):
		x_new = x - np.linalg.inv(H(x)).dot(df(x))
		if np.linalg.norm(x - x_new) < delta:
			break
		x = x_new
	return(x)

print(newton_multi(np.array([-5,-5]),H))