Jupyter Snippet NP 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).

%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'

import numpy as np

from scipy import integrate

import sympy
import mpmath

sympy.init_printing()

Simpson’s rule

a, b, X = sympy.symbols("a, b, x")
f = sympy.Function("f")

#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

w = [sympy.symbols("w_%d" % i) for i in range(len(x))] 

q_rule = sum([w[i] * f(x[i]) for i in range(len(x))])

q_rule

$\displaystyle w_{0} f{\left(a \right)} + w_{1} f{\left(\frac{a}{2} + \frac{b}{2} \right)} + w_{2} f{\left(b \right)}$

phi = [sympy.Lambda(X, X**n) for n in range(len(x))]

phi

$\displaystyle \left[ \left( x \mapsto 1 \right), \ \left( x \mapsto x \right), \ \left( x \mapsto x^{2} \right)\right]$

eqs = [q_rule.subs(f, phi[n]) - sympy.integrate(phi[n](X), (X, a, b)) for n in range(len(phi))]

eqs

$\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]$

w_sol = sympy.solve(eqs, w)

w_sol

$\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}$

q_rule.subs(w_sol).simplify()

$\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

def f(x):
    return np.exp(-x**2)

val, err = integrate.quad(f, -1, 1)

val

$\displaystyle 1.493648265624854$

err

$\displaystyle 1.6582826951881447e-14$

val, err = integrate.quadrature(f, -1, 1)

val

$\displaystyle 1.4936482656450039$

err

$\displaystyle 7.459897144457273e-10$

Extra arguments

def f(x, a, b, c):
    return a * np.exp(-((x-b)/c)**2)

val, err = integrate.quad(f, -1, 1, args=(1, 2, 3))

val

$\displaystyle 1.2763068351022229$

err

$\displaystyle 1.4169852348169507e-14$

Reshuffle arguments

from scipy.special import jv

val, err = integrate.quad(lambda x: jv(0, x), 0, 5)

val

$\displaystyle 0.7153119177847678$

err

$\displaystyle 2.47260738289741e-14$

Infinite limits

f = lambda x: np.exp(-x**2)

val, err = integrate.quad(f, -np.inf, np.inf)

val

$\displaystyle 1.7724538509055159$

err

$\displaystyle 1.4202636780944923e-08$

Singularity

f = lambda x: 1/np.sqrt(abs(x))

a, b = -1, 1

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.

$\displaystyle \left( inf, \ inf\right)$

integrate.quad(f, a, b, points=[0])

$\displaystyle \left( 3.9999999999999813, \ 5.684341886080802e-14\right)$

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

f = lambda x: np.sqrt(x)

a, b = 0, 2

x = np.linspace(a, b, 25)

y = f(x)

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")

val_trapz = integrate.trapz(y, x)

val_trapz

$\displaystyle 1.8808217160508505$

val_simps = integrate.simps(y, x)

val_simps

$\displaystyle 1.883665102448715$

val_exact = 2.0/3.0 * (b-a)**(3.0/2.0)

val_exact

$\displaystyle 1.8856180831641267$

val_exact - val_trapz

$\displaystyle 0.004796367113276245$

val_exact - val_simps

$\displaystyle 0.001952980715411723$

x = np.linspace(a, b, 1 + 2**6)

len(x)

$\displaystyle 65$

y = f(x)

val_exact - integrate.romb(y, dx=(x[1]-x[0]))

$\displaystyle 0.00037879842291310695$

val_exact - integrate.simps(y, dx=x[1]-x[0])

$\displaystyle 0.0004484855541582178$

Higher dimension

def f(x):
    return np.exp(-x**2)

%time integrate.quad(f, a, b)

CPU times: user 160 µs, sys: 40 µs, total: 200 µs
Wall time: 214 µs

$\displaystyle \left( 0.8820813907624215, \ 9.793070696178202e-15\right)$

