Chapter 8: Integration

Robert Johansson

Source code listings for Numerical Python - A Practical Techniques Approach for Industry (ISBN 978-1-484205-54-9).

The source code listings can be downloaded from http://www.apress.com/9781484205549

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl
In [2]:
import numpy as np
In [3]:
from scipy import integrate
In [4]:
import sympy
In [5]:
sympy.init_printing()

Simpson's rule

In [6]:
a, b, X = sympy.symbols("a, b, x")
f = sympy.Function("f")
In [7]:
#x = a, (a+b)/3, 2 * (a+b)/3, b # 3rd order quadrature rule
x = a, (a+b)/2, b # simpson's rule
#x = a, b # trapezoid rule
#x = ((b+a)/2,)  # mid-point rule
In [8]:
w = [sympy.symbols("w_%d" % i) for i in range(len(x))] 
In [9]:
q_rule = sum([w[i] * f(x[i]) for i in range(len(x))])
In [10]:
q_rule
Out[10]:
$$w_{0} f{\left (a \right )} + w_{1} f{\left (\frac{a}{2} + \frac{b}{2} \right )} + w_{2} f{\left (b \right )}$$
In [11]:
phi = [sympy.Lambda(X, X**n) for n in range(len(x))]
In [12]:
phi
Out[12]:
$$\left [ \left( x \mapsto 1 \right), \quad \left( x \mapsto x \right), \quad \left( x \mapsto x^{2} \right)\right ]$$
In [13]:
eqs = [q_rule.subs(f, phi[n]) - sympy.integrate(phi[n](X), (X, a, b)) for n in range(len(phi))]
In [14]:
eqs
Out[14]:
$$\left [ a - b + w_{0} + w_{1} + w_{2}, \quad \frac{a^{2}}{2} + a w_{0} - \frac{b^{2}}{2} + b w_{2} + w_{1} \left(\frac{a}{2} + \frac{b}{2}\right), \quad \frac{a^{3}}{3} + a^{2} w_{0} - \frac{b^{3}}{3} + b^{2} w_{2} + w_{1} \left(\frac{a}{2} + \frac{b}{2}\right)^{2}\right ]$$
In [15]:
w_sol = sympy.solve(eqs, w)
In [16]:
w_sol
Out[16]:
$$\left \{ w_{0} : - \frac{a}{6} + \frac{b}{6}, \quad w_{1} : - \frac{2 a}{3} + \frac{2 b}{3}, \quad w_{2} : - \frac{a}{6} + \frac{b}{6}\right \}$$
In [17]:
q_rule.subs(w_sol).simplify()
Out[17]:
$$- \frac{1}{6} \left(a - b\right) \left(f{\left (a \right )} + f{\left (b \right )} + 4 f{\left (\frac{a}{2} + \frac{b}{2} \right )}\right)$$

SciPy integrate

Simple integration example

In [18]:
def f(x):
    return np.exp(-x**2)
In [19]:
val, err = integrate.quad(f, -1, 1)
In [20]:
val
Out[20]:
$$1.49364826562$$
In [21]:
err
Out[21]:
$$1.65828269519e-14$$
In [22]:
val, err = integrate.quadrature(f, -1, 1)
In [23]:
val
Out[23]:
$$1.49364826565$$
In [24]:
err
Out[24]:
$$7.45989714446e-10$$

Extra arguments

In [25]:
def f(x, a, b, c):
    return a * np.exp(-((x-b)/c)**2)
In [26]:
val, err = integrate.quad(f, -1, 1, args=(1, 2, 3))
In [27]:
val
Out[27]:
$$1.2763068351$$
In [28]:
err
Out[28]:
$$1.41698523482e-14$$

Reshuffle arguments

In [29]:
from scipy.special import jv
In [30]:
val, err = integrate.quad(lambda x: jv(0, x), 0, 5)
In [31]:
val
Out[31]:
$$0.715311917785$$
In [32]:
err
Out[32]:
$$2.4726073829e-14$$

