# # Week 3 worksheet 2: Fixed point iteration and Newton's method

This notebook is modified from one created by Charlotte Desvages.

This week, we continue exploring **root-finding methods**, this time looking at **fixed-point iteration** and **Newton's method**.

The best way to learn programming is to write code. Don't hesitate to edit the code in the example cells, or add your own code, to test your understanding. You will find practice exercises throughout the notebook, denoted by 🚩 ***Exercise $x$:***.

#### Displaying solutions

Solutions will be released one week after the worksheets are released, as a new `.txt` file in the same GitHub repository. After pulling the file to your workspace, run the following cell to create clickable buttons under each exercise, which will allow you to reveal the solutions.


```python
%run scripts/create_widgets.py W03-W2
```

---
### 📚 Book sections

- **ASC**: sections **5.3**, **5.4**, **5.6**
- **PCP**: section 7.2

---

## 1. Fixed-point iteration

In Week 8, we implemented the bisection and regula falsi methods to find solutions of nonlinear equations; they are both *bracketing* methods. The methods we will see this week are different: they all start with an initial guess $x_0$, a single value (not an interval), and seek to refine this guess iteratively to eventually converge to the solution.

### 1.1. General construction

Fixed-point iteration seeks a root to a function $F \left( x \right)$ by defining a sequence

$$
  x_{k + 1} = G \left( x_k \right) \quad \textrm{ for } k = 0, 1, 2, \ldots
$$

where $G \left( x_k \right)$ is an **iteration function**, designed so that

$$
F \left( x_* \right) = 0 \quad \text{if and only if} \quad x_* = G \left( x_* \right).
$$

That is, the roots of $F \left( x \right)$ are *fixed-points* of the sequence, and conversely the fixed-points of the sequence are roots of $F \left( x \right)$.

Given an "initial guess" $x_0$, if the sequence converges, i.e. if

$$
  \lim_{k \rightarrow \infty} \left| x_{k + 1} - x_k \right| = 0,
$$

then it must, by construction, converge to a root of the original function $F \left( x \right)$.

---
🚩 ***Exercise 1***: Consider the quadratic polynomial

$$
F(x) = -x^2 - 3x + 2.
$$

Plot $F(x)$ over $[-5, 2]$. Compute and display the exact value of both roots $x_\ast$ of $F(x)$.


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex1
```

---
🚩 ***Exercise 2***: The roots $x_\ast$ both satisfy

$$
F(x_\ast) = -x_\ast^2 - 3x_\ast + 2 = 0.
$$

Two of the many possible ways to rearrange this equation to something of the form $x_\ast = G(x_\ast)$ are:

\begin{align}
x_\ast &= G_1(x_\ast) = \frac{2}{x_\ast + 3}, \\
x_\ast &= G_2(x_\ast) = -x_\ast^2 - 2x_\ast + 2.
\end{align}

Define two functions `G1(x)`, `G2(x)` corresponding to the above. Then, with a starting guess of $x_0 = -4$, compute $x_{k+1} = G_i(x_k)$ for $k=1, 2, \dots, 7$, and $i=1, 2$. Display all successive values of the guesses $x_k$ for both methods.

Do you observe convergence for both methods? If so, to which root?


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex2
```

---
### 1.2. Conditions for convergence

Although there are often many ways to write $F(x) = 0$ as $x = G(x)$, not all of these lead to convergence.

---
#### 🚩 Theorem: Conditions for convergence of fixed point iteration

Let $x_\ast$ be a solution of $F(x) = 0$, and an iteration function $G(x)$ be defined such that $x_\ast$ is a fixed point of $G$, i.e. we have $x_\ast = G(x_\ast)$.

If $G(x)$ is differentiable, the fixed-point iteration $x_{k+1} = G(x_k)$ will converge to $x_\ast$ if
- $|G'(x)| < 1$ in some neighbourhood of $x_\ast$,
- and either:
    - the initial guess $x_0$ is chosen in this neighbourhood, or
    - some later iteration value $x_k$ is in this neighbourhood.
    
---
The proof is given in **ASC section 5.3**, and can be established using the Mean Value Theorem (more about this in Video 1 this week) to write the error $e_{k+1} = x_{k+1} - x_\ast$ in terms of $e_k$, the error at the previous iteration:

$$
|e_{k+1}| = |G'(\eta)| |e_k|, \qquad \text{for some} \: \eta \in (\min(x_k, x_\ast), \max(x_k, x_\ast)).
$$

As long as $|G'(x)| < 1$ in a neighbourhood of $x_\ast$ (which includes $\eta$), then for any $x_k$ in that neighbourhood, we have $|e_{k+1}| < |e_k|$, meaning that the absolute error decreases monotonically for all subsequent iterations, guaranteeing convergence to $x_\ast$.

