---
author:
- "by [Eugeniy E. Mikhailov](http://physics.wm.edu/~evmik/) and [Greg
  Bentsen](https://physics.wm.edu/~gbentsen/)"
title: Root finding algorithms
---

## Logistics and Agenda

Homework 3 due Monday Sept 21 at 11:59pm

Agenda

-   This week: Root-finding, numerical integration
-   Next week: Random numbers and distributions

Today

-   Root finding methods

## Midterm 1:

Do it at home part

-   You will be given a task which would require to splice couple
    covered concepts

-   Organize yourself in groups of 2 or less

    -   give us group composition by Friday Sep 25

-   You have almost two weeks to work on assignment

    -   prepare code
    -   report

-   Midterm write up and code is due Sunday October 4th 11:59pm

-   You may use AI freely to write code, debug, and test.

-   **You are responsible for the code you submit.**

-   We expect you to understand what the code is doing.

-   You may **not** use AI to write the report.

    -   Report should be done by you and only one per group is needed

In class part

-   on Monday October 5th, with computers but no AI
    -   we will give you modified assignment which would be similar
    -   it will require modification or reuse of prepared code

# Root finding problem

## Definition and canonical form

We want to find root of equation ($x_0$) which satisfies

$$f(x_0) = 0$$

Often it would work a bit different

$$h(x_0) = g(x_0)$$

Which can be moved to canonical form

$$f(x_0) = h(x_0) - g(x_0) = 0$$

Example

$$2\sin(x_0) = 1 \to 2\sin(x_0) - 1 = 0$$

## Trial and error method

One can try to get the solution by just guessing with a hope to hit the
solution. This is not highly scientific.

However, each guess can provide some clues.

A general algorithm

-   decide on required precision (stop conditions)
-   make a guess ($x_i$)
-   make **intelligent** new guess ($x_{i+1}$) which uses available
    information
-   repeat as long as $|f(x_{i+1})| > \varepsilon_f$ and
    $|x_{i+1} - x_i| > \varepsilon_x$

## Bisection method pseudocode

::: columns
::: {.column width="30% "}
Works for any `\alert{continuous function}`{=tex} in vicinity of a
function root

-   make initial bracket for search, i.e. set $x_+$ and $x_-$ such that
    -   $f(x_+) >0$
    -   $f(x_-) <0$
-   loop begins
    -   make the new guess value $x_g=(x_+ + x_-)/2$
    -   if $|f(x_{g})| \le \varepsilon_f$ and
        $|x_{+}-x_{g}| \le  \varepsilon_x$\
        then stop, we found the solution with the desired precision
    -   otherwise if $f(x_{g})>0$ then $x_+=x_{g}$ else $x_-=x_{g}$
    -   continue the loop
:::

::: {.column width="70%"}
![Bisection method
illustration](./pic_sources/bisection_animation.mp4){width="100%"}
:::
:::

## Bisection python implementation

AI uses prompt in
[bisection_ai_assignment.txt](./code/bisection_ai_assignment.txt)

``` python
#!/usr/bin/env python3
"""
Bisection method implementation.

The function bis(f, x_m, x_p) finds a root of f in the interval [x_m, x_p]
assuming f(x_m) < 0 and f(x_p) > 0.

Parameters:
    f: callable, function whose root is to be found.
    x_m: float, lower bound of the bracket (f(x_m) < 0).
    x_p: float, upper bound of the bracket (f(x_p) > 0).
    eps_f: float, tolerance for function value (default 1e-6).
    eps_x: float, tolerance for interval width (default 1e-6).

Returns:
    x0: float, approximation of the root.
    f(x0): float, function value at the root approximation.
"""

def bis(f, x_m, x_p, eps_f=1e-6, eps_x=1e-6):
    """Bisection method to find a root of f in the bracket [x_m, x_p].

    The algorithm iteratively narrows the bracket until both the function
    value and the bracket width are within the specified tolerances.

    Args:
        f: callable, function whose root is to be found.
        x_m: float, lower bound of the bracket (f(x_m) < 0).
        x_p: float, upper bound of the bracket (f(x_p) > 0).
        eps_f: float, tolerance for function value (default 1e-6).
        eps_x: float, tolerance for interval width (default 1e-6).

    Returns:
        x0: float, approximation of the root.
        f(x0): float, function value at the root approximation.
    """
    x_m= x_m
    x_p= x_p
    while True:
        x_g = (x_p+ x_m) / 2.0
        f_g = f(x_g)
        if abs(f_g) <= eps_f and abs(x_p- x_g) <= eps_x:
            return x_g, f_g
        if f_g > 0:
            x_p= x_g
        else:
            x_m= x_g

if __name__ == "__main__":

# Simple test: find sqrt(2) as root of f(x) = x^2 - 2

def f(x):
        return x * x - 2.0
    root, val = bis(f, 1.0, 2.0)
    print(f"Root approximation: {root}")
    print(f"f(root) = {val}")
    print(f"x mismatch = {abs(root - 2**0.5)}")
```

# Solution convergence

## Solution convergence expression

We say that algorithm has defined convergence if it is possible to
express

$$\lim_{k\to\infty} (x_{k+1}-x_0)=c (x_k-x_0)^m$$

Where $x_0$ is true root of the equation, $c$ is some constant, and $m$
is the order of convergence.

The best algorithms have quadratic convergence, i.e. $m=2$

-   The bisection algorithm has the linear rate of convergence: ($m=1$)
    and $c=1/2$

Generally the speed of the algorithm is related to its convergence
order. However, other factors may affect the speed.

-   For some algorithm it is hard to define convergence parameters

## Newton-Raphson method

::: columns
::: {.column width="30% "}
$$x_{i+1}=x_i-\frac{f(x_i)}{f'(x_i)}$$

Need to provide a starting points $x_1$ and the derivative of the
function.

Newton-Raphson method converges quadratically ($m=2$).
:::

::: {.column width="70%"}
![Newton-Raphson method
illustration](./pic_sources/newton_raphson_animation.mp4){width="100%"}
:::
:::

# Numerical derivative of a function

## Numerical derivative of a function

Mathematical definition

$$f'(x)=\lim_{h \to 0}\frac{f(x+h)-f(x)}{h}$$

The initial intent is to calculate it at very small $h$.

Remember about round off errors

For computers with $h$ small enough $f(x+h)-f(x)=0$.

Let's be smarter. Recall Taylor series expansion

$$f(x+h)=f(x)+\frac{f'(x)}{1!}h+\frac{f''(x)}{2!}h^2+\cdots$$

So we can see

$$f_c'(x)=\color{blue}{ \frac{f(x+h)-f(x)}{h}} = f'(x) \color{red}{ +\frac{f''(x)}{2}h +\cdots}$$

Here computed approximation and algorithm error.

There is a range of optimal $h$ when both the round off and the
algorithm errors are small.

## Derivative via Forward Difference (fd)

$$\color{blue}{
    f_{fd}'(x) = \frac{f(x+h)-f(x)}{h}
    }$$

Algorithm error for small $h$

$$\varepsilon_{fd} \approx \frac{f''(x)}{2} h$$

This is quite bad since error is proportional to $h$.

::: example
$$f(x) = a+b x^2$$

$$f(x+h) = a+b (x+h)^2  = a + b x^2 + 2 b x h + b h^2$$

$$\color{blue}{f_{fd}'(x)} = \frac{f(x+h)-f(x)}{h} = { 2 b x } \color{red}{ + b h }$$

So for small $x$, the algorithm error dominates our approximation!
:::

## Derivative via Central Difference (cd)

This is much better way (though it use an extra computation for a
function)

$$\color{blue}{ f_c'(x) = \frac{f(x+h)-f(x-h)}{2h}}$$

Algorithm error $$\varepsilon_{cd} \approx \frac{f'''(x)}{6}h^2$$

# Root finding Python library

One option is to use `root_scalar`{.python} from
`scipy.optimize`{.python} library

``` python
def f(x): return (x**3-x-2)
from  scipy.optimize import root_scalar

# braketing method

sol=root_scalar(f, bracket=[1,2])
In [18]: sol.root
Out[18]: 1.5213797068045676

# single point guess

sol=root_scalar(f, x0=2)
In [22]: sol.root
Out[22]: np.float64(1.5213797068045676)
```