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

# Root finding algorithm gotchas

## Diverging Newton-Raphson

::: columns
::: {.column width="55%"}
![Example of diverging Newton-Raphson
case](./pics/newton_gotchas.svg){width="100%"}
:::

::: {.column width="45% "}
Newton-Raphson and other using estimate of derivative algorithms are
usually fast but starting points need to be close enough to the root.

Watch out for

-   roots which are surrounded by region of low derivative of the
    function
-   changing sign of the slope near root
:::
:::

## Example of a tricky function $\tanh(2.8-x)$

::: columns
::: {.column width="55%"}
![Plot of $\tanh(2.8-x)$](./code/tanh_x.png){width="100%"}
:::

::: {.column width="45% "}
Using `root_scalar` from `scipy.optimize` package

``` python
root_scalar(f, args=(), method=None, bracket=None, fprime=None, fprime2=None, x0=None, x1=None, xtol=None, rtol=None, maxiter=None, options=None)
```

See file [find_tanh_root.py](./code/find_tanh_root.py)

``` python
import numpy as np
from scipy.optimize import root_scalar

def f(x, *my_arg, **my_kwargs): # note dummy variable needed by root_scalar
    return np.tanh(2.8 - x)

```

``` python
In [78]: root_scalar(f, 'newton', x0=2.5)
Out[78]:
      converged: True
           flag: converged
 function_calls: 7
     iterations: 3
           root: 2.8
         method: newton

In [79]: root_scalar(f, 'newton', x0=1)
/usr/lib/python3/dist-packages/scipy/optimize/_root_scalar.py:327:
    RuntimeWarning: Derivative was zero.
  r, sol = methodc(f, x0, args=args, fprime=fprime, fprime2=None,
Out[79]:
      converged: False
           flag: convergence error
 function_calls: 6
     iterations: 3
           root: -596803.8420472746
         method: newton

```
:::
:::

## Tricky function $\tanh(2.8-x)$ and bracketing algorithms

::: columns
::: {.column width="55%"}
![Plot of $\tanh(2.8-x)$](./code/tanh_x.png){width="100%"}
:::

::: {.column width="45% "}
``` python
In [85]: root_scalar(f, method='bisect', bracket=[1,5])
Out[85]:
      converged: True
           flag: converged
 function_calls: 43
     iterations: 41
           root: 2.8000000000010914
         method: bisect

In [86]: root_scalar(f, method='ridder', bracket=[1,5])
Out[86]:
      converged: True
           flag: converged
 function_calls: 10
     iterations: 4
           root: 2.8
         method: ridder
```

Note that this methods take longer and more function calls. But they
always work if bracketed right.
:::
:::

# Root finding algorithm summary

::: columns
::: {.column width="50%"}
**Bracketing algorithm**

-   usually slower
-   robust and **always converge**
-   require setting the bracket
    -   need a priori knowledge

::: example
-   bisection
-   false position
-   Ridders
-   Brent's method
:::
:::

::: {.column width="40% "}
**Non bracketing algorithm**

-   Usually fast
-   Can diverge
    -   sensitive to starting point

::: example
-   Newton-Raphson
-   secant
:::
:::
:::

# Integration gotchas

## Ends of integration have infinity

One or more ends of integral have infinity

$$ \int_a^\infty f(x) dx $$

Solution: do a variable change.

For example $$ x = a + \frac{t}{1-t}$$

$$ \int_a^\infty f(x) dx  \to 
\int_0^1 f(a + \frac{t}{1-t}) \frac{dt}{(1-t)^2}
$$

## Undersampling

::: columns
::: {.column width="55%"}
![Peaky function](./code/peaky_f.png){width="100%"}
:::

::: {.column width="45% "}
Not enough points

-   function has narrow features
-   rapid oscillation

Solution

-   more points are needed
:::
:::

## Function have discontinuity

::: columns
::: {.column width="55%"}
![Step function](./code/step_f.png){width="100%"}
:::

::: {.column width="45% "}
Simpson and similar methods, which assumes smooth function, will have
larger error than seemingly inferior methods

``` python
from scipy.integrate import trapezoid, simpson
import numpy as np

x=np.linspace(-1,1, 11)
y=0*x
y[x>0] = 1

In [69]: trapezoid(y, x)  # note order of x and y
Out[69]: np.float64(0.8999999999999999)

In [70]: simpson(y, x)
Out[70]: np.float64(0.8666666666666667)
```
:::
:::