---
🚩 ***Exercise 3***: Plot $|G_1'(x)|$ and $|G_2'(x)|$ over $[-5, 5]$, with the exact roots of $F(x)$ indicated on the same plot. For either method, are there possible initial guesses which guarantee convergence towards either root?


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex3
```

### 1.3. Convergence rate for fixed point iteration

The next few exercises explore the rate of convergence of fixed point iteration.

---
🚩 ***Exercise 4***: Consider the function

\begin{equation}
  F \left( x \right) = x - x^2 \sin x.
\end{equation}

Plot this function in the interval $x \in \left[ -2, 2 \right]$ and identify the three roots in this interval. 


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex4
```

---
🚩 ***Exercise 5***: At a root $x_*$ of $F \left( x \right)$ we have

\begin{equation}
  x_* - x_*^2 \sin x_* = 0.
\end{equation}

Given an $\alpha \ne 0$ this is equivalent to

\begin{equation}
  x_* = \left( 1 + \alpha \right) x_* - \alpha x_*^2 \sin x_*.
\end{equation}

This suggests the fixed-point iteration

\begin{equation}\label{eqn:fp2}
  x_{n + 1} = \left( 1 + \alpha \right) x_n - \alpha x_n^2 \sin x_n \quad \textrm{ for } n = 0, 1, 2, \ldots.
\end{equation}


Use a `while` loop to implement fixed-point iteration using the scheme above with an initial guess $x_0 = 1$ and with $\alpha = 1.2$.  You will need to choose an appropriate condition on which to exit the loop.

Hint: you may find it convenient to define a function `G` representing the iteration function.


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex5
```

---
🚩 ***Exercise 6***: Use the function `fsolve()` from `scipy.optimize` to find the root of $F \left( x \right)$ near $x = 1$. Create a plot of error magnitude against iteration number, where the error is the difference between the solution obtained via `fsolve()` (taken as the "ground truth"), and the values for $x_k$ obtained via the fixed-point iteration used in Exercise 5. Use a logarithmic scale for the y-axis.


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex6
```

---
🚩 ***Exercise 7***: Find the root near $x = 1$ using bisection, with an initial interval defined by $a = 0.5$ and $b = 1.5$. In the same figure as the plot for Exercise 6, plot the magnitude of the difference between the root obtained via `fsolve()`, and the values of $\left( a + b \right) / 2$ obtained in each bisection iteration, plotting this error magnitude against iteration number. Which of bisection and fixed point iteration converges faster to the root?


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex7
```

---
🚩 ***Exercise 8***: Repeat Exercise 6 with $\alpha = 0.8$. Plot the error magnitude on the same graph as that from Exercises 6 and 7. What do you observe?


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex8
```

---
## 2. Newton's method

Newton's method proceeds by repeatedly finding the root $x_{k+1}$ of a linear approximation to the function $F$ around the current guess $x_k$. The algorithm proceeds with the following steps until convergence:
1. Start from a guess $x_k$.
2. Draw the tangent to $F(x)$ at $x = x_k$.
3. Find the point of intersection of this tangent with the x-axis. This is the next guess $x_{k+1}$.
4. Advance to the next iteration, i.e. set $x_k = x_{k+1}$ and go back to step 1.

| Initial guess ($x_0 = 0.5$) | Drawing the tangent |
|:-:|:-:|
|  |  | 
| **New guess $x_1$ - first iteration complete** | **Second iteration: draw the tangent at $x = x_1$ to find $x_2$** |
|  |  |

### 2.1. Derivation

Consider the Taylor expansion of $F(x)$ around $x_k$:

$$
F(x) = F(x_k) + (x - x_k) F'(x_k) + \dots
$$

The goal is to find a root of $F$ at a future iteration. Ignoring the higher order terms in the Taylor expansion lets us approximate $F$ with a linear function -- the tangent of $F$ at $x = x_k$:

$$
F(x) \approx F(x_k) + (x - x_k) F'(x_k).
$$

Newton's method then proceeds by finding the root of this linear approximation, i.e. the point of intersection between the tangent at $x = x_k$ and the x-axis, by setting

$$
F(x) \approx F(x_k) + (x - x_k) F'(x_k) = 0.
$$

Solving the above for $x$ gives us the next guess $x_{k+1}$:

