Numerical root finding and integration gotchas

by Eugeniy E. Mikhailov and Greg Bentsen

Root finding algorithm gotchas

Diverging Newton-Raphson

Example of diverging Newton-Raphson case
Example of diverging Newton-Raphson case

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)\tanh(2.8-x)

Plot of \tanh(2.8-x)
Plot of tanh(2.8−x)\tanh(2.8-x)

Using root_scalar from scipy.optimize package

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

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)

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)\tanh(2.8-x) and bracketing algorithms

Plot of \tanh(2.8-x)
Plot of tanh(2.8−x)\tanh(2.8-x)
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

Bracketing algorithm

  • usually slower
  • robust and always converge
  • require setting the bracket
    • need a priori knowledge
  • bisection
  • false position
  • Ridders
  • Brent’s method

Non bracketing algorithm

  • Usually fast
  • Can diverge
    • sensitive to starting point
  • Newton-Raphson
  • secant

Integration gotchas

Ends of integration have infinity

One or more ends of integral have infinity

∫a∞f(x)dx \int_a^\infty f(x) dx

Solution: do a variable change.

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

∫a∞f(x)dx→∫01f(a+t1−t)dt(1−t)2 \int_a^\infty f(x) dx \to \int_0^1 f(a + \frac{t}{1-t}) \frac{dt}{(1-t)^2}

Undersampling

Peaky function
Peaky function

Not enough points

  • function has narrow features
  • rapid oscillation

Solution

  • more points are needed

Function have discontinuity

Step function
Step function

Simpson and similar methods, which assumes smooth function, will have larger error than seemingly inferior methods

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

We know that integral exist but numerically it produce infinities, typically because we divide by zero

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

We can

  • skip the small vicinity around singularity ∫δ11xdx \int_\delta^1 \frac{1}{\sqrt{x}} dx

  • add tiny offset to denominator ∫011δ+xdx \int_0^1 \frac{1}{\delta+\sqrt{x}} dx

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

there are so called adaptive methods, which adjust number of point around singularities

  • for example quad from scipy.integrate
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=∫ab1f(x)dxI = \int_a^b \frac{1}{\sqrt{f(x)}} dx

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

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

Applying integration by parts and assuming that f′(x)≠0f'(x) \neq 0 on the interval. I=[2f(x)f′(x)]ab+2∫abf(x)f″(x)(f′(x))2dx, \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

I=∫ab1f(x)dx=[2f(x)f′(x)]ab+2∫abf(x)f″(x)(f′(x))2dx, \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 where f(x)=E−u(x)f(x) = E - u(x)

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=}")