Numerical integration

Simple example of numerical integration of a one-dimensional and multi-dimensional integrals using the quad() from the integrate sub-module of scipy.

Marty Ligare, August 2020

In [2]:
import numpy as np
from scipy import integrate   # importing sub-module of scipy

import matplotlib as mpl
import matplotlib.pyplot as plt
In [3]:
# Following is an Ipython magic command that puts figures in the  notebook.
%matplotlib notebook
  
# M.L. modifications of matplotlib defaults using syntax of v.2.0 
# More info at http://matplotlib.org/2.0.0/users/deflt_style_changes.html
# Changes can also be put in matplotlibrc file, or effected using mpl.rcParams[]
mpl.style.use('classic')
plt.rc('figure', figsize = (6, 4.5))            # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True)             # Adjusts supblot parameters for new size

Single variable

In [4]:
def func1(x):      # Continuous function
    return x**3

def func2(x):      # Discontinuous function
    if x< 2:
        return x**2
    if x>2:
        return x**2
    
def func3(x):    # Function with a singluarity
    return np.sin(x-2)/(x-2)
In [5]:
value, error = integrate.quad(func1,1,2)
value, error
Out[5]:
(3.7500000000000004, 4.1633363423443377e-14)

If there are discontinuities or singularities, quad will fail, eg.,

In [6]:
value, error = integrate.quad(func2,1,3)
value, error

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-6-da1962f33e71> in <module>
----> 1 value, error = integrate.quad(func2,1,3)
      2 value, error

/usr/remote/anaconda-3.7-2020-05-28/lib/python3.7/site-packages/scipy/integrate/quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    350     if weight is None:
    351         retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
--> 352                        points)
    353     else:
    354         if points is not None:

/usr/remote/anaconda-3.7-2020-05-28/lib/python3.7/site-packages/scipy/integrate/quadpack.py in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    461     if points is None:
    462         if infbounds == 0:
--> 463             return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
    464         else:
    465             return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

TypeError: must be real number, not NoneType

Specify troublesome points:

In [7]:
value, error = integrate.quad(func2,1,3,points=[2,])
value, error
Out[7]:
(8.666666666666668, 9.621932880084691e-14)
In [8]:
value, error = integrate.quad(func3,0,4,points=[2,])
value, error
Out[8]:
(3.210825953605389, 3.5647329017567276e-14)

Multivariable

Volume of a sphere I

$$ V = \int_0^R\int_0^\sqrt{1-x^2}\int_0^\sqrt{1-x^2-y^2} \, dz\, dy\, dx \equiv \int_0^R\int_{g(x)}^{h(x)}\int_{q(x,y)}^{r(x,y)}\, dz\, dy\, dx$$

With $g$, $h$, $q$, and $r$ defined normally:

In [9]:
def func4(z,y,x):    # ORDER OF ARGUMENTS IMPORTANT
    return 1

def g(x):
    return 0

def h(x):
    return np.sqrt(xulim**2-x**2)

def q(x,y):
    return 0

def r(x,y):
    return np.sqrt(xulim**2-x**2-y**2)
In [10]:
xllim = 0
xulim = 2
value, error = integrate.tplquad(func4, xllim, xulim, g, h, q, r)
8*value, error, 4*np.pi*xulim**3/3.
Out[10]:
(33.51032163829113, 3.533511261366584e-10, 33.510321638291124)

With $g$, $h$, $q$, and $r$ defined more concisely:

In [11]:
g = lambda x: 0
h = lambda x: np.sqrt(xulim**2-x**2)
q = lambda x,y: 0
r = lambda x,y: np.sqrt(xulim**2-x**2-y**2)
In [12]:
xllim = 0
xulim = 2
value, error = integrate.tplquad(func4, xllim, xulim, g, h, q, r)
8*value, error, 4*np.pi*xulim**3/3.
Out[12]:
(33.51032163829113, 3.533511261366584e-10, 33.510321638291124)

Or even more concisely:

In [13]:
xllim = 0
xulim = 2
value, error = integrate.tplquad(func4, xllim, xulim, \
                                    lambda x:0, lambda x:np.sqrt(xulim**2-x**2), \
                                    lambda x,y:0, lambda x,y:np.sqrt(xulim**2-x**2-y**2))
8*value, error, 4*np.pi*xulim**3/3.
Out[13]:
(33.51032163829113, 3.533511261366584e-10, 33.510321638291124)

Volume of a sphere I I

$$ V = \int_0^{2\pi}\int_0^\pi\int_0^R r^2 \sin\theta\, dr d\theta d\phi $$
In [14]:
def func5(phi,theta,r):
    return r**2*np.sin(theta)
In [15]:
rllim = 0
rulim = 2
value, error = integrate.tplquad(func5, rllim, rulim, \
                                    lambda theta:0, lambda theta:np.pi, \
                                    lambda theta,phi:0, lambda theta,phi:2.*np.pi)
value, error, 4*np.pi*xulim**3/3.
Out[15]:
(33.510321638291124, 5.556380318779972e-13, 33.510321638291124)

Version information

version_information is from J.R. Johansson (jrjohansson at gmail.com); see Introduction to scientific computing with Python for more information and instructions for package installation.

version_information is installed on the linux network at Bucknell

In [16]:
%load_ext version_information
In [17]:
version_information numpy, scipy, matplotlib
Out[17]:
SoftwareVersion
Python3.7.7 64bit [GCC 7.3.0]
IPython7.16.1
OSLinux 3.10.0 1062.9.1.el7.x86_64 x86_64 with centos 7.7.1908 Core
numpy1.18.5
scipy1.5.2
matplotlib3.3.0
Fri Aug 07 15:21:07 2020 EDT
In [ ]: