Data Science and Computing with Python for Pilots and Flight Test Engineers → Numerical Integration

Numerical Integration

In this section, we shall write ourselves a simple numerical integrator. Integrators are useful for two things:

to computer the area under a curve (i.e. integrating a function)
to solve a differential equation numerically.

While professionally coded integrators exist in libraries, it is useful to have your own simple integrator. It may not be as fast, accurate, and adaptive as those you will find in some libraries, but if a problem with an integration exists (e.g. a wildly fluctuating function), your own integrator may allow you to analyze the situation and debug the problem better, because it is a simple code that you can modify at will. The integrators from libraries are like a black box: they either work or they do not. While some have adjustable parameters, if the integration still does not work it is a challenge to enter the code and analyse what happens on each line individually and how the integration fails.

Integration of a Function (Computing the Area under a Curve)

In this section you will write a simple numerical integrator based on a Riemann sum to compute the area under a curve. You can use the professional integrator at the end of this section to check whether your own integrator works correctly.

Notation: The Summation Symbol

The summation symbol is shorthand for adding multiple terms together which bear an index. Let $a_1$, $a_2$, $a_3$, $a_4$, $a_5$ be five numbers.

Then we can write the following sum compactly using a summation symbol:

$$ a_1 + a_2 + a_3 + a_4 + a_5 = \sum_{i=1}^{5} a_i $$

In the right-hand expression, $i$ is called an index and it is understood to run from $1$ to $5$, as indicated below and above the summation symbol. The summation symbol is the capital letter Sigma of the Greek alphabet.

The Riemann Sum

Without showing the proof, we claim that the area $A$ under a curve $f(x)$ can be approximated by a Riemann sum:

$$ A \approx \sum_{i=0}^{n-1} f(x_i) \Delta x_i $$

where $x_i$ are discrete points along the horizontal axis if function $f(x)$ (enumerated by index $i$), $f(x_i)$ is the value of the function $f(x)$ at horizontal axis point $x_i$, and $\Delta x_i:=x_{i+1}-x_i$ is the difference (distance) between two neighboring points $x_i$ and $x_{i+1}$ along the horizontal axis.

The more points $x_i$ we choose, the smaller $\Delta x_i$ becomes and the more accurately the above Riemann sum approximates the actual area $A$ under the curve $f(x)$. In the limit of infinitesimally small $\Delta x_i \rightarrow 0$ the approximation becomes exact, resulting in the so called Riemann integral denoted by

$$ A = \int f(x) \, dx = \lim_{\Delta x_i\rightarrow 0\\n\rightarrow\infty} \sum_{i=0}^{n-1} f(x_i) \Delta x_i $$

A further accuracy improvements (for finite $\Delta x_i$) can be achieved in the above formula if we replace $f(x_i)$ by $(f(x_i) + f(x_{i+1}))/2$.

Writing Your Own Numerical Integrator

Using the formula for the Riemann sum and a while loop, you can now write your own simple numerical integrator, which will compute the area under the curve of a function. To help you get started, we have included some starter code below, and you just need to fill in the two lines in the while loop correctly.

In the next section further below we will show you how to use a professional numerical integrator. You can use the result from that section to compare it to your result here to see if you got the correct answer.

import numpy as np
import matplotlib.pyplot as plt

def function_integrator(x_start, x_end, dx, function, params=[]):
    """ A simple function integrator. Integrates a function f(x) from x_start to x_end. """
    # The argument params contains the parameters of the integrated function, which are held constant during integration.
    # If the function has no arguments, you can omit supplying them.
    
    # Initialize area under the curve to zero at start of integration. 
    area = 0.0
    x = x_start
    # Integrate:
    while (x<x_end):
        area = area + function(x, params) * dx
        x = x + dx   

    return area

This is the function we wish to integrate in this particular example (modify this to integrate other functions as you need):

def example_function(x, params=[]):
    
    y = params[0] * x**3 + params[1] * x**2 + params[2] * x + params[3] # change this as you like.
    
    return y

# Run your own numerical integrator:

# Starting and end point of integration:
x_start = 0.0
x_end = 2.0
# Time step chosen in Riemann sum (make this smaller for more accuracy):
dx = 0.0001

# Parameters of function to be integrated (if any):
params = [0, 1, 0, 0]

# Call our simple numerical integrator:
area = function_integrator(x_start, x_end, dx, example_function, params)

# Print resulting area under the curve:
print("Area under the curve:", area)

Copy and paste the above code cells into your Jupyter notebook and run your own simple numerical integrator.

Professional Numerical Integration

In real life as a programmer, scientist, or engineer, you would likely not write your own integrator as above and would instead use a professionally coded optimized numerical integrator instead. A variety of numerical integrators exist, and we show you how to use one below.

To use a professional integrator, we need to import it first into the notebook from a module. The module scipy, which we have encountered before, contains several integrators, and we shall import the integrator “quad”. After doing so, as you can see, the integration of the fit function (computation of the area under the curve) can be accomplished in just one line of code.

import numpy as np
import matplotlib.pyplot as plt

# Import the quad numerical integrator from the scipy module:
from scipy.integrate import quad

Define some arbitrary function of one variable $x$, which you wish to integrate. (Note that parameters $a, b, c, d$ are constants; only $x$ will be integrated over here.)

def example_function(x, a, b, c, d):
    
    y = a * x**3 + b * x**2 + c * x + d # change this as you like.
    
    return y

Now integrate this function with the quad integrator from the scipy module:

# Using the quad integrator from the scipy module:

# Function parameters:
a = 0
b = 1
c = 0
d = 0

x_start = 0.0
x_end = 2.0

# Integrate function "function2" with parameters a and b from start time to end time of data recording:
result = quad(example_function, x_start, x_end, args=(a,b,c,d))

# Extract information provided by the numerical integrator:
area = result[0]
numerical_error = result[1]

# Print result on screen:
print ("Area under curve = ", area, "computed with accuracy", numerical_error)

Copy and paste the above code cells into your Jupyter notebook and compare the accuracy of your own integrator from the previous section with the result from the quad integrator above. Also compare the run time, and see which integrator is faster. You can determine the runtime by paying attention to the time from the moment you hit Shift + Return to when the result appears on the screen.

You will likely see that quad probably works better than your own integrator. Nonetheless, writing your own integrator was useful. Not only did it reinforce the concept of integration (computing the area under the curve) using a Riemann sum for you, in some cases it is actually useful to use your own integrator. While it may not be the fastest and/or most accurate, it can help you analyze problems which may occur during the integration (e.g. for a wildly fluctuating function), because you have full control over the code and are not using a black box which either works or does not, with little that you can do about it.