# Week 3 worksheet 1: Root finding methods

This notebook is modified from one created by Charlotte Desvages.

This notebook investigates the first root finding algorithms to solve nonlinear equations: the **bisection method** and **fixed-point iteration**.

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-W1
```

---
### 📚 Book sections

- **ASC**: sections 5.1, **5.2**
- **PCP**: sections 7.1.3, 7.4

## Bracketing methods: bisection and regula falsi

**Bracketing methods** seek to find smaller and smaller **intervals** $[a, b]$ which contain a root.

They rely on the Intermediate Value Theorem: for a continuous function $F(x)$, if $F(a)$ has different sign than $F(b)$, then $F(x)$ has a root $x_\ast \in [a, b]$.

### Bisection

Key idea: start with an interval $[a, b]$ such that $F(a)$ and $F(b)$ have different signs.

Then, split the interval in two, and evaluate $F(\frac{a+b}{2})$ at the midpoint. Compare the sign of $F(\frac{a+b}{2})$ with the sign of $F(a)$ and $F(b)$, and choose the half-interval which still contains the root. Repeat the procedure with the new, smaller interval, until convergence.

### Regula falsi

Similar to bisection, but instead of choosing the midpoint, we choose the point of intersection between the x-axis and a straight line passing through the points $(a, F(a))$ and $(b, F(b))$.

---
## Solving nonlinear equations with root-finding algorithms

In this section, we consider equations of the form

$$
F(x) = 0,
$$

where $F(x)$ is a **nonlinear** function of $x$. Solving this equation is equivalent to finding the *root(s)* $x_\ast$ of the function $F$.

There are direct methods we can use to solve *linear* equations, even linear systems of equations as we saw previously; however, a nonlinear equation of the form given above doesn't always have a solution which can be found analytically.

The methods we will discuss to solve nonlinear equations are all **iterative**:
1. we start with a guess for where the root may be;
2. if we are close enough to the solution, for instance if the function is sufficiently close to zero for the current guess, we stop;
3. if not, we refine our guess using information we have about the function;
4. we go back to step 1 with our new guess.

Step 3 is what differentiates the methods we will use, as they each use a different process to refine the current best guess. For all these methods, the key idea is to reduce the nonlinear problem to smaller, simpler problems, which we solve repeatedly (iteratively).

### The bisection method

Given a continuous function $F \left( x \right)$, if $F \left( a \right) \le 0$ and $F \left( b \right) \ge 0$, the Intermediate Value Theorem tells us that there must be a root in the closed interval $\left[ a, b \right].$ The bisection method proceeds by testing the **sign** of $F \left( c \right)$ where $c$ is the **mid-point** of the interval $\left[ a, b \right]$, and uses this to halve the size of the interval in which a root is sought. The process is repeated with the new half-interval, until the interval is small enough to approximate the root with a given tolerance.

The next few exercises will guide you through implementing the bisection method yourself.

---
**📚 Learn more:**
- **ASC**: section 5.2
- **PCP**: section 7.4

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

\begin{equation}\label{eqn:F}
  F \left( x \right) = \sin \left( 2 \pi x \right) e^{4 x} + x.
\end{equation}

Plot this function in the interval $x \in \left[ 0, 1 \right]$ and identify the three roots in this interval. Check that $F \left( 0.2 \right)$ and $F \left( 0.6 \right)$ have opposite signs.

You may find it convenient to create a function `F`.


```python

```


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

---
🚩 ***Exercise 2:*** Define variables corresponding to $a = 0.2$ and $b = 0.6$. Define $c$ to be the midpoint $c = (a + b) / 2$.

Then, start a loop, which should iterate until the root is found. At each iteration:
- Depending on the sign of $F \left( a \right) F \left( c \right)$, decide whether to set $a$ or $b$ to be equal to $c$, so that there is a root of $F \left( x \right)$ in the *new* interval $\left[ a, b \right]$ (with half the width of the previous interval).
- Define $c$ to be the new midpoint $c = (a + b) / 2$.
- The loop should stop when you have found the root $x_\ast$ within an error of $10^{-10}$.

A possible convergence criterion is $|F(c)| < \varepsilon$, where $\varepsilon$ is the tolerance -- here, $10^{-10}$.

How many iterations are needed to find the root to within this error?


```python

