ch08-code-listing

Chapter 8: Integration

Robert Johansson

Source code listings for Numerical Python - Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib (ISBN 978-1-484242-45-2).

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl

mpl.rcParams['mathtext.fontset'] = 'stix'
mpl.rcParams['font.family'] = 'serif'
mpl.rcParams['font.sans-serif'] = 'stix'
In [2]:
import numpy as np
In [3]:
from scipy import integrate
In [4]:
import sympy
import mpmath
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]:
$\displaystyle 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]:
$\displaystyle \left[ \left( x \mapsto 1 \right), \ \left( x \mapsto x \right), \ \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]:
$\displaystyle \left[ a - b + w_{0} + w_{1} + w_{2}, \ \frac{a^{2}}{2} + a w_{0} - \frac{b^{2}}{2} + b w_{2} + w_{1} \left(\frac{a}{2} + \frac{b}{2}\right), \ \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]:
$\displaystyle \left\{ w_{0} : - \frac{a}{6} + \frac{b}{6}, \ w_{1} : - \frac{2 a}{3} + \frac{2 b}{3}, \ w_{2} : - \frac{a}{6} + \frac{b}{6}\right\}$
In [17]:
q_rule.subs(w_sol).simplify()
Out[17]:
$\displaystyle - \frac{\left(a - b\right) \left(f{\left(a \right)} + f{\left(b \right)} + 4 f{\left(\frac{a}{2} + \frac{b}{2} \right)}\right)}{6}$

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]:
$\displaystyle 1.493648265624854$
In [21]:
err
Out[21]:
$\displaystyle 1.6582826951881447e-14$
In [22]:
val, err = integrate.quadrature(f, -1, 1)
In [23]:
val
Out[23]:
$\displaystyle 1.4936482656450039$
In [24]:
err
Out[24]:
$\displaystyle 7.459897144457273e-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]:
$\displaystyle 1.2763068351022229$
In [28]:
err
Out[28]:
$\displaystyle 1.4169852348169507e-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]:
$\displaystyle 0.7153119177847678$
In [32]:
err
Out[32]:
$\displaystyle 2.47260738289741e-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]:
$\displaystyle 1.7724538509055159$
In [36]:
err
Out[36]:
$\displaystyle 1.4202636780944923e-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/miniconda3/envs/py3.6/lib/python3.6/site-packages/ipykernel_launcher.py:1: RuntimeWarning: divide by zero encountered in double_scalars
  """Entry point for launching an IPython kernel.
Out[39]:
$\displaystyle \left( inf, \ inf\right)$
In [40]:
integrate.quad(f, a, b, points=[0])
Out[40]:
$\displaystyle \left( 3.9999999999999813, \ 5.684341886080802e-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)
ax.set_xlim(-1, 1)

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]:
$\displaystyle 1.8808217160508505$
In [49]:
val_simps = integrate.simps(y, x)
In [50]:
val_simps
Out[50]:
$\displaystyle 1.883665102448715$
In [51]:
val_exact = 2.0/3.0 * (b-a)**(3.0/2.0)
In [52]:
val_exact
Out[52]:
$\displaystyle 1.8856180831641267$
In [53]:
val_exact - val_trapz
Out[53]:
$\displaystyle 0.004796367113276245$
In [54]:
val_exact - val_simps
Out[54]:
$\displaystyle 0.001952980715411723$
In [55]:
x = np.linspace(a, b, 1 + 2**6)
In [56]:
len(x)
Out[56]:
$\displaystyle 65$
In [57]:
y = f(x)
In [58]:
val_exact - integrate.romb(y, dx=(x[1]-x[0]))
Out[58]:
$\displaystyle 0.00037879842291310695$
In [59]:
val_exact - integrate.simps(y, dx=x[1]-x[0])
Out[59]:
$\displaystyle 0.0004484855541582178$

Higher dimension

In [60]:
def f(x):
    return np.exp(-x**2)
In [61]:
%time integrate.quad(f, a, b)
CPU times: user 160 µs, sys: 40 µs, total: 200 µs
Wall time: 214 µs
Out[61]:
$\displaystyle \left( 0.8820813907624215, \ 9.793070696178202e-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]:
$\displaystyle \left( 0.5577462853510337, \ 8.291374381535408e-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]:
$\displaystyle \left( 0.5577462853510337, \ 8.291374381535408e-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]:
$\displaystyle \left( 0.7320931000008094, \ 1.6564972931774035e-14\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]:
$\displaystyle \left( 0.4165383858866382, \ 8.291335287314424e-15\right)$
In [72]:
integrate.nquad(f, [(0, 1), (0, 1), (0, 1)])
Out[72]:
$\displaystyle \left( 0.4165383858866382, \ 8.291335287314424e-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 677 µs, sys: 79 µs, total: 756 µs
Wall time: 767 µs
Out[74]:
$\displaystyle \left( 0.7468241328124271, \ 8.291413475940725e-15\right)$
In [75]:
%time integrate.nquad(f, [(0,1)] * 2)
CPU times: user 8.53 ms, sys: 690 µs, total: 9.22 ms
Wall time: 8.75 ms
Out[75]:
$\displaystyle \left( 0.5577462853510337, \ 8.291374381535408e-15\right)$
In [76]:
%time integrate.nquad(f, [(0,1)] * 3)
CPU times: user 169 ms, sys: 3.38 ms, total: 172 ms
Wall time: 171 ms
Out[76]:
$\displaystyle \left( 0.4165383858866382, \ 8.291335287314424e-15\right)$
In [77]:
%time integrate.nquad(f, [(0,1)] * 4)
CPU times: user 3.31 s, sys: 14.7 ms, total: 3.32 s
Wall time: 3.32 s
Out[77]:
$\displaystyle \left( 0.31108091882287664, \ 8.291296193277774e-15\right)$
In [78]:
%time integrate.nquad(f, [(0,1)] * 5)
CPU times: user 1min 9s, sys: 237 ms, total: 1min 9s
Wall time: 1min 9s
Out[78]:
$\displaystyle \left( 0.23232273743438786, \ 8.29125709942545e-15\right)$

Monte Carlo integration

In [80]:
from skmonaco import mcquad
In [81]:
%time val, err = mcquad(f, xl=np.zeros(5), xu=np.ones(5), npoints=100000)
CPU times: user 2.01 s, sys: 160 ms, total: 2.17 s
Wall time: 2.05 s
In [82]:
val, err
Out[82]:
$\displaystyle \left( 0.2315775914499913, \ 0.00047512221420657915\right)$
In [83]:
%time val, err = mcquad(f, xl=np.zeros(10), xu=np.ones(10), npoints=100000)
CPU times: user 1.93 s, sys: 101 ms, total: 2.03 s
Wall time: 1.96 s
In [84]:
val, err
Out[84]:
$\displaystyle \left( 0.0540332249061147, \ 0.0001721803643623528\right)$

Symbolic and multi-precision quadrature

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

double and triple integrals

In [99]:
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 [100]:
integrate.dblquad(f2, 0, 1, lambda x : 0, lambda x : 1)
Out[100]:
$\displaystyle \left( 0.4305647943060992, \ 7.28494733311039e-15\right)$
In [101]:
integrate.tplquad(f3, 0, 1, lambda x : 0, lambda x : 1, lambda x, y : 0, lambda x, y : 1)
Out[101]:
$\displaystyle \left( 0.2825255795184269, \ 7.2848958098526e-15\right)$
In [102]:
x, y, z = sympy.symbols("x, y, z")
In [103]:
f2 = sympy.cos(x)*sympy.cos(y)*sympy.exp(-x**2-y**2)
In [104]:
f3 = sympy.cos(x)*sympy.cos(y)*sympy.cos(z) * sympy.exp(-x**2 - y**2 - z**2)
In [105]:
sympy.integrate(f3, (x, 0, 1), (y, 0, 1), (z, 0, 1))  # this does not succeed
Out[105]:
$\displaystyle \left(\int\limits_{0}^{1} e^{- x^{2}} \cos{\left(x \right)}\, dx\right) \left(\int\limits_{0}^{1} e^{- y^{2}} \cos{\left(y \right)}\, dy\right) \int\limits_{0}^{1} e^{- z^{2}} \cos{\left(z \right)}\, dz$
In [106]:
f2_numpy = sympy.lambdify((x, y), f2, 'numpy')
In [107]:
integrate.dblquad(f2_numpy, 0, 1, lambda x: 0, lambda x: 1)
Out[107]:
$\displaystyle \left( 0.4305647943060992, \ 7.28494733311039e-15\right)$
In [108]:
f3_numpy = sympy.lambdify((x, y, z), f3, 'numpy')
In [109]:
integrate.tplquad(f3_numpy, 0, 1, lambda x: 0, lambda x: 1, lambda x, y: 0, lambda x, y: 1)
Out[109]:
$\displaystyle \left( 0.2825255795184269, \ 7.2848958098526e-15\right)$
In [110]:
mpmath.mp.dps = 30
In [111]:
f2_mpmath = sympy.lambdify((x, y), f2, 'mpmath')
In [112]:
res = mpmath.quad(f2_mpmath, (0, 1), (0, 1))
res
Out[112]:
mpf('0.430564794306099099242308990195783')
In [113]:
f3_mpmath = sympy.lambdify((x, y, z), f3, 'mpmath')
In [114]:
res = mpmath.quad(f3_mpmath, (0, 1), (0, 1), (0, 1))
In [115]:
sympy.sympify(res)
Out[115]:
$\displaystyle 0.282525579518426896867622772405$
In [116]:
%time res = sympy.sympify(mpmath.quad(f3_mpmath, (0, 1), (0, 1), (0, 1)))
CPU times: user 2min 44s, sys: 294 ms, total: 2min 44s
Wall time: 2min 44s

Line integrals

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

Integral transformations

Laplace transforms

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

Fourier Transforms

In [132]:
w = sympy.symbols("omega")
In [133]:
f = sympy.exp(-a*t**2)
In [134]:
F = sympy.fourier_transform(f, t, w)
In [135]:
F
Out[135]:
$\displaystyle \frac{\sqrt{\pi} e^{- \frac{\pi^{2} \omega^{2}}{a}}}{\sqrt{a}}$
In [136]:
sympy.inverse_fourier_transform(F, w, t)
Out[136]:
$\displaystyle e^{- a t^{2}}$
In [137]:
sympy.fourier_transform(sympy.cos(t), t, w)  # not good
Out[137]:
$\displaystyle 0$

Versions

In [138]:
%reload_ext version_information
In [139]:
%version_information numpy, matplotlib, scipy, sympy, mpmath, skmonaco
Out[139]:
SoftwareVersion
Python3.6.8 64bit [GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]
IPython7.5.0
OSDarwin 18.2.0 x86_64 i386 64bit
numpy1.16.3
matplotlib3.0.3
scipy1.2.1
sympy1.4
mpmath1.1.0
skmonaco0.2.1
Mon May 06 14:55:36 2019 JST