What They Don’t Tell You About Numerical Integration

This is what you can do if numerical integration fails. It happens more often than you may think.

In elementary courses on numerical methods, people are taught some standard ways how one can integrate functions numerically, like Simpson’s rule, Gaussian integration, etc. While those methods are not too hard to implement, in practice, you should rarely be in a position where you have to do that. There is always a solution already waiting for you. However, this leaves the impression, that for numerical integration, all you have to do is apply a function from a package. But all too often, this is not the case, and people don’t know what to do then! Basically, there are two situations, where any numerical integration scheme has problems: apparent poles in the integrand and infinite integration boundaries. I will show you, based on the nasty integral of yesterday’s post, what you can do, if you are faced with problems like that. It turns out, that this integral is not only hard for SymPy to solve analytically, but also hard to solve numerically. As a reminder, the integral that we want to solve now is

and we call the integrand 𝑓(𝑥). From yesterday’s article, we know that the analytic solution is 2𝜋/3. Let’s see how we can compute that numerically and get around the problems we will encounter.

The problem with SciPy at apparent poles

In the Python world, the first address for solving integrals numerically would be the SciPy package, which is part of the official Python scientific stack. And the general-purpose function for integration is scipy.integrate.quad , like "quadrature". A first try at the problem is straightforward. For now, let's ignore the infinite boundaries and just take big numbers hoping that this will give a good enough result (it won't!).

import numpy as np
from scipy.integrate import quad

def f(x):
    return np.sin(x)**4 / x**4

quad(f, -100, 100)

Out: (nan, nan)

SciPy complains that there is a problem with the integrand and proposes to check, if there are singularities in the integrand. Well, actually there is no singularity, because as 𝑥→0, we have 𝑓(𝑥)→1, but numerically, the 𝑥 in the denominator is a problem because numerically, we don’t have a limiting process. So when quad tries to evaluate the integrand at 𝑥=0, it gets an undefined number. SciPy proposes to split up the integration interval and treat each interval separately. We will do that, but what then? The problem remains for the subinterval around 𝑥=0! One easy and quick-and-dirty solution would be to add a tiny positive number like that:

def f_dirty(x):
    eps = 1.E-20
    return (np.sin(x)**4 + eps) / (x**4 + eps)

sol, resid = quad(f_dirty, -100, 100)

Out: (2.0943948562777077, 2.713381952505866e-08)

The quad function returns a tuple with the zeroth entry being the integral value and the second one an estimated error. So this works in principle.

But as I said, that’s quick and dirty. Another thing that is worth trying is to exploit the symmetry of the integrand and use

It seems like we have gained nothing because we still have to integrate over the numerically problematic point 𝑥=0. But this time, the problem is at the boundary of the integration domain and not inside! This makes a difference because the implementation of the integration algorithms is more robust on the boundaries. Let’s give it a try:

2 * quad(f, 0, 100)[0]

Out: 2.0943948562777077

Yes, that worked! But in other cases, it may not. So it’s worthwhile to see, what else one can do. We could also split up the interval into three parts and isolate the problematic region:

for some small positive constant 𝜀. In the interval (−𝜀,𝜀) we then replace 𝑓(𝑥) by its (truncated) Taylor expansion. You can get that by using sympy.series like that:

import sympy
x = sympy.symbols('x')
f_ = sympy.sin(x)**4 / x**4
sympy.pycode(f_.series(n=12).removeO())

Out: '-4/10395*x**10 + (62/14175)*x**8 - 34/945*x**6 + (1/5)*x**4 - 2/3*x**2 + 1'

Where I have used sympy.pycodeto produce Python code. And now use it:

def f_expansion(x):
    return -4/10395*x**10 + (62/14175)*x**8 - 34/945*x**6 + (1/5)*x**4 - 2/3*x**2 + 1

eps = 0.001
quad(f, -100, -eps)[0] + quad(f_expansion, -eps, eps)[0] + \
     quad(f, eps, 100)[0]

Out: 2.0943948562777077

So that worked, too.

Infinite integration domains

We know from yesterday’s article that the exact result should be 2𝜋/3 which is approximately 2.0943951023931953. Most of the error in our numerical integral at this stage is that we don’t integrate out to infinity but only to 100. How does the error get smaller when take higher values for the boundaries?

import matplotlib.pyplot as plt

lims = np.logspace(1, 5, 20)
vals = np.array([2*quad(f, 0, lim)[0] for lim in lims])
errs = np.abs(2*np.pi / 3 - vals)

plt.loglog(lims, errs)
plt.grid()

First, with larger integration regions, things get better, then we get stuck around an error of 1E-9. In fact, SciPy already complains with a warning about problems around a limit of 1000. So what can we do? The usual way to approach this is to map the infinite integration domain to a finite interval. One possibility is to make the substitution:

which maps all 𝑥 from minus infinity to infinity to the interval (−𝜋/2,𝜋/2).

Then we have

and

So let’s split up the integral once more:

And put it in code:

def f2(theta):
    return np.sin(np.tan(theta))**4 * np.tan(theta)**-4 \
            *(1 + np.tan(theta)**2)

a = 100
eps = 1.E-12
2 * quad(f, 0, a)[0] + 2 * quad(f2, np.arctan(a), np.pi/2 - eps)[0]

Out: 2.09439510234672

With this approach, we have finally reduced the error to 1E-11.

That’s it for today. I hope it was useful for you and thanks for reading!

If you liked this article, you may want to consider becoming a Medium member to get unlimited access to all Medium articles. If you register using this link you can even support me as a writer, as I will get a share of $2.27 (at the time of this writing) from the Medium fee.