SciPy - 2 (scipy.integrate module)

SciPy has a range of different functions to integrate equations and data. Integration is a crucial tool in math and science, as differentiation and integration are the two key components of calculus. Given a curve from a function or a dataset, we can calculate the area below it. In the traditional classroom setting we would integrate a function analytically, but data in the research setting is rarely given in this form, and we need to approximate its definite integral.

Scientifc Python provides a number of integration routines. A general purpose tool to solve integrals
I of the kind


 

is provided by the quad() function of the scipy.integrate module. Numerical evaluation of a function of I type is called numerical quadrature, or simply quadature.

It takes as input arguments the function f(x) to be integrated (the \integrand"), and the lower and upper limits a and b. It returns two values (in a tuple): the first one is the computed results and the second one is an estimation of the numerical error of that result. SciPy provides a series of functions for different kind of quadrature, for example the quad, dblquad and tplquad for single, double and triple integrals, respectively.

The following program shows the integration:

import math
from math import cos , exp , pi
from scipy . integrate import quad

# function we want to integrate
def f(x):
    return exp( math .cos (-2 * x * pi )) + 3.2

# call quad to integrate f from -2 to 2
res , err = quad (f, -2, 2)

print ("The numerical result is {:f} (+ -{:g})". format (res , err ))

We have imported the quad() function from scipy . integrate module. The program solves the integral numerically for . The output of the program is shown below:

The numerical result is 17.864264 (+ -1.55117e-11)
------------------
(program exited with code: 0)

Press any key to continue . . .



The quad() function also takes optional parameters epsabs and epsrel to increase or decrease the accuracy of its computation. The default values are epsabs=1.5e-8 and epsrel=1.5e-8.

Suppose we are given data instead of some known equation and numerical integration is needed.  The figure illustrates what type of data sample can be used to approximate acceptable indefinite integrals. The original function is the line and the randomly sampled data points are shown as dots.



The following program performs the integration:

import numpy as np
from scipy.integrate import quad, trapz

# Setting up fake data
x = np.sort(np.random.randn(150) * 4 + 4).clip(0,5)
func = lambda x: np.sin(x) * np.cos(x ** 2) + 1
y = func(x)
# Integrating function with upper and lower
# limits of 0 and 5, respectively
fsolution = quad(func, 0, 5)
dsolution = trapz(y, x=x)

print('fsolution = ' + str(fsolution[0]))
print('dsolution = ' + str(dsolution))
print('The difference is ' + str(np.abs(fsolution[0] - dsolution)))

The output is as follows:

fsolution = 5.100345067540932
dsolution = 5.1336794982929606
The difference is 0.03333443075202869
------------------
(program exited with code: 0)

Press any key to continue . . . 

The quad integrator can only work with a callable function, whereas trapz is a numerical integrator that utilizes data points.

The next program shows the basic usage of the quad function:

import numpy as np
import math
import scipy.special
from scipy.integrate import quad, dblquad, tplquad


# define a simple function for the integrand
def f(x):
  
    return x

x_lower = 0 # the lower limit of x
x_upper = 1 # the upper limit of x

val, abserr = quad(f, x_lower, x_upper)
print ("integral value =", val, ", absolute error =", abserr)
print('\n')

print("If we need to pass extra arguments to integrand function we can use the args keyword argument\n")

def integrand(x, n):
    """
    Bessel function of first kind and order n.
    """
    return scipy.special.jn(n, x)
  
x_lower = 0 # the lower limit of x
x_upper = 10 # the upper limit of x
val, abserr = quad(integrand, x_lower, x_upper, args=(3,))
print (val, abserr)

#we can use a lambda function (name-less function) instead of explicitly defining a function for the integrand

print("\nUsing `inf' or `-inf' as integral limits\n")

# positive infinity
p_inf = float("inf")
# negative infinity
n_inf = float("-inf")

val, abserr = quad(lambda x: math.exp(-x ** 2), n_inf , p_inf)
print ("numerical =", val, abserr)
analytical = math.sqrt(math.pi)
print ("analytical =", analytical)
print('\n')

print("Passing lambda functions for the limits for the y integration\n")

def integrand(x, y):
    return math.exp(-x**2-y**2)
x_lower = 0
x_upper = 10
y_lower = 0
y_upper = 10
val, abserr = dblquad(integrand, x_lower, x_upper, lambda x : y_lower, lambda x: y_upper)
print (val, abserr)

The output of the program is shown below:

integral value = 0.5 , absolute error = 5.551115123125783e-15


If we need to pass extra arguments to integrand function we can use the args key
word argument

0.7366751370811073 9.389126882496413e-13

Using `inf' or `-inf' as integral limits

numerical = 1.7724538509055159 1.4202636756659625e-08
analytical = 1.7724538509055159


Passing lambda functions for the limits for the y integration

0.7853981633974476 1.3753098510218537e-08


------------------
(program exited with code: 0)

Press any key to continue . . .

In order to understand the above program a basic understanding of Bessel functions, lambda function is required.

Here I am ending this post regarding SciPy integration. Next post will see how to solve an ordinary differential equation. Till we meet next keep learning and practicing Python as Python is easy to learn!