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

Question

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?

Answer 1

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)

Answer 2

A common way of approximating the derivative of f at x is using a finite difference:

f'(x) = (f(x+h) - f(x))/h                   Forward difference
f'(x) = (f(x+h) - f(x-h))/2h                Symmetric

The best choice of h depends on x and f: mathematically the difference approaches the derivative as h tends to 0, but the method suffers from loss of accuracy due to catastrophic cancellation if h is too small. Also x+h should be distinct from x. Something like h = x*1e-15 might be appropriate for your application. See also implementing the derivative in C/C++.

You can avoid approximating f' by using the secant method. It doesn't converge as fast as Newton's, but it's computationally cheaper and you avoid the problem of having to calculate the derivative.

Answer 3

You can approximate fprime any number of ways. One of the simplest would be something like:

lambda fprime x,dx=0.1: (f(x+dx) - f(x-dx))/(2*dx)

the idea here is to sample f around the point x. The sampling region (determined by dx) should be small enough that the variation in f over that region is approximately linear. The algorithm that I've used is known as the midpoint method. You could get more accurate by using higher order polynomial fits for most functions, but that would be more expensive to calculate.

Of course, you'll always be more accurate and efficient if you know the analytical derivative.