How to Compute Numerical integration in Numpy (Python)?

The definite integral over a range (a, b) can be considered as the signed area of X-Y plane along the X-axis. The formula to compute the definite integral is:

where F() is the antiderivative of f().

We can then differential the range from a to b into as many steps (rectangles) as possible and sum up the area of the rectangles.

We now need to write our numerical integration function.

Say we want to integrate f() between a and b. We are first going to evaluate f(x) at lots of positions between a and b, (say N positions). To do this we can first generate a number line with N points between a and b stored in the vector x. We can do this using numpy’s linspace

function x = np.linspace(a, b, N)

We then pass the vector x to our function f() in the usual way.

fx = f(x)

This gives us the heights of the rectangles on the previous page. There will be N rectangles between a and b so we can work out that their width will be (b−a)/N. So the area of each one is fx(i).(b−a)/N. So our integral, which is the total area of all the rectangles combined, is given by:

np.sum(fx)*(b-a)/N

The more slithers we use the better the accuracy will be. The function signature for our integrate function should look like

integrate(f, a, b, N)

Putting this all together we have the code below:

def integrate(f, a, b, N):
    x = np.linspace(a, b, N)
    fx = f(x)
    area = np.sum(fx)*(b-a)/N
    return area

We can now use the function to integrate a sine curve between 0 and and pi/2. This should produce 1. Let’s run it using 100 steps.

integrate(np.sin, 0, np.pi/2, 100)

This gives:

0.99783321217729803

The answer is off by about 0.002. This is not bad. In fact the integrate function above is simple but it is not quite right. If you look back at the figure and consider the points generated by linspace then you might be able to spot the error. Try fixing it in rerunning the test. If you get it right the error will be about 200 times smaller (Rewrite the integrate function in the cell below)

def integrate(f, a, b, N):
    x = np.linspace(a+(b-a)/(2*N), b-(b-a)/(2*N), N)
    fx = f(x)
    area = np.sum(fx)*(b-a)/N
    return area

integrate(np.sin, 0, np.pi/2, 100)

This gives:

1.0000102809119051

In general, the more steps, the higher accuracy of the integral but at the cost of the computation time and complexity.

–EOF (The Ultimate Computing & Technology Blog) —

586 words
Last Post: How to Pass Function as Parameter in Python (Numpy)?
Next Post: How to Compute the mean of a distribution using Python and Numpy?