Week 4

Examples for finding roots of equations

In [6]:
# illustrates one possible implementation of the Newton-Raphson algorithm
#
# fabs() is the floating-point ABSOLUTE VALUE function
#


from math import *

def func1(x):
    """
    This function is for the eqn:  x - exp(-x) = 0
    """
    df = 1.0 + exp(-x)         # change these lines
    return x-exp(-x), df       # if you want to solve a different equation

#--------------------------------------------------------------------------

def Newton(func, xstart, xend, x):
    """
    Finds an approximation for the real root x of the function 'func' bracketed 
    in the interval [xstart, xend] using the Newton-Raphson algorithm
    input:    func, initial approximation
    output:   the approximation of the root
    """
    
    EPS = 1e-10       # rel precision of the approximation
    ITMAX = 100       # max no of iterations
    
    for i in range(1, ITMAX+1):
        f, df = func(x)              # get function & its derivative at x
        if fabs(df) > EPS:           # the actual Newton-Raphson iteration
            dx = -f/df
        else:
            -f                       # the "tricky bit"
        x += dx
        if  ((x < xstart) or (x > xend)):                 # no root in [xstart, xend] 
            print('Interval does not contain a root')
            return x
        if fabs(dx) <= EPS*fabs(x):                       # check the convergence
            print('Everything\'s fine')
            return x
    
    print('Max no. iterations exceeded')                  # if it gets here, 
    return x                                              # too many iterations....
    
#-----------------------------------------------------------------                  
                  
def main():
      # here we test the above function:
      a, b = -1.0, 1.0        # interval in which the root is bracketed
      x = 4.0                 # initial guess
    
      # the above values are equation-specific...
      # change them accordingly for other equations
    
      x = Newton(func1, a, b, x)      # get approximation
    
      print('The solution is: ', format(x, '.6f'))     # display result on the screen
                  
main()
Everything's fine
The solution is:  0.567143
In [7]:
# Another (basic) implementation of Newton-Raphson:

# Constants used in the program:
TOL = 1.0E-7
ITMAX = 100

def f(x): return (x-1.0)*(x-5.0)         # we are solving f(x) = 0

def dfdx(x): return 2.0*x-6.0

def solve(f, x):
    n=0                                      # keeps track of no iterations
    # The actual Newton-Raphson iteration
    while abs(f(x)) > TOL and n < ITMAX:  
        x = x - f(x)/dfdx(x)
        n += 1
    return x, n, f(x)            # note that we return a tuple
    
def main():
    xguess = 0.3                 # initial guess
    res = solve(f,xguess)        # catch the answer in 'res' (tuple)
    
    # display info/results on the screen:
    print('The root is:', res[0])      
    print('Number of iterations:', res[1])
       
main()
The root is: 0.9999999999997655
Number of iterations: 4
In [10]:
# Derivative approximation for Newton-Raphson:
# This is the same as above, but instead of the exact formula
# for the derivative, we use the definition (of the derivative) to 
# find an approximation

# constants used in the program:
TOL = 1.0E-7
ITMAX = 100
h = 1.0E-8                               # needed for dfdx


def f(x): return (x-1.0)*(x-5.0)         # we are solving f(x) = 0
                                         # change the returned expression 
                                         #             for other equations
#----------------------------------------------------------------------------
def dfdx(f, x):
    # remember the definition of the derivative:
    return (f(x+h)-f(x))/float(h)
#----------------------------------------------------------------------------
def solve(f, dfdx, x):
    n=0                                      # keeps track of no iterations
    # The actual Newton-Raphson iteration
    while abs(f(x)) > TOL and n < ITMAX:  
        x = x - f(x)/dfdx(f,x)
        n += 1
    return x, n, f(x)            # note that we return a tuple
    
def main():
    xguess = 16.3                    # (ridiculous) initial guess
    res = solve(f, dfdx, xguess)     # "catch" the answer in 'res' (tuple)
    
    # display info/results on the screen:
    print('The root is:', res[0])      
    print('Number of iterations:', res[1])
       
main()
The root is: 5.000000015094005
Number of iterations: 6