Using scipy.optimize.newton to find a root of the following functions (with the given starting point, $x_0$) fails. Explain why and find the roots either by modifying the call to newton or by using a different method.
(a) $$ f(x) = x^3 - 5x, \quad x_0 = 1 $$
(b) $$ f(x) = x^3 - 3x + 1, \quad x_0 = 1 $$
(c) $$ f(x) = 2-x^5, \quad x_0 = 0.01 $$
(d) $$ f(x) = x^4 - (4.29)x^2 -5.29, \quad x_0 = 0.8 $$
Some examples of root-finding for which the Newton-Raphson algorithm fails and how to solve this.
(a)
In [x]: newton(lambda x: x**3 - 5*x, 1, lambda x: 3*x**2 - 5)
...
RuntimeError: Failed to converge after 50 iterations, value is 1.0
The Newton-Raphson algorithm enters an endless cycle of values for $x$: \begin{align*} x_0 = 1:& \quad x_1 = x_0 - f(x_0)/f'(x_0) = -1\\ x_1 = -1: & \quad x_2 = x_1 - f(x_1)/f'(x_1) = 1\\ x_2 = 1: & \quad x_1 = x_0 - f(x_0)/f'(x_0) = -1\\ \cdots \end{align*} Alternative starting points converge correctly on a root. Even a very small displacement from $x=0$ ensures convergence:
In [x]: newton(lambda x: x**3 - 5*x, 1.0001, lambda x: 3*x**2 - 5)
Out[x]: 2.23606797749979
In [x]: newton(lambda x: x**3 - 5*x, 1.1, lambda x: 3*x**2 - 5)
Out[x]: -2.23606797749979
In [x]: newton(lambda x: x**3 - 5*x, 0.5, lambda x: 3*x**2 - 5)
Out[x]: 0.0
(b)
In [x]: f, fp = lambda x: x**3 - 3*x+1, lambda x: 3*x**2 - 3
In [x]: newton(f, 1, fp)
Out[x]: 1.0
In [x]: f(1.0)
Out[x]: -1
The algorithm converged, but not on a root! Unfortunately, the gradient of the function is zero at the chosen starting point and because of round-off error this has not led to a ZeroDivisionError. To find the roots, choose different starting points such that $f'(x_0) \ne 0$, or use a different method after bracketing the roots by inspection of a plot of the function:
In [x]: brentq(f, -0.5, 0.5), brentq(f, -2, -1.5), brentq(f, 1, 2)
Out[x]: (0.34729635533386066, -1.879385241571423, 1.532088886237956)
(c) The function $f(x) = 2-x^5$ has a flat plateau around $f(0)=2$ and the small gradient here leads to slow convergence on the root:
In [x]: newton(f, 0.01, fp)
...
RuntimeError: Failed to converge after 50 iterations, value is ...
To find it using newton either move the starting point closer to the root, or increase the maximum number of iterations:
In [x]: newton(f, 0.01, fp, maxiter=100)
Out[x]: 1.148698354997035
(d) This is another example of a function which generates an endless cycle of values from the Newton-Raphson method:
In [x]: f = lambda x: x**4 - 4.29 * x**2 - 5.29
In [x]: fp = lambda x: 4*x**3 - 8.58 * x
In [x]: newton(f, 0.8, fp)
...
RuntimeError: Failed to converge after 50 iterations, value is ...
Unlike the function in (a), the region $0.6 \le x_0 \le 1.1$ attracts this cyclic behaviour, so one needs to initialize the algorithm outside this range to obtain the roots $\pm 2.3$. For example,
In [x]: newton(f, 1.2, fp)
Out[x]: -2.3