Infinite limits

In [33]:
f = lambda x: np.exp(-x**2)
In [34]:
val, err = integrate.quad(f, -np.inf, np.inf)
In [35]:
val
Out[35]:
$$1.77245385091$$
In [36]:
err
Out[36]:
$$1.42026367809e-08$$

Singularity

In [37]:
f = lambda x: 1/np.sqrt(abs(x))
In [38]:
a, b = -1, 1
In [39]:
integrate.quad(f, a, b)
/Users/rob/miniconda/envs/py27-npm/lib/python2.7/site-packages/IPython/kernel/__main__.py:1: RuntimeWarning: divide by zero encountered in double_scalars
  if __name__ == '__main__':
Out[39]:
$$\left ( inf, \quad inf\right )$$
In [40]:
integrate.quad(f, a, b, points=[0])
Out[40]:
$$\left ( 4.0, \quad 5.68434188608e-14\right )$$
In [41]:
fig, ax = plt.subplots(figsize=(8, 3))

x = np.linspace(a, b, 10000)
ax.plot(x, f(x), lw=2)
ax.fill_between(x, f(x), color='green', alpha=0.5)
ax.set_xlabel("$x$", fontsize=18)
ax.set_ylabel("$f(x)$", fontsize=18)
ax.set_ylim(0, 25)

fig.tight_layout()
fig.savefig("ch8-diverging-integrand.pdf")

Tabulated integrand

In [42]:
f = lambda x: np.sqrt(x)
In [43]:
a, b = 0, 2
In [44]:
x = np.linspace(a, b, 25)
In [45]:
y = f(x)
In [46]:
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(x, y, 'bo')
xx = np.linspace(a, b, 500)
ax.plot(xx, f(xx), 'b-')
ax.fill_between(xx, f(xx), color='green', alpha=0.5)
ax.set_xlabel(r"$x$", fontsize=18)
ax.set_ylabel(r"$f(x)$", fontsize=18)
fig.tight_layout()
fig.savefig("ch8-tabulated-integrand.pdf")
In [47]:
val_trapz = integrate.trapz(y, x)
In [48]:
val_trapz
Out[48]:
$$1.88082171605$$
In [49]:
val_simps = integrate.simps(y, x)
In [50]:
val_simps
Out[50]:
$$1.88366510245$$
In [51]:
val_exact = 2.0/3.0 * (b-a)**(3.0/2.0)
In [52]:
val_exact
Out[52]:
$$1.88561808316$$
In [53]:
val_exact - val_trapz
Out[53]:
$$0.00479636711328$$
In [54]:
val_exact - val_simps
Out[54]:
$$0.00195298071541$$
In [55]:
x = np.linspace(a, b, 1 + 2**6)
In [56]:
len(x)
Out[56]:
$$65$$
In [57]:
y = f(x)
In [58]:
val_exact - integrate.romb(y, dx=(x[1]-x[0]))
Out[58]:
$$0.000378798422913$$
In [59]:
val_exact - integrate.simps(y, dx=x[1]-x[0])
Out[59]:
$$0.000448485554158$$

Higher dimension

In [60]:
def f(x):
    return np.exp(-x**2)
In [61]:
%time integrate.quad(f, a, b)
CPU times: user 98 µs, sys: 46 µs, total: 144 µs
Wall time: 142 µs
Out[61]:
$$\left ( 0.882081390762, \quad 9.79307069618e-15\right )$$
In [62]:
def f(x, y):
    return np.exp(-x**2-y**2)
In [63]:
a, b = 0, 1
In [64]:
g = lambda x: 0
In [65]:
h = lambda x: 1
In [66]:
integrate.dblquad(f, a, b, g, h)
Out[66]:
$$\left ( 0.557746285351, \quad 6.19222767896e-15\right )$$
In [67]:
integrate.dblquad(lambda x, y: np.exp(-x**2-y**2), 0, 1, lambda x: 0, lambda x: 1)
Out[67]:
$$\left ( 0.557746285351, \quad 6.19222767896e-15\right )$$
In [68]:
fig, ax = plt.subplots(figsize=(6, 5))

x = y = np.linspace(-1.25, 1.25, 75)
X, Y = np.meshgrid(x, y)

c = ax.contour(X, Y, f(X, Y), 15, cmap=mpl.cm.RdBu, vmin=-1, vmax=1)

bound_rect = plt.Rectangle((0, 0), 1, 1,
                           facecolor="grey")
ax.add_patch(bound_rect)

ax.axis('tight')
ax.set_xlabel('$x$', fontsize=18)
ax.set_ylabel('$y$', fontsize=18)

fig.tight_layout()
fig.savefig("ch8-multi-dim-integrand.pdf")
In [69]:
integrate.dblquad(f, 0, 1, lambda x: -1 + x, lambda x: 1 - x)
Out[69]:
$$\left ( 0.732093100001, \quad 8.1278661579e-15\right )$$
In [70]:
def f(x, y, z):
    return np.exp(-x**2-y**2-z**2)
In [71]:
integrate.tplquad(f, 0, 1, lambda x : 0, lambda x : 1, lambda x, y : 0, lambda x, y : 1)
Out[71]:
$$\left ( 0.416538385887, \quad 4.62450506652e-15\right )$$
In [72]:
integrate.nquad(f, [(0, 1), (0, 1), (0, 1)])
Out[72]:
$$\left ( 0.416538385887, \quad 8.29133528731e-15\right )$$

nquad

In [73]:
def f(*args):
    return  np.exp(-np.sum(np.array(args)**2))
In [74]:
%time integrate.nquad(f, [(0,1)] * 1)
CPU times: user 361 µs, sys: 94 µs, total: 455 µs
Wall time: 440 µs
Out[74]:
$$\left ( 0.746824132812, \quad 8.29141347594e-15\right )$$
In [75]:
%time integrate.nquad(f, [(0,1)] * 2)
CPU times: user 10.5 ms, sys: 1.79 ms, total: 12.3 ms
Wall time: 11.6 ms
Out[75]:
$$\left ( 0.557746285351, \quad 8.29137438154e-15\right )$$
In [76]:
%time integrate.nquad(f, [(0,1)] * 3)
CPU times: user 130 ms, sys: 15.2 ms, total: 146 ms
Wall time: 135 ms
Out[76]:
$$\left ( 0.416538385887, \quad 8.29133528731e-15\right )$$
In [77]:
%time integrate.nquad(f, [(0,1)] * 4)
CPU times: user 2.59 s, sys: 43.5 ms, total: 2.63 s
Wall time: 2.67 s
Out[77]:
$$\left ( 0.311080918823, \quad 8.29129619328e-15\right )$$
In [78]:
%time integrate.nquad(f, [(0,1)] * 5)
CPU times: user 51.8 s, sys: 1.93 s, total: 53.7 s
Wall time: 52.8 s
Out[78]:
$$\left ( 0.232322737434, \quad 8.29125709943e-15\right )$$

Monte Carlo integration

In [79]:
from skmonaco import mcquad
In [80]:
%time val, err = mcquad(f, xl=np.zeros(5), xu=np.ones(5), npoints=100000)
CPU times: user 1.23 s, sys: 50 ms, total: 1.28 s
Wall time: 1.28 s
In [81]:
val, err
Out[81]:
$$\left ( 0.232011116953, \quad 0.000472864897239\right )$$
In [82]:
%time val, err = mcquad(f, xl=np.zeros(10), xu=np.ones(10), npoints=100000)
CPU times: user 1.23 s, sys: 49.1 ms, total: 1.28 s
Wall time: 1.27 s
In [83]:
val, err
Out[83]:
$$\left ( 0.0540787152508, \quad 0.000170681220268\right )$$

Symbolic and multi-precision quadrature

