Sunday, March 8, 2015

py34. Newton's method in Python

Newton's method is based on the approximating a function with only its first Taylor Series term:


y(x+dx) = y(x) + dx y'(x)


We are interested in dx, which will result in y(x + dx) = 0.


0 = y(x) + dx y'(x).


-y(x) = dx y'(x).


dx = -y(x)/y'(x).


A function newton is written implementing this, and used to calculate cube root of 3. Later, we will use it for more complicated cases.

# ex34.py
from __future__ import print_function, division

def fit(x):
    return 3-x**3

def newton(f, tol=1E-10, x0=1, N=50, h=0.2):
    for i in range(N):
        df_dx = (f(x0+h/2)-f(x0-h/2))/h
        x1 = x0 - f(x0)/df_dx
        error = abs(x1-x0)
        print("i = %d,x1 = %.4f, error = %.4e" % (i,x1,error))
        if error < tol: break
        x0 = x1
    return x1

cube_root_3 = newton(fit)
print('cube_root_3 =',cube_root_3,
      '\nfit(cube_root_3) = ',fit(cube_root_3))
#i = 0,x1 = 1.6645, error = 6.6445e-01
#i = 1,x1 = 1.4708, error = 1.9363e-01
#i = 2,x1 = 1.4428, error = 2.7982e-02
#i = 3,x1 = 1.4423, error = 5.9354e-04
#i = 4,x1 = 1.4422, error = 1.1936e-06
#i = 5,x1 = 1.4422, error = 1.9106e-09
#i = 6,x1 = 1.4422, error = 3.0567e-12
#cube_root_3 = 1.44224957031 
#fit(cube_root_3) =  -3.15303338994e-14

No comments:

Post a Comment