Newton-Raphson method using Python SymPy

Newton-Raphson method

Newton-Raphson is a very popular method for the numerical calculation of an equation’s root. Its’ basic concepts for formulation originate from the Taylor theorem and of course the fact that function value becomes zero at the root point. You can read more about the method in the Wikipedia entry.

SymPy

SymPy is a Python library which enables the use of symbolic mathematical operations. It can solve equations, differentiate or integrate, simplify complex expressions or evaluate mathematical functions.

In various distributions you can install the package directly with the help of the corresponding package manager. For example in Ubuntu:

apt-get install python-sympy

 

 

For the latest version, you can use pip:


pip install sympy

Code

Our “solver” will receive the precision as an optional parameter, with default value of 5⋅10-6, the function representation in a SymPy compatible format and letting x be the symbolic variable (check the docs) and the starting guess value.

#!/usr/bin/env python
import argparse
import math
import sys
import time
import sympy

def get_parameters():
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "-f", "--function",
        help="Define a function"
    )
    parser.add_argument("-s", "--starting", 
                        help="Starting point value", type = float, default = 0.0)
    parser.add_argument("-p", 
                        "--precision", 
                        help="Convergence precision", 
                        type = float, default = 5*10**(-6))
    return parser.parse_args()

if __name__ == "__main__":
    parameters = get_parameters()

    sym_x = sympy.Symbol('x')
    # convert the given function to a symbolic expression
    try:
        fx = sympy.S(args.function)
    except sympy.core.sympify.SympifyError:
        sys.exit('Unable to convert function to symbolic expression.')

    # calculate the differential of the function
    dfdx = diff(fx, Symbol('x'))

    # e is the relative error between 2 consecutive 
    # estimations of the root
    e = 1
    x0 = parameters.starting
    iterations = 0
start_time = time.time()
while e > args.precision:
        # new root estimation
        try:
            r = x0 - fx.subs({sym_x: x0}) / dfdx.subs({sym_x: x0})
        except ZeroDivisionError:
            print "Function derivative is zero. Division by zero, program will terminate."
            raise
        # relative error
        e = abs((r - x0) / r)
        iterations += 1
        x0 = r

total_time = time.time() - start_time

print 'Function:'
sympy.pprint(fx)
print 'Derivative:'
sympy.pprint(dfdx)
print 'Root %10.6f calculated after %d iterations' % (r, iterations)
print 'Function value at root %10.6f' % fx.subs({sym_x : r})
print 'Finished in %10.6f seconds'% total_time

Example

We want to solve the equation

6x2 = 400

near the initial value 910 (the positive root).
This is the output of a test run.

$ ./nr.py -f '6*x**2 - 400' -s 910 

Function:

6⋅x - 400

Derivative:

12⋅x

Root 8.164966 calculated after 11 iterations

Function value at root 0.000000

Finished in 0.035108 seconds

 

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.