In [84]:
x = sympy.symbols("x")
In [85]:
f = 2 * sympy.sqrt(1-x**2)
In [86]:
a, b = -1, 1
In [87]:
sympy.plot(f, (x, -2, 2));
In [88]:
val_sym = sympy.integrate(f, (x, a, b))
In [89]:
val_sym
Out[89]:
$$\pi$$
In [90]:
sympy.mpmath.mp.dps = 75
In [91]:
f_mpmath = sympy.lambdify(x, f, 'mpmath')
In [92]:
val = sympy.mpmath.quad(f_mpmath, (a, b))
In [93]:
sympy.sympify(val)
Out[93]:
$$3.14159265358979323846264338327950288419716939937510582097494459230781640629$$
In [94]:
sympy.N(val_sym, sympy.mpmath.mp.dps+1) - val
Out[94]:
$$6.90893484407555570030908149024031965689280029154902510801896277613487344253 \cdot 10^{-77}$$
In [95]:
%timeit sympy.mpmath.quad(f_mpmath, [a, b])
100 loops, best of 3: 12.3 ms per loop
In [96]:
f_numpy = sympy.lambdify(x, f, 'numpy')
In [97]:
%timeit integrate.quad(f_numpy, a, b)
1000 loops, best of 3: 996 µs per loop

double and triple integrals

In [98]:
def f2(x, y):
    return np.cos(x)*np.cos(y)*np.exp(-x**2-y**2)

def f3(x, y, z):
    return np.cos(x)*np.cos(y)*np.cos(z)*np.exp(-x**2-y**2-z**2)
In [99]:
integrate.dblquad(f2, 0, 1, lambda x : 0, lambda x : 1)
Out[99]:
$$\left ( 0.430564794306, \quad 4.78022948232e-15\right )$$
In [100]:
integrate.tplquad(f3, 0, 1, lambda x : 0, lambda x : 1, lambda x, y : 0, lambda x, y : 1)
Out[100]:
$$\left ( 0.282525579518, \quad 3.13666403427e-15\right )$$
In [101]:
x, y, z = sympy.symbols("x, y, z")
In [102]:
f2 = sympy.cos(x)*sympy.cos(y)*sympy.exp(-x**2-y**2)
In [103]:
f3 = sympy.exp(-x**2 - y**2 - z**2)
In [104]:
f2_numpy = sympy.lambdify((x, y), f2, 'numpy')
In [105]:
integrate.dblquad(f2_numpy, 0, 1, lambda x: 0, lambda x: 1)
Out[105]:
$$\left ( 0.430564794306, \quad 4.78022948232e-15\right )$$
In [106]:
f3_numpy = sympy.lambdify((x, y, z), f3, 'numpy')
In [107]:
integrate.tplquad(f3_numpy, 0, 1, lambda x: 0, lambda x: 1, lambda x, y: 0, lambda x, y: 1)
Out[107]:
$$\left ( 0.416538385887, \quad 4.62450506652e-15\right )$$
In [108]:
sympy.mpmath.mp.dps = 30
In [109]:
f2_mpmath = sympy.lambdify((x, y), f2, 'mpmath')
In [110]:
sympy.mpmath.quad(f2_mpmath, (0, 1), (0, 1))
Out[110]:
mpf('0.430564794306099099242308990195783')
In [111]:
f3_mpmath = sympy.lambdify((x, y, z), f3, 'mpmath')
In [112]:
res = sympy.mpmath.quad(f3_mpmath, (0, 1), (0, 1), (0, 1))
In [114]:
sympy.sympify(res)
Out[114]:
$$0.416538385886638169609660243601$$

Line integrals

In [115]:
t, x, y = sympy.symbols("t, x, y")
In [116]:
C = sympy.Curve([sympy.cos(t), sympy.sin(t)], (t, 0, 2 * sympy.pi))
In [117]:
sympy.line_integrate(1, C, [x, y])
Out[117]:
$$2 \pi$$
In [118]:
sympy.line_integrate(x**2 * y**2, C, [x, y])
Out[118]:
$$\frac{\pi}{4}$$

