Toggle navigation

.ipynb .pdf

Root-finding solver based on Brent’s method

Now we start with our Brent root-finding algorithm example code.

We start with untyped, script-like code, and modularize it bit by bit.

This document is a Jupyter notebook! Thus it shows how to use MyST-NB to include notebooks in a Sphinx site.

from math import cos

x_rel_tol = 1e-12  # X relative tolerance
x_abs_tol = 1e-12  # X absolute tolerance
y_tol = 1e-12  # Y tolerance
verbose = True  # Whether to display progress

f = cos
a = 0.0  # first point that brackets the root
b = 3.0  # second point that brackets the root
fa = f(a)  # f(a)
fb = f(b)  # f(b)

if abs(fa) < abs(fb):
    # force abs(fa) >= abs(fb), make sure that b is the best root approximation known so far
    # and a is the contrapoint
    b, a = a, b
    fb, fa = fa, fb

c = a  # Previous iterate
fc = fa  # f(c)
d = a  # Iterate before the previous iterate
fd = fa  # f(d)
last_step = None
step = None
iter = 1  # Current iteration number (1-based)

# display the initial state
(iter, b, fb, abs(a-b), abs(fa-fb), last_step)

(1, 3.0, -0.9899924966004454, 3.0, 1.9899924966004454, None)

while True: # loop forever
    iter = iter + 1
    last_step = step
    if fa != fc and fb != fc:
        # perform a quadratic interpolation step
        # https://en.wikipedia.org/wiki/Inverse_quadratic_interpolation
        L0 = (a * fb * fc) / ((fa - fb) * (fa - fc))
        L1 = (b * fa * fc) / ((fb - fa) * (fb - fc))
        L2 = (c * fb * fa) / ((fc - fa) * (fc - fb))
        s = L0 + L1 + L2
        step = "quadratic"
    else:
        # perform a secant step
        s = b - fb * (b - a) / (fb - fa)
        step = "secant"
    perform_bisection = False
    if a <= b and not ((3 * a + b) / 4 <= s <= b):
        perform_bisection = True
    elif b <= a and not (b <= s <= (3 * a + b) / 4):
        perform_bisection = True
    elif last_step == "bisection" and abs(s - b) >= abs(b - c) / 2:
        perform_bisection = True
    elif last_step != "bisection" and abs(a - b) >= abs(c - d) / 2:
        perform_bisection = True
    elif last_step == "bisection" and abs(b - c) < x_abs_tol:
        perform_bisection = True
    elif last_step != "bisection" and abs(c - d) < x_abs_tol:
        perform_bisection = True
    if perform_bisection:
        # perform a bisection step
        s = min(a, b) + abs(b - a) / 2
        step = "bisection"
    fs = f(s)
    d = c
    c = b
    fc = fb
    # check which point to replace to maintain (a,b) have different signs
    if f(a) * f(s) < 0:
        b = s
        fb = fs
    else:
        a = s
        fa = fs
    # keep b as the best guess
    if abs(fa) < abs(fb):
        b, a = a, b
        fb, fa = fa, fb

    # Checks convergence
    if fb == 0:
        print("Exact root found")
        break
    x_delta = abs(a - b)
    if x_delta <= x_abs_tol:
        print("Met x_abs_tol criterion")
        break
    if x_delta / max(a, b) <= x_rel_tol:
        print("Met x_rel_tol criterion")
        break
    y_delta = abs(fb)
    if y_delta <= y_tol:
        print("Met y_tol criterion")
        break

    dx = abs(a - b)
    dy = abs(fa - fb)
    print(f"{iter}\t{b:.3e}\t{fb:.3e}\t{dx:.3e}\t{dy:.3e}\t{last_step}")

print(f"The approximate root is {b}")

2	1.500e+00	7.074e-02	1.500e+00	1.061e+00	None
3	1.600e+00	-2.923e-02	1.000e-01	9.997e-02	bisection
4	1.571e+00	1.434e-05	2.925e-02	2.924e-02	secant
5	1.571e+00	-2.042e-09	1.434e-05	1.434e-05	secant
Met y_tol criterion
The approximate root is 1.5707963267948966