Optimisation with Newton Raphson#

Root finding scripts can be used just by changing the testing method, the newton raphson method for root finding remains unchanged except that it is used on the first derivative instead of the original function which means that the second derivative replaces the first derivative in root finding. If the first derivative can be found it generally is better to use the newton raphson method straight away.

Using the following script, which uses both the first and second derivatives, the following results:

+------+-----------+-------------+-------------+------------+
| step |     x0    |      df     |     d2f     |   error    |
+------+-----------+-------------+-------------+------------+
|  0   |     1     | -87.4336294 | 212.5663706 | -0.4113239 |
|  1   | 1.4113239 | -32.4697385 |  83.7122749 | -0.3878731 |
|  2   |  1.799197 |  -8.2823778 |  46.9058612 | -0.1765745 |
|  3   | 1.9757715 |  -0.788624  |  38.4974066 | -0.0204851 |
|  4   | 1.9962566 |  -0.0081486 |  37.7072757 | -0.0002161 |
|  5   | 1.9964727 |    -9e-07   |  37.6991127 |    -0.0    |
+------+-----------+-------------+-------------+------------+
The approximate value of x is: 1.9964726889295161
Made in 5 iterations

The answer was made in 5 iterations, and only required one starting point, in this case 1. If the starting point was 5 it required 6 iterations:

+------+-----------+--------------+-------------+------------+
| step |     x0    |      df      |     d2f     |   error    |
+------+-----------+--------------+-------------+------------+
|  0   |     5     |  58.8318531  |  14.1663706 | 4.1529235  |
|  1   | 0.8470765 | -128.7206394 | 341.6164924 | -0.3767987 |
|  2   | 1.2238752 | -51.3817768  | 121.6648252 | -0.4223224 |
|  3   | 1.6461976 | -16.2140978  |  57.3979599 | -0.2824856 |
|  4   | 1.9286832 |  -2.6464825  |  40.4434532 | -0.0654366 |
|  5   | 1.9941198 |  -0.0888081  |  37.7881818 | -0.0023502 |
|  6   | 1.9964699 |  -0.0001046  |  37.6992167 |  -2.8e-06  |
+------+-----------+--------------+-------------+------------+
The approximate value of x is: 1.9964699371211263
Made in 6 iterations

Show/Hide Code newton_optimize.py

'''
uses The central problem of optimization is minimization of functions. Let us
first consider the case of univariate functions, i.e., functions of a single real
variable.

xk+1 = xk - f´(xk)/f´´(xk)

using first and second derivatives
'''

# find minimum area given volume of a cylinder

'''
V = π r² h
A = 2 π r h + 2 π r²
A = 2 π r (h + r) = 2 π r (r + V / π r²) = 2 π r² + 2 V / r
'''
from math import pi
from prettytable import PrettyTable

f = lambda r: 2 * (pi*r**2 + 50/r)
df = lambda r: 4*pi*r -100/r**2
d2f = lambda r: 4*pi + 200/r**3
t = PrettyTable(['step', 'x0', 'df', 'd2f', 'error'])

def newtonOpt(df, d2f, x0, tol=1e-5, stepsmax=20):
    df0 = df(x0)
    d2f0 = d2f(x0)
    h = df0/d2f0
    steps=0
    t.add_row([steps, round(x0,7), round(df0,7), round(d2f0,7), round(h,7)])
    while abs(h) >= tol:
        x0 = x0 - h
        df0 = df(x0)
        d2f0 = d2f(x0)
        h = df0/d2f0

        steps += 1
        t.add_row([steps, round(x0,7), round(df0,7), round(d2f0,7), round(h,7)])
        if steps > stepsmax:
            return x0, steps

    return x0, steps

fn, steps = newtonOpt(df, d2f, 1)

print(t)
print("The approximate value of x is: "+str(fn))
print('Made in '+str(steps)+' iterations')

Applying a Newton-Raphson method to a simple complex function, one obtains a fractal:-

../_images/newton_fractal.png

Fractal using \(z = x + i\cdot y\)#

This Page