## Function can be integrated but is infinite at some locations

::: columns
::: {.column width="55%"}
We know that integral exist but numerically it produce infinities,
typically because we divide by zero

::: example
For example: $$ \int_0^b \frac{1}{\sqrt{x}} dx = 2\sqrt{b} $$
:::

We can

-   skip the small vicinity around singularity
    $$ \int_\delta^1 \frac{1}{\sqrt{x}} dx$$

-   add tiny offset to denominator
    $$ \int_0^1 \frac{1}{\delta+\sqrt{x}} dx$$

The above not always work, since as $\delta \to 0$ we will have some
problem with the numerical evaluation. For example, Simpson method might
not converge.
:::

::: {.column width="45% "}
there are so called adaptive methods, which adjust number of point
around singularities

-   for example `quad` from `scipy.integrate`

``` python
In [13]: from scipy.integrate import quad
In [14]: def f(x): return 1/np.sqrt(x)

In [15]: quad(f, 0, 1)
Out[15]: (1.9999999999999984, 5.773159728050814e-15)
```
:::
:::

## Function can be integrated but is infinite at some locations continued

Alternatively, we can recall that "computer is not substitute for the
brain" and do a variable substitution trick

$$I = \int_a^b \frac{1}{\sqrt{f(x)}} dx$$

A more useful identity comes from differentiating $\sqrt{f(x)}$: $$
\frac{d}{dx}\sqrt{f(x)}
=
\frac{f'(x)}{2\sqrt{f(x)}}.
$$ Hence, provided $f'(x)\neq0$, $$
\boxed{
\frac1{\sqrt{f(x)}} = \frac{2}{f'(x)} \frac{d}{dx}\sqrt{f(x)}.
}
$$

Substituting this into the integral, $$
I = \int_a^b \frac{2}{f'(x)} d\left(\sqrt{f(x)}\right).
$$

Applying integration by parts and assuming that $f'(x) \neq 0$ on the
interval. $$
\boxed{
I
=
\left[\frac{2\sqrt{f(x)}}{f'(x)}\right]_a^b
+
2\int_a^b
\sqrt{f(x)}
\frac{f''(x)}{(f'(x))^2}\,dx,
}
$$

## Sketch of the integral code

$$
\boxed{
I
= \int_a^b \frac{1}{\sqrt{f(x)}} dx
=
\left[\frac{2\sqrt{f(x)}}{f'(x)}\right]_a^b
+
2\int_a^b
\sqrt{f(x)}
\frac{f''(x)}{(f'(x))^2}\,dx,
}
$$

see file [one_over_sqrt_integral.py](./code/one_over_sqrt_integral.py)
where $f(x) = E - u(x)$

``` python
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

from p3_simpson_integration import mysimpson

def u(x):
    return x*x

def f_with_E(x, En):
    return (En - u(x))

En= 4
f = lambda x:  f_with_E(x, En)
sqrtF = lambda x:  1/np.sqrt(f(x))

a= -np.sqrt(En)  # technically we need to run root solver
b=  np.sqrt(En)

print(f"we expect integral to be pi, i.e. {np.pi}")

resLib = quad(sqrtF, a,  b)
print(f"scipy library gives {resLib=}")

Np = 10001
delta=1e-10
resMy  =  mysimpson(sqrtF, a+delta, b-delta, Np )
print(f"with small limit change Simpson gives {resMy=} for {Np=} and {delta=}")

"""
Now we will do with variable change
note that full integral is 2 * Integral from 0 to b (due to symmetry)
we split this integral 0 to b/2 where there is nothing bad with integrand
"""
int1 = mysimpson(sqrtF, 0, b/2, Np)

# and b/2 to b, where we will do substitution of a variable

def dfdx(f, x, h=1e-6):  return (f(x+h) - f(x-h))/(2*h)
def df2dx2(f, x, h=1e-6): return (dfdx(f, x+h, h) - dfdx(f, x-h, h)) /(2*h) # brute force

def g(x):
    return np.sqrt(f(x))/ (dfdx(f, x)*dfdx(f, x)) * df2dx2(f, x)


int2 = 2*np.sqrt(f(b))/dfdx(f, b) - 2*np.sqrt(f(b/2))/dfdx(f, b/2) + 2 *mysimpson(g, b/2, b, Np)

resByParts = 2 * (int1 + int2)
print(f"with a variable change trick Simpson gives {resByParts=} for {Np=}")




```