ニュートンラフソン法で極値を求める

1変数関数

 f(x) = x^3 -2x^2+x+3極値

# coding:utf-8

def f_dash(x):
    return float(3*x**2 -4*x + 1)

def f_2dash(x):
    return float(6*x - 4)

x = 100 
while True:
    x_nxt = x - f_dash(x)/f_2dash(x)
    delta = abs(x_nxt-x)
    if (delta<0.000001):
        break
    x = x_nxt

print x
1.00000071536

2変数関数

 f(x,y)= x^3+y^3-9xy+27極値

import numpy as np

def nabla_f(p):
    x = p[0].item()
    y = p[1].item()
    J = np.matrix([3*x**2-9*y,3*y**2-9*x]).T 
    return J

def Hesse(p):
    x = p[0].item()
    y = p[1].item() 
    H = np.matrix([[6*x,-9],[-9,6*y]])
    return H

x = np.matrix([3,2]).T
while True:
    df = nabla_f(x)
    H = Hesse(x)
    delta_x = np.linalg.solve(H,-df)
    if (np.linalg.norm(delta_x) < 0.001):
        break
    x = x + delta_x

print x
[[ 3.00022318]
 [ 3.00029142]]