Numerical integration

In [1]:
import scipy as sp
from scipy import integrate   # not included in basic scipy

import matplotlib as mpl       # As of July 2017 Bucknell computers use v. 2.x 
import matplotlib.pyplot as plt

# Following is an Ipython magic command that puts figures in the  notebook.
# For figures in separate windows, comment out following line and uncomment
# the next line
# Must come before defaults are changed.
%matplotlib notebook
#%matplotlib

# As of Aug. 2017 reverting to 1.x defaults.
# In 2.x text.ustex requires dvipng, texlive-latex-extra, and texlive-fonts-recommended, 
# which don't seem to be universal
# See https://stackoverflow.com/questions/38906356/error-running-matplotlib-in-latex-type1cm?
mpl.style.use('classic')
        
# 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[]
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 [2]:
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 sp.sin(x-2)/(x-2)
In [3]:
value, error = sp.integrate.quad(func1,1,2)
value, error
Out[3]:
(3.7500000000000004, 4.1633363423443377e-14)

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

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

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

/usr/remote/anaconda-3.6/lib/python3.6/site-packages/scipy/integrate/quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    321     if (weight is None):
    322         retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
--> 323                        points)
    324     else:
    325         retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel,

/usr/remote/anaconda-3.6/lib/python3.6/site-packages/scipy/integrate/quadpack.py in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    386     if points is None:
    387         if infbounds == 0:
--> 388             return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
    389         else:
    390             return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

TypeError: must be real number, not NoneType

Specify troublesome points:

In [5]:
value, error = sp.integrate.quad(func2,1,3,points=[2,])
value, error
Out[5]:
(8.666666666666668, 9.621932880084691e-14)
In [6]:
value, error = sp.integrate.quad(func3,0,4,points=[2,])
value, error
Out[6]:
(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 [7]:
def func4(z,y,x):    # ORDER OF ARGUMENTS IMPORTANT
    return 1

def g(x):
    return 0

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

def q(x,y):
    return 0

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

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

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

Or even more concisely:

In [11]:
xllim = 0
xulim = 2
value, error = sp.integrate.tplquad(func4, xllim, xulim, \
                                    lambda x:0, lambda x:sp.sqrt(xulim**2-x**2), \
                                    lambda x,y:0, lambda x,y:sp.sqrt(xulim**2-x**2-y**2))
8*value, error, 4*sp.pi*xulim**3/3.
Out[11]:
(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 [12]:
def func5(phi,theta,r):
    return r**2*sp.sin(theta)
In [13]:
rllim = 0
rulim = 2
value, error = sp.integrate.tplquad(func5, rllim, rulim, \
                                    lambda theta:0, lambda theta:sp.pi, \
                                    lambda theta,phi:0, lambda theta,phi:2.*sp.pi)
value, error, 4*sp.pi*xulim**3/3.
Out[13]:
(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: http://nbviewer.jupyter.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-0-Scientific-Computing-with-Python.ipynb for more information and instructions for package installation. If version_information has been installed system wide (as it has been on linuxremotes), continue with next cell as written. If not, comment out top line in next cell and uncomment the second line.

In [14]:
%load_ext version_information

#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py

Loading extensions from ~/.ipython/extensions is deprecated. We recommend managing extensions like any other Python packages, in site-packages.
In [15]:
version_information scipy, matplotlib
Out[15]:
SoftwareVersion
Python3.6.1 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython6.1.0
OSLinux 3.10.0 327.36.3.el7.x86_64 x86_64 with redhat 7.2 Maipo
scipy0.19.1
matplotlib2.0.2
Tue Aug 01 11:18:58 2017 EDT
In [ ]: