WS09: Robust nonlinear solvers - partial solutions

WS09: Robust nonlinear solvers - partial solutions#

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.

These are partial solutions. Please create Issues and Pull requests with your solutions.

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 curve 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)
for i in range(len(Y)):
    Y[i] = f(X[i])

import matplotlib.pyplot as plt

fig = plt.figure()
plt.plot(X, Y)
plt.grid()
plt.show()

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
        x = x - f(x1) * (x1 - x0) / (f(x1) - f(x0))  # apply one Newton iteration
        x0 = x1
        x1 = x
        it = it + 1

    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")
np.testing.assert_allclose(abs(f(x)), 0.0, atol=1.0e-6)

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)

root = optimize.newton(f, 1.5, full_output=True)
print(root)
Contents