Integral transformations

Laplace transforms

In [119]:
s = sympy.symbols("s")
In [120]:
a, t = sympy.symbols("a, t", positive=True)
In [121]:
f = sympy.sin(a*t)
In [122]:
sympy.laplace_transform(f, t, s)
Out[122]:
$$\left ( \frac{a}{a^{2} + s^{2}}, \quad -\infty, \quad 0 < \Re{s}\right )$$
In [123]:
F = sympy.laplace_transform(f, t, s, noconds=True)
In [124]:
F
Out[124]:
$$\frac{a}{a^{2} + s^{2}}$$
In [125]:
sympy.inverse_laplace_transform(F, s, t, noconds=True)
Out[125]:
$$\sin{\left (a t \right )}$$
In [126]:
[sympy.laplace_transform(f, t, s, noconds=True) for f in [t, t**2, t**3, t**4]]
Out[126]:
$$\left [ \frac{1}{s^{2}}, \quad \frac{2}{s^{3}}, \quad \frac{6}{s^{4}}, \quad \frac{24}{s^{5}}\right ]$$
In [127]:
n = sympy.symbols("n", integer=True, positive=True)
In [128]:
sympy.laplace_transform(t**n, t, s, noconds=True)
Out[128]:
$$\frac{\Gamma{\left(n + 1 \right)}}{s^{n + 1}}$$
In [129]:
sympy.laplace_transform((1 - a*t) * sympy.exp(-a*t), t, s, noconds=True)
Out[129]:
$$\frac{s}{\left(a + s\right)^{2}}$$

Fourier Transforms

In [130]:
w = sympy.symbols("omega")
In [131]:
f = sympy.exp(-a*t**2)
In [132]:
help(sympy.fourier_transform)
Help on function fourier_transform in module sympy.integrals.transforms:

fourier_transform(f, x, k, **hints)
    Compute the unitary, ordinary-frequency Fourier transform of `f`, defined
    as
    
    .. math:: F(k) = \int_{-\infty}^\infty f(x) e^{-2\pi i x k} \mathrm{d} x.
    
    If the transform cannot be computed in closed form, this
    function returns an unevaluated :class:`FourierTransform` object.
    
    For other Fourier transform conventions, see the function
    :func:`sympy.integrals.transforms._fourier_transform`.
    
    For a description of possible hints, refer to the docstring of
    :func:`sympy.integrals.transforms.IntegralTransform.doit`.
    Note that for this transform, by default ``noconds=True``.
    
    >>> from sympy import fourier_transform, exp
    >>> from sympy.abc import x, k
    >>> fourier_transform(exp(-x**2), x, k)
    sqrt(pi)*exp(-pi**2*k**2)
    >>> fourier_transform(exp(-x**2), x, k, noconds=False)
    (sqrt(pi)*exp(-pi**2*k**2), True)
    
    See Also
    ========
    
    inverse_fourier_transform
    sine_transform, inverse_sine_transform
    cosine_transform, inverse_cosine_transform
    hankel_transform, inverse_hankel_transform
    mellin_transform, laplace_transform

In [133]:
F = sympy.fourier_transform(f, t, w)
In [134]:
F
Out[134]:
$$\frac{\sqrt{\pi}}{\sqrt{a}} e^{- \frac{\pi^{2} \omega^{2}}{a}}$$
In [135]:
sympy.inverse_fourier_transform(F, w, t)
Out[135]:
$$e^{- a t^{2}}$$

Versions

In [136]:
%reload_ext version_information
In [137]:
%version_information numpy, matplotlib, scipy, sympy
Out[137]:
SoftwareVersion
Python2.7.10 64bit [GCC 4.2.1 (Apple Inc. build 5577)]
IPython3.2.1
OSDarwin 14.1.0 x86_64 i386 64bit
numpy1.9.2
matplotlib1.4.3
scipy0.16.0
sympy0.7.6