Secant and Newton-Raphson method - Executable Engineering: Interactive Numerical Methods in Python

import numpy as np
import plotly.graph_objects as go

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

from executable_engineering import RootFinderOpen

Secant method¶

Let’s reconsider the method of False Position but now we always disregard bracket updating and simply take $c$ as our new guess. Our algorithm is now:

Take 2 initial guesses: $x^i$
Calculate the next guess:

x^{i+1} = {x^i} - f(x^{i}) \frac{x^i-x^{i-1}}{f(x^i)-f(x^{i-1})}

(1)

Check if tolerance is met

(Which tolerance?)

def f(x):
  return x**2 - 2

exe.plot_secant(f, 1, 2, tolerance=1e-6, max_iterations=100)

The Newton-Raphson method¶

The Secant method maintains superlinear convergence... but we can do better with a little more information...

Take another look at the fraction in the Secant method update equation:

x^{i+1} = {x^i} - f(x^{i}) \frac{x^i-x^{i-1}}{f(x^i)-f(x^{i-1})}

(2)

This is an (inverse) approximation of $\frac{\partial f}{\partial x}$ , the derivative of $f$ ! Typically, we are able to find this quantity, and the algorithm becomes:

x^{i+1} = {x^i} - \frac{f(x^{i})}{f'(x^i)}

(3)

or, solving for the increment $\Delta x = x^{i+1}-x^i$ and dropping the indicies,

\begin{align} \Delta x &= - \frac{f(x)}{f'(x)} \\ f'(x) \Delta x &= - f(x) \end{align}

(4)

def f(x):
  return x**2 - 2

def df(x):
  return 2*x

exe.plot_newton_raphson(f, df, 1, tolerance=1e-6, max_iterations=100)

Interactive Comparison¶

The Newton-Raphson method has quadratic convergence ( $k=2$ ) near the root which is a great result! It does so, however at the cost of calculating the Jacobian and solving a linear system.

f = lambda x: x**3 - 27
fprime = lambda x: 3 * x**2

x0 = 27.0
x1 = 25.0
tol = 1e-6
max_iter = 20

solver = RootFinderOpen(f=f, fprime=fprime, x0=x0, x1=x1, tol=tol, max_iter=max_iter)
solver.show_toggle_open(interval_ms=750)

Non-converging example¶

f = lambda x: x**2 + 2
fprime = lambda x: 2 * x

x0 = 5.0
x1 = 1.0
tol = 1e-6
max_iter = 15

solver = RootFinderOpen(f=f, fprime=fprime, x0=x0, x1=x1, tol=tol, max_iter=max_iter)
solver.show_toggle_open(interval_ms=750)