Numerical Integration

Introduction

This week, we will learn how to perform the numerical integration techniques you have been taught in PX1120. We will start with a simple example and write a code to do the integration. Finally, we will look at the various tools for doing numerical integration that exist within NumPy and SciPy packages.

Useful Links

Numpy trapz - Integration using the trapezium rule

SciPy Integration Routines - We will be integrating functions given fixed samples for now

Syntax Summary

Array Creation

Integration

The trapz() and polyint() functions are available in NumPy. More complex integration methods are coded in the scipy.integrate package. Note that this also has a version of trapz(), which appears to be the same as what's in NumPy.

Importing Libraries

from scipy import integrate

Worksheet

Numerically integrating a function at a given set of points

In this section, we will explore a number of methods for integrating a function. We will work with a polynomial, because this can be evaluated exactly (either analytically or using the poly1d package in NumPy). This will allow us to check that we are getting it right. The polynomial we will consider is:

p(x) = -10x⁴ + 4x³ + 30x² + 6x

Throughout, we will focus on integrating between 0 and 2.

As you go through the worksheet, write a program that performs all of the integration techniques introduced below on the polynomial p(x).

Polynomial Integration

We have made use of the poly1d functions in previous weeks, but not specifically the polyint() function. It works in a similar way to polyder() that we have used. For the example, you need to do

p = poly1d([-10, 4, 30, 6, 0])
q = polyint(p)
print(q)

Note that the constant of integration defaults to zero, although you can set the value if you wish. Then, to evaluate the integral, simply calculate q(x) at x=2 and x=0

integral = q(2) - q(0) 

This is the value that we will compare the results of our numerical simulation against.

Rectangular Integration

Let us start out with the simplest numerical integration technique -- rectangular integration. We will do both left and right sided rectangular integration. Although we will begin by only using a small number of points, we will make use of arrays so that we can generalize easily to more general situations.

First, define an array of points in the range between x=0 and x=2 and the corresponding values of the polynomial as

#create 3 x values between 0 and 2
x = linspace(0,2,3)
#Calculate the corresponding y values
y = p(x)

For rectangular integration, we must evaluate the function at either the left or right hand side of each interval and then multiply by the width of the intervals, which we denote dx. This is done as follows:

#Special case - 2 rectangles
#dx is equalt to the difference of 2 consecutive values in the x array
dx = x[1] - x[0]

#Approxiamte area as sum of 2 rectangles: dx*y[0]  and dx*y[1]
rect_left = dx * (y[0] + y[1])

#or
rect_right = dx * (y[1] + y[2])

Figure 1:The area under the blue line (polynomial p, between x=0 and x=2) can be approximated (not very accurately) by the area of two rectangles A and B, of width δx=1. For the first rectangle A (x between 0 and 1) there are two way to set the height, either setting it to p(x=0) as shown in the left graph, or setting it to p(x=1) as shown in the graph on the right.

p(0) vs y[0]:

This tends to cause a bit of confusion.

• p() is a function (a polynomial) and like all functions (eg sin() or cos()) it is followed by parentheses. The input to p() can be any number, eg. p(0), p(1.3), p(-3.4).
• y is an array and like all arrays you can access any element of this array using square brackets. So for example y[0] will give you the first number in y, y[1] will give you the second number in y or, as we learned in Week 3, y[-1] will return the last number in array y. Since the number in the square brackets denotes the position in the array, this number can only be an integer

The above code example is for the case of approximating an integral with 2 rectangles. Hwoever, this can be extended quite easily to a greater number of points. The easiest way would be to simply add up all the points except the last (or first) for left (or right) rectangular integration. We will instead introduce a slightly more general method that can then be easily extended to other integration methods:

1. Calculate dx
2. Generate a "weighting" array, w, of the same size as y that tells us how much of each y point we need. For left rectangular integration, this array will be [1,1,0]; for right rectangular integration [0,1,1].
3. Calculate the approximate value of the integral as dx Σ_i y_i W_i

Note: Here, we are assuming that the intervals are of a constant width (i.e. dx is the same for each interval). It's not hard to generalize to the situation where this isn't true by introducing dxi, but we will not cover that in this class.

In code, this is

#General Method
dx = x[1] - x[0]
w_left = ones_like(y)
w_left[-1] = 0
rect_left2 = dx * sum(w_left * y)

