Introduction to Numerical Integration

Introduction

Numerical integration (sometimes called numerical quadrature) is a term referring to a variety of procedures (or algorithms) for approximating the numerical value of a definite integral. There are a number of scenarios when using numerical integration on a function \(f(x)\), rather than finding the value of the integral through an anti-derivative, can be advantageous.¹  First, you may only have access to some sampled values of \(f(x)\), meaning that you do not have a functional expression for \(f(x)\).  This is very common in applications.  Second, even if you do have an expression in the variable that defines the function, it may be impossible to find an anti-derivative in terms of elementary functions.  We saw examples where this was the case when discussing arclength and surface area.  Finally, there are situations when you can compute an anti-derivative, but doing so is much more computationally costly than using numerical integration.

The Trapezoid Rule

One method used in numerical integration is the trapezoid rule.  To understand the trapezoid rule, we recall how Riemann sums approximate the area under a curve.

When \(f(x) \ge 0\) for \(x\) within the interval \([a, b]\), then the integral \(\int_a^b f(x)\ dx\) can be interpreted as the area between the graph of \(f(x)\) and the \(x\)-axis.  Riemann sums approximate that area using areas of rectangles and dividing the interval \([a,b]\) into some number of subintervals \([x_i, x_{i+1}]\), with \(0\le i\le n-1\), where \(x_0 = a\).  In most introductions to Riemann sums, each subinterval has width \(\Delta x = (b-a)/n\), and a rectangle for \([x_i, x_{i+1}]\) has width \(\Delta x\) and height equal to either \(f(x_i)\) or \(f(x_{i+1})\).

In the trapezoid rule, we replace each rectangle with a trapezoid.  The trapezoid for the interval \([x_i, x_{i+1}]\) has corners at the four points \((x_i,0)\), \((x_i, f(x_i))\), \((x_{i+1}, f(x_{i+1}))\), and \((x_{i+1}, 0)\).  In other words, the top of the trapezoid region is on the secant line for \(f(x)\) between \(x=x_i\) and \(x=x_{i+1}\).  The area of that trapezoid is

\( \Delta x\left(\frac{f(x_i) + f(x_{i+1})}{2}\right). \)

Note what happens with the sum of the areas, one over \([x_0,x_1]\) and the other over \([x_1,x_2]\): we get

\( \Delta x\left(\frac{f(x_0) + f(x_{1})}{2}\right) + \Delta x\left(\frac{f(x_1) + f(x_{2})}{2}\right), \)

which equals

\(\frac{\Delta x}{2}\left(f(x_0)+f(x_1)+f(x_1)+f(x_2)\right) = \frac{\Delta x}{2}\left(f(x_0)+2f(x_1)+f(x_2)\right). \)

The function value \(f(x_1)\) appears twice since \(x_1\) is an endpoint in two of the intervals.  In the approximation from a sum of \(n\) trapezoid areas, the same will occur for values \(f(x_i)\) with \(1\le i\le n-1\).  So, to approximate the integral of \(f(x)\) over \([a,b]\), we compute

\(\int_a^b f(x)\ dx \approx \frac{\Delta x}{2}\left(f(x_0) + f(x_n) + \sum_{i=1}^{n-1}2f(x_i)\right).\)

Implementing the Trapezoid Rule

In order for SageMath to make this computation, it needs to know the function \(f(x)\), the endpoints \(a\) and \(b\), and how many subintervals will be used — the value of \(n\).  Hence, we should be able to write a custom procedure that takes these four items as input and then makes and returns the trapezoid rule computation.  Within the procedure, we will need to determine each \(x_i = a + i\Delta x\).

Say that we want to call the procedure (or function) trap_rule.  To make the definition of trap_rule in SageMath we should have code that looks like the following.²

def trap_rule(f, a, b, n):
    # lines that determine Deltax and the list of subinterval endpoints X[i]
    # lines assigning a variable (call it: trap_sum) 
    #   the value of the sum in the trapezoid rule, times Deltax/2
    return trap_sum

Recall that \(\Delta x = (b-a)/n\).  I also hinted at how to get the list of \(x_i\) values in a paragraph above.  Each \(x_i\) is equal to \(a + i*\Delta x\). We want to get these into a list.  The Python way to do this is to put the following, once Deltax is assigned.

X = [a + i*Deltax for i in range(n+1)]

