Numerical Integration

Here is an introduction to Numerical Integration with Python!

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

Let’s look at a simple function:

\[g(x) = \sin(x)\]

which we know the analytical integral as:

\[G(x) = - \cos(x)\]

and let’s integrate that from 0 to \(\frac{5\pi}{4}\). Analytically, this is:

\[\int_0^{\frac{5\pi}{4}} g(x) \, dx = - \cos(x) |_0^{\frac{5\pi}{4}} = - \cos(\frac{5\pi}{4}) + \cos(0) = \frac{\sqrt{2}}{2} + 1 = 1.7071\]
np.sqrt(2)/2 + 1

Out:

1.7071067811865475

We can visualize this as the area under the curve - which is all a numerical integral is; an approximation of the area under the curve. To do a numerical integration, we simply discretize our x space and then evaluate our function. Once we have that, we can sum the evaluation and multiply by the discretization factor to yield the effective area under the curve.

So let’s discretize \(g(x)\) into several rectangles which we can use to solve for the area under this curve.

g = lambda x: np.sin(x)
a, b = 0, 5*np.pi/4

y = np.linspace(a, b)

dy = np.pi/16
yy = np.arange(a, b, dy)

def plot_rects(xx, f, delta):
    for v in xx:
        value = f(v)
        if value >= 0.0:
            loc = 0.0
            h = value
        else:
            loc = value
            h = 0 - value
        r = mpl.patches.Rectangle((v-delta/2., loc),
                                  width=delta, height=h,
                                  edgecolor="k", linewidth=1)
        plt.gca().add_patch(r)
    return

def plot_g():
    plt.plot(y, g(y))
    plt.plot(yy, g(yy), "ro")

    plot_rects(yy, g, dy)

    plt.xlabel("$x$")
    plt.ylabel("$g(x)$")

plot_g()
plt.show()
../_images/sphx_glr_numerical-integration_001.png

Out:

/home/travis/build/banesullivan/mywebsite/python-blog/numerical-integration.py:77: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

So now we can sum the areas of all of those rectangles to get an approximation of the area under the curve (the integral)

np.sum(g(yy-(dy/2)) * dy)


# That's not too bad considering our rectangles do not perfectly follow the
# curve. Let's now try choosing a :math:`\Delta x` value that will increase the
# sampling creating finer spaced rectangles to improve our approximation.

dy = np.pi/32
yy = np.arange(a, b, dy)

plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plot_g()
plt.title("Area = {:.3f}".format(np.sum(g(yy) * dy)))

plt.subplot(1,2,2)

dy = np.pi/128
yy = np.arange(a, b, dy)

plot_g()
plt.title("Area = {:.3f}".format(np.sum(g(yy) * dy)))

plt.tight_layout()
plt.show()
../_images/sphx_glr_numerical-integration_002.png

Out:

/home/travis/build/banesullivan/mywebsite/python-blog/numerical-integration.py:106: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()
np.sum(g(yy) * dy)

Out:

1.7156985903219493

Now that’s much better! as we decrease the \(\Delta x\) value, we see the apprximationm start to converge on the analtical answer for the integral. The rectangles appear to be consistently overestimating the integral.

Integrating a Complex Function

Now let’s apply that concept to numerically integrate a more complex function such as:

\[f(x) = \sin(\frac{1}{x})\]

from 0 to \(2\pi\). This allows us to really showcase where numerical integration benefits as I don’t even know where to begin when it comes to analytically integrating that.

f = lambda x: np.sin(1/x)

To do so, we’ll need to define a range of x values on which to evaluate the function. Let’s try out a few different sets of x values to plot the values of \(f(x)\) to gain insight into how this function behaves.