The function ones_like(y) generates an array of ones, of the same size as the array y. The code for calculating the right rectangular integration is similar, although the weighting array now starts (rather than ends) with a 0. Check that this gives the same answer as you got before.

It might seem like overkill to introduce the weighting array, but this method generalizes nicely to trapezium integration and Simpson's rule.

Increase the number of points

It is straightforward to recalculate the integral using a different number of points. All that you need to change is the definition of the points x:

x = linspace(0,2,n_points)

Everything else stays the same. Increase the number of points to 101 (i.e. 100 rectangles) and calculate the integral.

Trapezium Rule

In the trapezium rule, the average value of the function at the two edges of the interval are used in calculating the result. It is straightforward to calculate the trapezium rule numerically. The way we have done it just requires a different weighting function: w_trap = [1/2,1,1,1, ... ,1/2]. Add this to your code and use it to evaluate the integral with 101 points in x (between 0 and 2). How does the accuracy compare with the rectangle rule?

Simpson's Rule

Again, it is straightforward to generalize to Simpson's rule, by introducing another weighting function. The weighting function should be w_simp = 1/3 [1,4,2,4,2,4,...,4,1]. This one is a bit harder to generate, but is possible using some of the array manipulation techniques that were introduced in weeks 2 and 3. In particular, recall that you can set the value of the odd entries of an array using:

start = 1
stop = -1
step = 2
w[start:stop:step] = 4

Generate the Simpson's rule weighting function and calculate the integral using 101 points. How does the accuracy compare with the rectangle and trapezoidal rules?

Numpy/Scipy Integration Routines

The code that we have written above is capable of integrating a function that we have evaluated at a fixed number of (evenly spaced) points. There are routines in the NumPy and SciPy packages to perform similar calculations.

Numpy has a built in trapezium rule integrator called trapz(). Read the documentation and use this function to integrate the polynomial given above, and show that it gives the same answer as your trapezium rule integration.

Importing Modules from a Library

SciPy has additional numerical integrators. Since the functions in SciPy are not available by default, we need to first load the scipy.integrate module into python. There are several methods of doing this. The one we suggest is

from scipy import integrate

This will make all of the functions in scipy.integrate available. In order to run them, you will need to type (for example)

integrate.simps()

Alternatively, you can import the integrate module using

import scipy.integrate
scipy.integrate.simps()

Finally, you can import all of the functions in a module using

from scipy.integrate import *
simps()

This last method seems easiest, but it can lead to overwriting existing functions. For example, there is a function trapz in numpy. By running the code above, you will have overwritten the NumPy trapz() function with the SciPy version. In this case, it doesn't really matter as they do essentially the same thing. However, there will be cases where you overwrite a function with another function of the same name that performs a different operation.

Using Integration methods from SciPy

You should now have access to the simps() integration method. Using the help message in the object explorer, use this function to integrate the polynomial given above, and show that it gives the same answer as your Simpson's rule integration.

SciPy also provides a function cumtrapz() that cumulatively integrates the function for you. Rather than returning a single answer, it returns an array (that is 1 entry shorter than what you pass in) containing the value of the integral up to the corresponding x-value.

Note: By default, cumtrapz() will return an array that is 1 entry shorter than the arrays you pass into it. So, if you want to make a plot involving a function that you have integrated with cumtrapz, you will need to drop the first point of the x-array to ensure it is the same length as the array returned by cumtrapz.

Subplots

Sometimes it is useful to compare several similar data sets, or display one data set in different ways. Multiple plots can be displayed on one figure window using subplots. The subplot() function requires 3 numbers, the first two give number of subplots in the figure (number of rows and number of columns) and the third gives which subplot you wish to be active. For example:

subplot(2,2,1)

This will create a figure with 4 subplots (2 rows x 2 columns) and make the first one active. Note that numbering of subplots increases left to right then top to bottom, like reading a book. Here is an example code to create 3 subplots (3 rows, 1 column) showing powers of x.

subplot(3,1,1)
plot(x, x**2)
ylabel("x^2")

subplot(3,1,2)
plot(x, x**3)
ylabel("x^3")

subplot(3,1,3)
plot(x, x**4)
ylabel("x^4")
xlabel("x")

Excercises