$$
x_{k+1} = x_k - \frac{F(x_k)}{F'(x_k)}.
$$

---

Note that **Newton's method** is a particular type of **fixed-point iteration**. Starting from an initial guess $x_0$, we have

$$
x_{k+1} = G(x_k) \quad \textrm{ for } k = 0, 1, 2, \ldots
$$

where the iteration function $G(x)$ takes the form

$$
G(x) = x - \frac{F(x)}{F'(x)}.
$$

---
🚩 ***Exercise 9***:
Consider the function

$$
  F \left( x \right) = e^{0.1 x} \sin \left( 4 \pi x \right) + x^2 + 0.5.
$$

Define a function `F` to represent $F \left( x \right)$. Use this function `F` to plot $F(x)$ in the interval $x \in \left[ 0, 1 \right]$, and identify the approximate values of the two roots in this interval.


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex9
```

---
🚩 ***Exercise 10***: Define a second function `Fp` which can be used to evaluate the derivative $F' \left( x \right)$.

Then, define a function `G` which can be used to evaluate the iteration function associated with using Newton's method to find a root of $F \left( x \right)$.

Plot $G \left( x \right) - x$. By visual comparison with the plot from Exercise 9, verify graphically that, for $x \in \left[ 0, 1 \right]$, $G \left( x \right) = x$ at the same values of $x$ at which $F \left( x \right) = 0$ (i.e. that the fixed points of $G(x)$ are the same as the roots of $F(x)$).


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex10
```

---
🚩 ***Exercise 11***: For each of the two roots, choose an appropriate initial guess, and use Newton's method to compute them.


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex11
```

---
### 2.2. Convergence of Newton's method

Newton's method is a particular type of fixed point iteration. As we've seen in Section 1.2 earlier, the absolute error at iterations $k$ and $k+1$ are related by

$$
|e_{k+1}| = |G'(\eta)| |e_k|, \qquad \text{for some} \: \eta \in (\min(x_k, x_\ast), \max(x_k, x_\ast)).
$$

Consider the case where $F$ is twice differentiable (at least in a neighbourhood of $x_\ast$). What is the value of $G'(x)$ at the root $x = x_\ast$?

$$
G'(x_\ast) = 1 - \frac{F'(x_\ast)}{F'(x_\ast)} + \frac{F(x_\ast) F''(x_\ast)}{\left(F'(x_\ast)\right)^2}
= 0, \qquad \text{since by definition,} \: F(x_\ast) = 0.
$$

Therefore, if $G(x)$ is *continuously differentiable* in some neighbourhood of $x_\ast$, we have

$$
\lim_{x \to x_\ast} G'(x) = 0.
$$

This means that for $\eta$ in this neighbourhood, $|G'(\eta)|$ decreases more and more as we approach the root, and convergence **accelerates**!

---
#### Conditions for convergence of Newton's method

To guarantee convergence for fixed-point iteration, we saw earlier that it is sufficient to choose an initial guess $x_0$ in a neighbourhood of $x_\ast$ where $G(x)$ is differentiable and $|G'(x)|<1$.

For Newton's method, we have $G'(x_\ast) = 0$. Therefore, if $G'(x)$ is continuously differentiable in some neighbourhood of $x_\ast$, then by the Intermediate Value Theorem applied to $G'(x)$, there must be a neighbourhood (possibly smaller) near the root for which $|G'(x)| < 1$.

Therefore, for Newton's method to converge, it is sufficient that
- $G(x)$ is continuously differentiable,
- and either
    - $x_0$ is chosen to be sufficiently close to $x_\ast$, or
    - some later iteration value $x_k$ is sufficiently close to $x_\ast$.

This ensures that all subsequent iterations after $x_0$ (or after $x_k$) will reduce the absolute error.

---
#### Order of convergence of Newton's method

In Week 8 we saw that the order of convergence $p$ of a root-finding method verified

$$
\lim_{k\to \infty} \frac{|e_{k+1}|}{|e_k|^p} = \alpha,
$$

where $\alpha$ is a constant.

---
🚩 ***Exercise 12***: Consider a function

$$
  F \left( x \right) = e^x - a,
$$

where $a$ is a given positive real number. The root of this equation is $x_\ast = \log a$.

Implement Newton's method to compute $\log 5$ as the root of $F(x)$ for $a = 5$, starting from an initial guess of $x_0 = 1$. How many iterations does Newton's method take until $\left| x_{k + 1} - x_k \right| < 10^{-14}$?


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex12
```

---
🚩 ***Exercise 13***: Modify your code from Exercise 12 to store the successive values of the absolute error $|e_k| = |x_k - x_\ast|$. Use these values to determine numerically the order of convergence of Newton's method.


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex13
```

---
🚩 ***Exercise 14***: Compute and display the successive values of the absolute error using different values of the initial guess $x_0$. How many iterations are needed to converge to the root for different $x_0$? What happens when $x_0$ is set further away from $x_\ast$?


```python

```


```python
%run scripts/show_solutions.py W03-W2_ex14
```

147e159ef4cb82138f84cee675b593564565708a
ipynb
Jupyter Notebook
Workshops/W03-W2_NMfCE_Newtons_method.ipynb
DrFriedrich/nmfce-2021-22
2ccee5a97b24bd5c1e80e531957240ffb7163897
MIT
Workshops/W03-W2_NMfCE_Newtons_method.ipynb
DrFriedrich/nmfce-2021-22
2ccee5a97b24bd5c1e80e531957240ffb7163897
MIT
Workshops/W03-W2_NMfCE_Newtons_method.ipynb
DrFriedrich/nmfce-2021-22
2ccee5a97b24bd5c1e80e531957240ffb7163897
MIT
Qwen/Qwen-72B
1. YES
2. YES