Skip to content
scipython.com

P8.4.2Root-finding with newton and brentq

Find all of the roots of

\[f(x) + \frac{1}{5} + x\cos\left(\frac{3}{x}\right).\]

using (a) scipy.optimize.brentq and (b) scipy.optimize.newton

Show solution

To bracket the roots, we can evaluate the function on a grid of points and find where it changes sign. In the code below, lower_bounds are the values of \(x=a\) on this grid immediately before a sign change, so the root is bracketed by the tuple (a, a+dx) where \(dx\) is the grid point spacing.

The Newton-Raphson method works best if the analytical form for the first derivative can be given:

\[f'(x) = \cos\left(\frac{3}{x}\right) + \frac{3}{x}\sin\left(\frac{3}{x}\right).\]

This function is defined in the lambda expression fp.

rootfinding
import numpy as np
from scipy.optimize import brentq, newton

f = lambda x: 0.2 + x * np.cos(3 / x)

x, dx = np.linspace(-1, 1, 1000, retstep=True)

# Find where sign changes between consecutive points.
fvals = f(x)
signs = np.sign(fvals)
sign_changes = np.where(np.diff(signs) != 0)[0]
# Use these indices to get the x-values bracketing the roots.
lower_bounds = x[sign_changes]

# Brent's method
brent_roots = np.array([brentq(f, a, a + dx) for a in lower_bounds])

# Newton's method: requires first derivative of f(x)
fp = lambda x: np.cos(3 / x) + 3 / x * np.sin(3 / x)
newton_roots = np.array([newton(f, a, fp) for a in lower_bounds])
print("Roots of f(x) = 1/5 + x.cos(3/x):")
print(f"{'Brent':11s} {'Newton':11s}")
r = np.vstack((brent_roots, newton_roots)).T
for br, nr in r:
    print(f"{br:11.8f} {nr:11.8f}")

Output:

Roots of f(x) = 1/5 + x.cos(3/x):
Brent       Newton     
-0.59333063 -0.59333063
-0.40858197 -0.40858197
-0.25181946 -0.25181946
-0.22939031 -0.22939031
 0.29285950  0.29285950
 0.35492010  0.35492010
 0.67969870  0.67969870