newton.py

# Simple root finding using Newton's method
#
# Function to be optimized: f(x) = x^2 - 4 * sin(x)

import math

def objfun(x):
    '''Compute the value and gradient of our function at x.'''
    f = x**2 - 4 * math.sin(x)
    fp = 2 * x - 4 * math.cos(x)
    return f, fp

tol = 1.0e-14
x0 = -10
f0, fp0 = objfun(x0)
x, f, fp = x0, f0, fp0

while abs(f) > tol * abs(f0):
    x = x - f / fp
    f, fp = objfun(x)
    print x

# Two dimensional problem optimizing the following objective
# f(x1, x2) = (x1^3 - x2 + gamma, -x1 + x2^2 + gamma)

import numpy as np
from numpy.linalg import norm, solve

def objfun2(x, gamma=0.5):
    f = np.array([x[0] ** 3 - x[1] + gamma, -x[0] + x[1] ** 2 + gamma])
    fp = np.array([[3 * x[0] ** 2, -1], 
                   [-1, 2 * x[1]]])
    return f, fp

i = 0
maxiter = 50
tol = 1.0e-14
x0 = np.array([0, 0])
f0, fp0 = objfun2(x0)
x, f, fp = x0, f0, fp0

while norm(f) > tol * norm(f0) and i <= maxiter:
    s = solve(fp, -f)
    x = s + x
    f, fp = objfun2(x)
    i += 1
    print x
    print f