def f(x, y):
    return np.exp(-x**2-y**2)

a, b = 0, 1

g = lambda x: 0

h = lambda x: 1

integrate.dblquad(f, a, b, g, h)

$\displaystyle \left( 0.5577462853510337, \ 8.291374381535408e-15\right)$

integrate.dblquad(lambda x, y: np.exp(-x**2-y**2), 0, 1, lambda x: 0, lambda x: 1)

$\displaystyle \left( 0.5577462853510337, \ 8.291374381535408e-15\right)$

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")

integrate.dblquad(f, 0, 1, lambda x: -1 + x, lambda x: 1 - x)

$\displaystyle \left( 0.7320931000008094, \ 1.6564972931774035e-14\right)$

def f(x, y, z):
    return np.exp(-x**2-y**2-z**2)

integrate.tplquad(f, 0, 1, lambda x : 0, lambda x : 1, lambda x, y : 0, lambda x, y : 1)

$\displaystyle \left( 0.4165383858866382, \ 8.291335287314424e-15\right)$

integrate.nquad(f, [(0, 1), (0, 1), (0, 1)])

$\displaystyle \left( 0.4165383858866382, \ 8.291335287314424e-15\right)$

nquad

def f(*args):
    return  np.exp(-np.sum(np.array(args)**2))

%time integrate.nquad(f, [(0,1)] * 1)

CPU times: user 677 µs, sys: 79 µs, total: 756 µs
Wall time: 767 µs

$\displaystyle \left( 0.7468241328124271, \ 8.291413475940725e-15\right)$

%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

$\displaystyle \left( 0.5577462853510337, \ 8.291374381535408e-15\right)$

%time integrate.nquad(f, [(0,1)] * 3)

CPU times: user 169 ms, sys: 3.38 ms, total: 172 ms
Wall time: 171 ms

$\displaystyle \left( 0.4165383858866382, \ 8.291335287314424e-15\right)$

%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

$\displaystyle \left( 0.31108091882287664, \ 8.291296193277774e-15\right)$

%time integrate.nquad(f, [(0,1)] * 5)

CPU times: user 1min 9s, sys: 237 ms, total: 1min 9s
Wall time: 1min 9s

$\displaystyle \left( 0.23232273743438786, \ 8.29125709942545e-15\right)$

Monte Carlo integration

from skmonaco import mcquad

%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

val, err

$\displaystyle \left( 0.2315775914499913, \ 0.00047512221420657915\right)$

%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

val, err

$\displaystyle \left( 0.0540332249061147, \ 0.0001721803643623528\right)$

Symbolic and multi-precision quadrature

x = sympy.symbols("x")

f = 2 * sympy.sqrt(1-x**2)

a, b = -1, 1

sympy.plot(f, (x, -2, 2));

val_sym = sympy.integrate(f, (x, a, b))

val_sym

$\displaystyle \pi$

mpmath.mp.dps = 75

f_mpmath = sympy.lambdify(x, f, 'mpmath')

val = mpmath.quad(f_mpmath, (a, b))

sympy.sympify(val)

$\displaystyle 3.14159265358979323846264338327950288419716939937510582097494459230781640629$

sympy.N(val_sym, mpmath.mp.dps+1) - val

$\displaystyle 6.90893484407555570030908149024031965689280029154902510801896277613487344253 \cdot 10^{-77}$

%timeit mpmath.quad(f_mpmath, [a, b])

7.26 ms ± 82.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

f_numpy = sympy.lambdify(x, f, 'numpy')

%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

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)

integrate.dblquad(f2, 0, 1, lambda x : 0, lambda x : 1)

$\displaystyle \left( 0.4305647943060992, \ 7.28494733311039e-15\right)$

integrate.tplquad(f3, 0, 1, lambda x : 0, lambda x : 1, lambda x, y : 0, lambda x, y : 1)

$\displaystyle \left( 0.2825255795184269, \ 7.2848958098526e-15\right)$

