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