First, let’s see what happens when we plot this function at linear intervals between 0 and \(2\pi\) (note that we cannont have a division by zero, so our lower boound must be some small value, not 0: we’ll use \(10^{-9}\).

num = 100000
zero = 1e-9

# Create an array of x values to evaluate
x = np.linspace(zero, 2*np.pi, num)

# Evaluate and plot the values
plt.plot(x, f(x))
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
plt.show()
../_images/sphx_glr_numerical-integration_003.png

Out:

/home/travis/build/banesullivan/mywebsite/python-blog/numerical-integration.py:153: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

Huh, that’s weird. The function behaves quite different from a typical \(\sin\) function and we appear to see some sort of decay - let’s look into that a bit more. At what x value are we seeing a change from oscillatory behavior to a decay?

From visual inspection, it looks like some value above \(\frac{1}{\pi}\) (0.318). Let’s reevaluate the function from that x value onward. Then find the maximum and see what interval of \(\pi\) is causing the change in behavior.

subset = np.linspace(0.3, 2*np.pi, num)
subset[np.argmax(f(subset))] * np.pi

Out:

1.9999924576268395

Aha! It appears that the x value of \(\frac{2}{\pi}\) marks the start of decay. This makes sense because \(\sin(1/2/\pi) = \sin(\pi/2) = 1\)

plt.plot(x, f(x), label="$f(x)$")
plt.plot(2./np.pi, f(2./np.pi), 'ro', label=r"$x=\frac{2}{\pi}$")
plt.legend()
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
plt.show()
../_images/sphx_glr_numerical-integration_004.png

Out:

/home/travis/build/banesullivan/mywebsite/python-blog/numerical-integration.py:178: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

So now we know that have a speration point in our function - this might come in handy later. as we may have to treat those two regions of the function separately when numerically integrating.

Now, let’s explore the jumbled up portion of the plot for very low x values. There appears to be some high frequency oscillation between -1 and 1, but we can’t really tell. To do this, let’s create an exponentially varying x space to evaluate our function which will bring out finer sampling for lower values and see how the function behaves as x approaches zero.

x = np.geomspace(1e-2, 2*np.pi, num)

plt.plot(x, f(x))
plt.gca().set_xscale('log')
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
plt.show()
../_images/sphx_glr_numerical-integration_005.png

Out:

/home/travis/build/banesullivan/mywebsite/python-blog/numerical-integration.py:197: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

Now that’s cool! \(f(x)\) appears to have increasing frequency as x approaches zero. This is quite different from a normal \(\sin\) function where we don’t see any frequency change at all. This frequency variation will definitely pose a challenge for our numerical integration as we will need to make sure we properly sample this signal without aliasing.

Let’s see if we can come up with an equation for capturing that frequency variation as something we can use to build out our values that we will use for the numerical integration.

First, let’s evaluate our function between integer changes in \(\pi\) to see how \(f(x)\) oscilates between -1 and 1. We know from above that the last value of 1 occurred at \(\frac{2}{\pi}\), so let’s use decreasing values from there.

f(3./np.pi)

Out:

0.8660254037844386
f(3./(2.*np.pi))

Out:

0.8660254037844387
f(2./(3*np.pi))

Out:

-1.0
f(2./(4*np.pi))

Out:

-2.4492935982947064e-16
f(2./(5*np.pi))

Out:

1.0

So this is interesting. I’m seeing a weird pattern where \(\frac{2}{n \pi}\) appears to oscillate \(f(x)\) between -1 and 1. Let’s test this out.

xx = lambda n: 2/(n*np.pi)
nn = 65
n = np.arange(1, nn)
xn = xx(n)

plt.plot(xn, f(xn), "ro")
plt.plot(x, f(x))
plt.gca().set_xscale('log')
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
plt.show()
../_images/sphx_glr_numerical-integration_006.png

Out:

/home/travis/build/banesullivan/mywebsite/python-blog/numerical-integration.py:247: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

It appears that the equation \(\frac{2}{n \pi}\) captures the oscilatory behavior of \(f(x)\) quite well! We see that it marks the minimums, maximums, and zero-cross points for every period of \(f(x)\). This is just what we need to effectively integrate \(f(x)\).

To start integrating this function, we simply need to sample several times between values of \(n\). Doing so is as simple as bumping up the number of setting a step size in numpy’s arange method. Let’s use a step size of 0.25 as that will give us 4 values between each cross-over and the min/max point. Let’s see if that gives us a visually pleasing approximation of our function.

def plot_f():
    plt.gca().set_xscale('log')
    plt.plot(x, f(x))
    plt.plot(xn, f(xn), "ro")

    for i, v in enumerate(xn):
        if i == len(xn)-1:
            continue
        value = f(v)
        if value >= 0.0:
            loc = 0.0
            h = value
        else:
            loc = value
            h = 0 - value
        delta = abs(xn[i] - xn[i+1])
        r = mpl.patches.Rectangle((v - delta/2, loc),
                                  width=delta, height=h,
                                  edgecolor="k", linewidth=1)
        plt.gca().add_patch(r)
    plt.xlabel("$x$")
    plt.ylabel("$f(x)$")
    return
dx = 0.25
n = np.arange(1, nn, step=dx)
xn = xx(n)

plot_f()
plt.show()
../_images/sphx_glr_numerical-integration_007.png

Out:

/home/travis/build/banesullivan/mywebsite/python-blog/numerical-integration.py:295: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

That looks decent, but I think we can do better. Let’s try 0.1 for 10 samples between the cross overs and the min/max

dn = 0.1
# 1/np.pi**2
n = np.arange(1, nn, step=dn)
xn = xx(n)

plot_f()
plt.show()
../_images/sphx_glr_numerical-integration_008.png

Out:

/home/travis/build/banesullivan/mywebsite/python-blog/numerical-integration.py:307: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

Now that looks pretty darn good! How about we try integrating \(f(x)\) now! To do so, we simply take our values and multiply by the \(\Delta x\) value (the rectangle widths). It’s important to note here that \(\Delta x\) is not a constant this time, so well need to compute an array of all of our \(\Delta x\) values first.

dx = np.pi/(128*4)
x = np.arange(2/np.pi, 2*np.pi, step=dx)
right = np.sum(f(x) * dx)
right

Out:

2.1014656967879217

The right side there is pretty good (check with wolfram alpha)

def compute_left(dn, nn=1000):
    n = np.arange(1, nn, step=dn)
    xn = xx(n)
    dx = np.abs(np.diff(xn))
    return np.sum(f(xn)[:-1] * dx)

left = compute_left(dn)
left

Out:

0.17968499172741753

Well that’s not very good as we know the true value for the left side to be 0.164619 from Wolfram Alpha. Let’s see if incresing the discretization helps.

left = compute_left(0.0001)
left

Out:

0.16463442561346053
left+right

Out:

2.266100122401382

The full integral is known to be 2.26277 from Wolfram Alpha so that’s not too bad of an approximation!

And there you have it! The integral of \(\sin(\frac{1}{x})\) has been numerically integrated!

Total running time of the script: ( 0 minutes 3.727 seconds)

Gallery generated by Sphinx-Gallery