Then, we can reference \(x_i\) by typing X[i]. (An equivalent way to do this, with more lines of code, is the following for loop.)

X = []
for i in range(n+1):
    X.append(a + i*Deltax)

In SageMath and Python, say that a variable (let’s call it part_sum) has been assigned a numerical value.  To add a number to part_sum, say 3 (so that its value is 3 larger the next time you use part_sum), you can type the following.

part_sum += 3
# this adds 3 to part_sum, then reassigns that larger number to part_sum

So, to get our value of trap_sum that will be the approximation to the integral we can put the following in the trap_rule definition (after assigning Deltax and X).

part_sum = f(X[0])
for i in range(1,n):
    part_sum += 2*f(X[i]) # 2*f(X[n-1]) is the last term in the loop
part_sum += f(X[n])
trap_sum = (Deltax/2)*part_sum

Exercise 1. Using what we have discussed in this section, fill out the definition of trap_rule (to get a working procedure) using 7 lines of code that are between the first line def trap_rule(f,a,b,n): and the last line return trap_sum.  After running the code block containing the definition of trap_rule, evaluate the following two lines in a new code block (getting an approximation to the integral \(\int_0^3\frac{1}{1+x^2}\ dx\) that uses 10 trapezoids).  When you do, you should get an approximation equal to 1.24859642022472.

f(x) = 1/(1+x^2)
trap_rule(f, 0.0, 3.0, 10)

Exercise 2. With your procedure trap_rule, use the trapezoid rule with 50 trapezoids to approximate the following two integrals.

• • • \(\int_0^2 e^{(-x^2)}\ dx\)
    • \(\int_1^3 (\sin(3\ln(x)) + 0.3)\ dx\)

Implementing the Midpoint Rule

The midpoint rule approximates the integral \(\int_a^b f(x)\ dx\) with areas of rectangles over the subintervals \([x_i, x_{i+1}]\), like a right-hand or left-hand Riemann sum.  However, the heights of the rectangles do not use values of \(f(x)\) at one of the endpoints; instead we compute the midpoint \(\overline x_i = \frac{x_i+x_{i+1}}{2}\) and set the height of the rectangle over \([x_i, x_{i+1}]\) equal to \(f(\overline x_i)\).

Written out in full, for the midpoint rule we compute

\(\int_a^b f(x)\ dx \approx \Delta x\left(f(\overline x_0) + f(\overline x_1) + \ldots + f(\overline x_{n-1})\right).\)

To implement this, we can define a procedure in a similar way to how we did with trap_rule. One difference is that, after assigning the list X as we did before, we will want a list of the midpoints. In other words, after the definition of X insert the line below to make a list of midpoints M.

M = [0.5*(X[i]+X[i+1]) for i in range(n)]

In the definition of trap_rule, we made it so that part_sum was equal to \(f(x_0)+2f(x_1)+\ldots+2f(x_{n-1})+f(x_n)\), before multiplying it by \(\Delta x/2\) and returning that as the approximation.  In order to get the midpoint rule approximation, let’s make part_sum equal to \(f(\overline x_0) + f(\overline x_1) +\ldots+ f(\overline x_{n-1})\). Doing so only requires a slight adjustment to the lines that deal with part_sum:

part_sum = f(M[0])
for i in range(1,n):
    part_sum += f(M[i])

In addition to those changes, make sure that you do not add f(X[n]) to part_sum after the for loop. Also, for the midpoint approximation you multiply this part_sum by \(\Delta x\) rather than \(\Delta x/2\).

Exercise 3. Within SageMath, in a new code block, define a procedure mid_rule that computes the midpoint rule approximation.  This can be done by adapting parts of the definition of trap_rule as explained above.  After running the code block with your definition of mid_rule, evaluate the following two lines in a new code block to use the midpoint rule to approximate the integral \(\int_0^3\frac{1}{1+x^2}\ dx\) with 10 rectangles.  When you do, you should get an approximation equal to 1.24927020548728.

 
f(x) = 1/(1+x^2) 
mid_rule(f, 0.0, 3.0, 10) 

Exercise 4. With your procedure mid_rule, use the midpoint rule with 50 rectangles to approximate the following two integrals.

• • • \(\int_0^2 e^{(-x^2)}\ dx\)
    • \(\int_1^3 (\sin(3\ln(x)) + 0.3)\ dx\)