WS09: Robust nonlinear solvers#

These exercises are indented to give you practice at using the material on numerical approximation and are intended to reinforce the material that was covered in lectures.

Please attempt the worksheet before your tutorial. Support is available in your tutorial or in the Class Team.

Partial solutions are available to check your understanding. Please create Issues and Pull requests with your solutions or share your answers for peer feedback in the class team Worksheet solutions channel

Part a (pen and paper warm up)#

1. Convergence of Newton’s method#

For

\[ f(x) = x(x-3)^2 \]

with

\[ f^\prime(x) = 3(x-1)(x-3), \]

we know \(f(x) = 0\) has two solution \(x=0\) and \(x=3\). You might plot this function in the following Q4 to help you visualise it.

Starting form initial guesses \(x_0 = 0.5\) and \(x_0=1.5\) respectively, and apply two Newton’s iterations to see which solution you will get.

2. Problems of Newton’s method#

Create an equation so that Newton’s method does not work, and prove it using your implementation in WS08.

3. The secant method#

Use the secant method to solve the above equation \(f(x)=0\), and iterate three times from two initial guesses \(x_0=4\), \(x_1=5\).

Part b (code implementations and testing)#

4. Plot the function f(x) in Q1#

import numpy as np


def f(x):
    return x * x * x - 6.0 * x * x + 9.0 * x


X = np.linspace(-1.0, 4.0, num=100, endpoint=True)
Y = np.zeros(100)
"""
Add a loop to compute array Y=f(X)
"""

"""
Add code to plot Y=f(X)
"""

5. Implement the secant method#

import numpy as np

def secant(f, x0, x1, tol):
    x = x1
    it = 0
    while abs(f(x)) > tol:  # iterate until less than or eq tol
        """
        TODO: implement the secant method
        """

    return x, it

6. Test your implementation#

x, it = secant(f, 4.0, 5.0, 1.0e-6)
print(f"The secant method: {x} after {it} iterations")

Part c: Extension#

Run the following SciPy solver, which struggles to find the solution \(x=3\) for equation \(f(x) = x(x-3)^2 =0\).

from scipy import optimize

root = optimize.brentq(f, -1, 10, full_output=True)
print(root)