Newton.py

def Newton(f, x, dfdx, epsilon=1.0E-7, N=100, store=False):
    f_value = f(x)
    n = 0
    if store: info = [(x, f_value)]
    while abs(f_value) > epsilon and n <= N:
        dfdx_value = float(dfdx(x))
        if abs(dfdx_value) < 1E-14:
            raise ValueError("Newton: f'(%g)=%g" % (x, dfdx_value))

        x = x - f_value/dfdx_value

        n += 1
        f_value = f(x)
        if store: info.append((x, f_value))
    if store:
        return x, info
    else:
        return x, n, f_value


def _g(x):
    return exp(-0.1*x**2)*sin(pi/2*x)

def _dg(x):
    return -2*0.1*x*exp(-0.1*x**2)*sin(pi/2*x) + \
           pi/2*exp(-0.1*x**2)*cos(pi/2*x)

def _test():
    from scitools.std import sin, cos, exp, linspace, plot, pi
    import sys

    x0 = float(sys.argv[1])
    x, info = Newton(_g, x0, _dg, store=True)
    print 'root: %.16g' % x
    for i in range(len(info)):
        print 'Iteration %2d: f(%g)=%g' % \
              (i, info[i][0], info[i][1])

    x = linspace(-7, 7, 401)
    y = _g(x)
    plot(x, y)

if __name__ == '__main__':
    _test()