Search code examples
pythonmathpython-3.xsolvernewtons-method

Given f, is there an automatic way to calculate fprime for Newton's method?


The following was ported from the pseudo-code from the Wikipedia article on Newton's method:

#! /usr/bin/env python3
# https://en.wikipedia.org/wiki/Newton's_method
import sys

x0 = 1
f = lambda x: x ** 2 - 2
fprime = lambda x: 2 * x
tolerance = 1e-10
epsilon = sys.float_info.epsilon
maxIterations = 20

for i in range(maxIterations):
    denominator = fprime(x0)
    if abs(denominator) < epsilon:
        print('WARNING: Denominator is too small')
        break
    newtonX = x0 - f(x0) / denominator
    if abs(newtonX - x0) < tolerance:
        print('The root is', newtonX)
        break
    x0 = newtonX
else:
    print('WARNING: Not able to find solution within the desired tolerance of', tolerance)
    print('The last computed approximate root was', newtonX)

Question

Is there an automated way to calculate some form of fprime given some form of f in Python 3.x?


Solution

  • Answer

    Define the functions formula and derivative as the following directly after your import.

    def formula(*array):
        calculate = lambda x: sum(c * x ** p for p, c in enumerate(array))
        calculate.coefficients = array
        return calculate
    
    def derivative(function):
        return (p * c for p, c in enumerate(function.coefficients[1:], 1))
    

    Redefine f using formula by plugging in the function's coefficients in order of increasing power.

    f = formula(-2, 0, 1)
    

    Redefine fprime so that it is automatically created using functions derivative and formula.

    fprime = formula(*derivative(f))
    

    That should solve your requirement to automatically calculate fprime from f in Python 3.x.

    Summary

    This is the final solution that produces the original answer while automatically calculating fprime.

    #! /usr/bin/env python3
    # https://en.wikipedia.org/wiki/Newton's_method
    import sys
    
    def formula(*array):
        calculate = lambda x: sum(c * x ** p for p, c in enumerate(array))
        calculate.coefficients = array
        return calculate
    
    def derivative(function):
        return (p * c for p, c in enumerate(function.coefficients[1:], 1))
    
    x0 = 1
    f = formula(-2, 0, 1)
    fprime = formula(*derivative(f))
    tolerance = 1e-10
    epsilon = sys.float_info.epsilon
    maxIterations = 20
    
    for i in range(maxIterations):
        denominator = fprime(x0)
        if abs(denominator) < epsilon:
            print('WARNING: Denominator is too small')
            break
        newtonX = x0 - f(x0) / denominator
        if abs(newtonX - x0) < tolerance:
            print('The root is', newtonX)
            break
        x0 = newtonX
    else:
        print('WARNING: Not able to find solution within the desired tolerance of', tolerance)
        print('The last computed approximate root was', newtonX)