```


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

---
🚩 ***Exercise 3:*** Choose different $a$ and $b$ values in order to find the root near $x = 1$, to within an error of $10^{-10}$.

You may wish to write your code from Exercise 2 into a function `bisection(F, a, b, tol)`, which finds the root of a function `F` in the interval `[a, b]` to within an error of `tol`, and returns the value of the roots and the number of iterations.


```python

```


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

---
### Regula falsi

The "regula falsi" method is similar to the bisection method, with an important difference: instead of selecting the midpoint of the interval $[a, b]$ at each iteration, we trace a straight line between the points $(a, F(a))$ and $(b, F(b))$, and select the point $c$ where this line intersects the x-axis. In other words, we interpolate $F$ linearly between $a$ and $b$, and we find the root of this interpolating polynomial (of degree 1) at each iteration.

---
🚩 ***Exercise 4:*** Show that $c$ is given by

$$
c = \frac{a F(b) - b F(a)}{F(b) - F(a)}.
$$

*Hint: a line with slope $\alpha$ which passes through the point $(x_0, y_0)$ has equation*

$$
y - y_0 = \alpha (x - x_0).
$$


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

---
🚩 ***Exercise 5:*** Consider the same function $F$ as in section 2.1.

Define variables corresponding to $a = 0.2$ and $b = 0.6$.

Then, proceed as you did for the bisection method, but instead of defining $c$ to be the midpoint of $[a, b]$, define $c$ as above.

How many iterations are needed to find the root to within a tolerance of $10^{-10}$?

You may wish to define a function `regula_falsi(F, a, b, tol)` to find a root of a function `F` within an interval `[a, b]` to within an error `tol`, which returns the computed root and the number of iterations.


```python

```


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

---
## Convergence of root-finding methods

The bisection and regula falsi methods are guaranteed to converge to a root, provided $F$ is sufficiently smooth and the starting interval $[a, b]$ is chosen appropriately.

But different methods may converge to a root at different *speeds*. The **order of convergence** for root-finding algorithms is defined in terms of successive values of the error $e_k := x_k - x_\ast$ between the true solution $x_\ast$ and the guess $x_k$ obtained at the $k$th iteration.

---
#### 🚩 Definition: Order of convergence of root-finding methods

A convergent root-finding algorithm converges **at $p$th order** if

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

where $\alpha \in \mathbb{R}$ is a constant.

---

For a $p$th order convergent method, we expect the error at the $k+1$th iteration to be roughly proportional to the $p$th power of the error at the $k$th iteration, for sufficiently large $k$ -- that is, when we are in a close enough neighbourhood of $x_\ast$.

Note that $p$ is not always an integer.

---
🚩 ***Exercise 6:*** Modify your code from Exercise 5 so that all the successive guesses $x_k$ are stored in a Numpy array. Perform the task from Exercise 6 again -- use the regula falsi method to find the same root of `F`, using the same starting interval and tolerance. You should obtain the same result, but now you should also have a vector `x` with length $k_\max + 1$ containing all the guesses.

Consider $x_{k_\max}$, the last guess obtained by the method, to be the "true solution". Compute the magnitude of the error $e_k$ between each of the previous guesses $x_k$ and the true solution.

For $p=1, 1.5, 2, 2.5$, compute the ratio $\frac{|e_{k+1}|}{|e_k|^p}$ for $k=0, 1, \dots, k_\max - 1$, and plot it against $k$. Set your y-axis limits to $[0, 1000]$ to start with, and reduce the range as necessary.

Given the definition above, what do you expect is the order of convergence of regula falsi? How do you explain the appearance of the graph?


```python

```


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

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