Some functions cannot be integrated analytically (on paper).

The integrals above cannot be solved analytically, that is, with pencil and paper or by hand. They can't be expressed exactly as some combination of elementary functions.

Another way to say it is that for each of those integrals above (and many others), there is no combination of trigonometric, polynomial, root, exponential or logarithmic functions that can be found which, when differentiated, give back the integrand.

Those examples are just a few among a large number of integrals we can't solve exactly, so it's worth figuring out another way to do it.

1. Riemann integration

You've seen Riemann integration as you developed the concept of the definite integral. It's a conceptually-easy method of finding the area under a curve when we run out of other options.

The scheme of Riemann integration is shown on the right. The yellow boxes are the left Riemann sum, so named because the height of each box is determined by where its top left corner intersects the function. By the same logic, the green boxes (behind) represent the right Riemann sum.

The area under the curve is broken into n regions (n = 4 here) of equal width, Δx = (b-a)/n, where   a   and   b   are the limits of the integral. Here are the summation expressions for L_N and R_N, the left and right Riemann sums.

Example — Riemann integration

Let's use Riemann integration on a definite integral we know how to calculate:

$$f(x) = -(x - 2)^2 + 4$$

That function crosses the x-axis at x=0 and x=4, so if we find the area under the curve, we get

Implementation – some Python code

Here is a short bit of Python code that performs a left Riemann sum on a function that you can define in the function func(). Paste it into your favorite IDE (integrated developing environment) or just run it from the command line. Play with the number of rectangles to see convergence toward a more precise value of the integral. This particular block of code is set up to find the area under the cosine function, but you could change it to any other function.

If you're not familiar with python coding, just skip this example.

2. Trapezoidal integration

For reasons that are obvious when you compare the figure on the right with Riemann integration, we would expect that fitting trapezoids under a curve would provide a better estimate of the area beneath it than rectangles. And like Riemann sums, the more trapezoids the better (but the gain in accuracy of both methods diminishes as the number of rectangles or trapezoids grows).

To refresh your memory, the area of a trapezoid is the width times the average of the lengths of the two parallel sides. The lengths of those sides can be calculated using the function: f(x_o), f(x₁) for the left-most trapezoid, and so on. The area of the left-most trapezoid, where (b-a)/n is the constant width of each, is

$$\int_a^b f(x) \, dx \approx \frac{b - a}{n} \sum_{i = 0}^{n - 1} \, [f(x_i) + f(x_{i + 1})]$$

If we add up the areas of all n trapezoids, we get

$$= \frac{b - a}{2n}[f(x_0) + 2f(x_1) + \dots + 2 f(x_n) + f(x_{n+1})]$$

Notice that apart from f(x_o) and f(x_n) every other term appears twice, so we can simplify a bit:

$$\int_a^b f(x) \, dx \approx \frac{b - a}{n} \left[ \frac{f(x_0) + f(x_1)}{2} + \frac{f(x_1) + f(x_2)}{2} + \dots + \frac{f(x_{n - 1}) + f(x_n)}{2} \right]$$

Then in summation form, we define trapezoidal integration of a function divided into n equal-width trapezoids with the trapezoid rule:

$$= \frac{b - a}{n} \left[ \frac{f(x_0) + f(x_n)}{2} \right] + f(x_1) + f(x_2) + \dots + f(x_{n - 1})$$

Example – Trapezoidal integration

Let's apply trapezoidal integration to the example function we integrated above, just to compare the two methods.

The solution using sixteen trapezoids is 10.63 square units, identical to our Riemann sum with 16 divisions. This brings up a couple of interesting points.

First, The integrand in this example is a symmetric function, so the right and left Riemann sums will be identical. Therefore the right and left sums, and their average are all equal.

Error analysis

Here is a plot of the errors in several trapezoid-rule integrations of our example function above (error as a percent of the known value of 32/3). Integrations using 10, 20, 40, 80 and 160 trapezoids were performed using a spreadsheet program.

Notice that while the error becomes smaller as the number of trapezoids increases, there is a diminishing return. The difference in improvement in the percent error decreases from one trial to another. Adding more iterations has makes less of a difference as we go on.

Nevertheless, modern computers are blindingly fast at doing integrations using even 1000 trapezoids, so adequate accuracy can usually be obtained.

3. Monte Carlo integration

The city of Monte Carlo in Monaco, famous in part for its central gambling casino, was founded on the laws of probability, thus the name for this probability-based integration technique.

Some integrals are not amenable to Riemann-type integration, including trapezoidal integration, so we need another way. Monte Carlo integration is one such way.

It works like this: The area to be integrated is surrounded by a bounding box, leaving about the same area above the curve as below. Because the bounding box is a rectangle, its area is known. The yellow bounding box shown here has area (b-a) x (d-c).

Pairs of random numbers, one between a and b, the other between c and d, are generated using a high-quality random number generator (it's not as easy as you might think to generate truly random numbers).

Implementation – some Python code

Here is a short bit of Python code that performs a Monte Carlo integration of the function f(x) = x⁴ - 5x³ - 9x² + 45x. This is a function for which we can obviously find the analytical integral – just a well-behaved polynomial function – but that allows us to examine the accuracy of the approximation. If you run it, you'll find that it's pretty good with a million "shots" or so.

If you're not familiar with python coding, just skip this example.