x, y, z = sympy.symbols("x, y, z")

f2 = sympy.cos(x)*sympy.cos(y)*sympy.exp(-x**2-y**2)

f3 = sympy.cos(x)*sympy.cos(y)*sympy.cos(z) * sympy.exp(-x**2 - y**2 - z**2)

sympy.integrate(f3, (x, 0, 1), (y, 0, 1), (z, 0, 1))  # this does not succeed

$\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$

f2_numpy = sympy.lambdify((x, y), f2, 'numpy')

integrate.dblquad(f2_numpy, 0, 1, lambda x: 0, lambda x: 1)

$\displaystyle \left( 0.4305647943060992, \ 7.28494733311039e-15\right)$

f3_numpy = sympy.lambdify((x, y, z), f3, 'numpy')

integrate.tplquad(f3_numpy, 0, 1, lambda x: 0, lambda x: 1, lambda x, y: 0, lambda x, y: 1)

$\displaystyle \left( 0.2825255795184269, \ 7.2848958098526e-15\right)$

mpmath.mp.dps = 30

f2_mpmath = sympy.lambdify((x, y), f2, 'mpmath')

res = mpmath.quad(f2_mpmath, (0, 1), (0, 1))
res

mpf('0.430564794306099099242308990195783')

f3_mpmath = sympy.lambdify((x, y, z), f3, 'mpmath')

res = mpmath.quad(f3_mpmath, (0, 1), (0, 1), (0, 1))

sympy.sympify(res)

$\displaystyle 0.282525579518426896867622772405$

%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

t, x, y = sympy.symbols("t, x, y")

C = sympy.Curve([sympy.cos(t), sympy.sin(t)], (t, 0, 2 * sympy.pi))

sympy.line_integrate(1, C, [x, y])

$\displaystyle 2 \pi$

sympy.line_integrate(x**2 * y**2, C, [x, y])

$\displaystyle \frac{\pi}{4}$

Integral transformations

Laplace transforms

s = sympy.symbols("s")

a, t = sympy.symbols("a, t", positive=True)

f = sympy.sin(a*t)

sympy.laplace_transform(f, t, s)

$\displaystyle \left( \frac{a}{a^{2} + s^{2}}, \ -\infty, \ \operatorname{re}{\left(s\right)} > 0\right)$

F = sympy.laplace_transform(f, t, s, noconds=True)

F

$\displaystyle \frac{a}{a^{2} + s^{2}}$

sympy.inverse_laplace_transform(F, s, t, noconds=True)

$\displaystyle \sin{\left(a t \right)}$

[sympy.laplace_transform(f, t, s, noconds=True) for f in [t, t**2, t**3, t**4]]

$\displaystyle \left[ \frac{1}{s^{2}}, \ \frac{2}{s^{3}}, \ \frac{6}{s^{4}}, \ \frac{24}{s^{5}}\right]$

n = sympy.symbols("n", integer=True, positive=True)

sympy.laplace_transform(t**n, t, s, noconds=True)

$\displaystyle s^{- n - 1} n!$

sympy.laplace_transform((1 - a*t) * sympy.exp(-a*t), t, s, noconds=True)

$\displaystyle \frac{s}{\left(a + s\right)^{2}}$

Fourier Transforms

w = sympy.symbols("omega")

f = sympy.exp(-a*t**2)

F = sympy.fourier_transform(f, t, w)

F

$\displaystyle \frac{\sqrt{\pi} e^{- \frac{\pi^{2} \omega^{2}}{a}}}{\sqrt{a}}$

sympy.inverse_fourier_transform(F, w, t)

$\displaystyle e^{- a t^{2}}$

sympy.fourier_transform(sympy.cos(t), t, w)  # not good

$\displaystyle 0$

Versions

%reload_ext version_information

%version_information numpy, matplotlib, scipy, sympy, mpmath, skmonaco