Skip to content

Latest commit

 

History

History
952 lines (466 loc) · 14.7 KB

ch08-integration.md

File metadata and controls

952 lines (466 loc) · 14.7 KB

Integration

  • Numerical Methods
  • Numerical Integration with SciPy
  • Tabulated Integrand
  • Multiple Integration
  • Symbolic & Arbitrary-Precision Integration
  • Line Integrals
  • Integral Transforms
  • Main focus: numerical integration (aka "quadrature").
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
from scipy import integrate
import sympy
sympy.init_printing()

Simpson's rule

  • Definite integrals - interpreted as area between integrand curve and x axis:

interpretation-as-area-between-curve-and-x-axis

  • Strategy for evaluating integral I(f): write it as discrete sum that approximates value of integral. ("n-point quadrature rule".)
  • Quadrature rules derived from f(x) interpolations along interval [a,b].
  • Midpoint rule: using 0th-order polynomial (constant value) of midpoint.
  • Trapezoid rule: using 1st-order polynomial, evalated at endpoints.
  • Simpson's rule: using 2nd-order polynomial, midpoint & endpoints.
a, b, X = sympy.symbols("a, b, x")
f = sympy.Function("f")
x = a, (a+b)/2, b # simpson's rule
# weight factors
w = [sympy.symbols("w_%d" % i) 
     for i in range(len(x))] 
w

$\displaystyle \left[ w_{0}, \ w_{1}, \ w_{2}\right]$

# quadrature rule
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)}$

# using lambda to create symbolics for basis functions
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]$

# solves for weight factors
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}$

Numerical Integration w/ SciPy

  • Two types of quadrature routines:

  1. accepts integrand as Python function using a Gaussian quadrature - quad(), quadrature(), fixed_quad())

  2. Newton-Cotes methods: accepts arrays of integrand samples at given points - trapz(), simps(), romb())

Simple integration example

def f(x):
    return np.exp(-x**2)
# integrate.quad(function, upper limit, lower limit)
# return estimated integral & abs error
# error tolerances can be set using epsabs, epsrel keywords
val, err = integrate.quad(f, -1, 1)
val, err

$\displaystyle \left( 1.49364826562485, \ 1.65828269518814 \cdot 10^{-14}\right)$

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, err

$\displaystyle \left( 1.27630683510222, \ 1.41698523481695 \cdot 10^{-14}\right)$

Reshuffle arguments

from scipy.special import jv # zeroth-order bessel function
val, err = integrate.quad(lambda x: jv(0, x), 0, 5)
val, err

$\displaystyle \left( 0.715311917784768, \ 2.47260738289741 \cdot 10^{-14}\right)$

Infinite limits

f = lambda x: np.exp(-x**2)
val, err = integrate.quad(f, -np.inf, np.inf)
val, err

$\displaystyle \left( 1.77245385090552, \ 1.42026367809449 \cdot 10^{-8}\right)$

Singularities

  • quadrature() & fixed_quad() functions only support finite integration limits.
f = lambda x: 1/np.sqrt(abs(x))
a, b = -1, 1
# this integral diverges at x=0. quad() can't handle it.
integrate.quad(f, a, b)
/tmp/ipykernel_65767/525825373.py:1: RuntimeWarning: divide by zero encountered in scalar divide
  f = lambda x: 1/np.sqrt(abs(x))
/tmp/ipykernel_65767/14331687.py:2: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  integrate.quad(f, a, b)

$\displaystyle \left( \text{NaN}, \ \text{NaN}\right)$

# to avoid problem points, use points keyword
integrate.quad(f, a, b, points=[0])

$\displaystyle \left( 3.99999999999998, \ 5.6843418860808 \cdot 10^{-14}\right)$

# visualize to see what's going on
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")

png

Tabulated integrand

  • use case: integrand that is specified only at predetermined points. (ie, not using a Python function that can be evaluated at any point.)
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")

png

# trapezoid method
val_trapz = integrate.trapz(y, x)
val_trapz

$\displaystyle 1.88082171605085$

# simpson's rule
val_simps = integrate.simps(y, x)
val_simps

$\displaystyle 1.88366510244871$

# find error estimates by computing integral & comparising to above estimate
val_exact = 2.0/3.0 * (b-a)**(3.0/2.0)

val_exact, (val_exact-val_trapz), (val_exact-val_simps)

$\displaystyle \left( 1.88561808316413, \ 0.00479636711327625, \ 0.00195298071541172\right)$

  • We can't ask trapz or simps to find more accurate answers unless more samples are given, or a higher-order method is used.
  • romb() function helps -- implementation of Romberg method
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.000378798422913107$

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

$\displaystyle 0.000448485554158218$

Multiple integration

  • dblquad and tplquad functions available in SciPy.
def f(x):
    return np.exp(-x**2)
%time integrate.quad(f, a, b)
CPU times: user 37 µs, sys: 21 µs, total: 58 µs
Wall time: 60.1 µs

$\displaystyle \left( 0.882081390762422, \ 9.7930706961782 \cdot 10^{-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.557746285351034, \ 8.29137438153541 \cdot 10^{-15}\right)$

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

$\displaystyle \left( 0.557746285351034, \ 8.29137438153541 \cdot 10^{-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")

png

# because g&h are functions,
# we can find integrals with x-dependent integration limits along y.
integrate.dblquad(
    f, 
    0, 1, 
    lambda x: -1+x, 
    lambda x: 1-x)

$\displaystyle \left( 0.732093100000809, \ 1.6564972931774 \cdot 10^{-14}\right)$

triple integrals

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.416538385886638, \ 8.29133528731442 \cdot 10^{-15}\right)$

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

$\displaystyle \left( 0.416538385886638, \ 8.29133528731442 \cdot 10^{-15}\right)$

Arbitrary number of integrations

  • Computational complexity grows very quickly.
def f(*args):
    return  np.exp(-np.sum(np.array(args)**2))
%time integrate.nquad(f, [(0,1)] * 1)
CPU times: user 110 µs, sys: 60 µs, total: 170 µs
Wall time: 173 µs

$\displaystyle \left( 0.746824132812427, \ 8.29141347594073 \cdot 10^{-15}\right)$

%time integrate.nquad(f, [(0,1)] * 2)
CPU times: user 1.55 ms, sys: 0 ns, total: 1.55 ms
Wall time: 1.55 ms

$\displaystyle \left( 0.557746285351034, \ 8.29137438153541 \cdot 10^{-15}\right)$

%time integrate.nquad(f, [(0,1)] * 3)
CPU times: user 40.2 ms, sys: 134 µs, total: 40.3 ms
Wall time: 40.3 ms

$\displaystyle \left( 0.416538385886638, \ 8.29133528731442 \cdot 10^{-15}\right)$

%time integrate.nquad(f, [(0,1)] * 4)
CPU times: user 593 ms, sys: 1.2 ms, total: 594 ms
Wall time: 593 ms

$\displaystyle \left( 0.311080918822877, \ 8.29129619327777 \cdot 10^{-15}\right)$

%time integrate.nquad(f, [(0,1)] * 5)
CPU times: user 11.4 s, sys: 2.3 ms, total: 11.4 s
Wall time: 11.4 s

$\displaystyle \left( 0.232322737434388, \ 8.29125709942545 \cdot 10^{-15}\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));

png

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

$\displaystyle \pi$

Sympy mpath library

  • alternate method of numerical quadrature using arbitrary-precisin math
  • much slower than floating-point math, but better precision than SciPy quadrature funcs.
import mpmath as mpmath
mpmath.mp.dps = 75 # number of digits of precision
# need mpmath-compatible python function as integrand.
f_mpmath = sympy.lambdify(x, f, 'mpmath')
# compute integral & display result
val = mpmath.quad(f_mpmath, (a, b))
sympy.sympify(val)

$\displaystyle 3.14159265358979323846264338327950288419716939937510582097494459230781640629$

# compare result with known value. error sufficiently small?
sympy.N(val_sym, mpmath.mp.dps+1) - val

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

%timeit mpmath.quad(f_mpmath, [a, b])
2.27 ms ± 501 µ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)
242 µs ± 54.3 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Line integrals

  • 1st argument = integrand (symPy expression)
  • 2nd argument = symPy Curve object
  • 3rd argument = list of integration variables
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

  • integral transform: accepts a function as input, returns another function.
  • best when computed symbolically
  • two examples: laplace (differential equation into algebraic equation) and fourier (time domain to frequency domain).

Laplace transforms

s    = sympy.symbols("s")
a, t = sympy.symbols("a, t", positive=True)
f    = sympy.sin(a*t)
# compute transform
# returns tuple: result transform, value 'A' from convergence condition, any add'l conditions.

sympy.laplace_transform(f, t, s)

$\displaystyle \left( \frac{a}{a^{2} + s^{2}}, \ 0, \ \text{True}\right)$

# if you want just the result function:
F = sympy.laplace_transform(f, t, s, noconds=True); F

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

# inverse transform (swap 's','t' arguments)
sympy.inverse_laplace_transform(F, s, t, noconds=True)

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

# sympy can find transforms for many elementary functions automatically
# ex: laplace transforms of polynomials

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

# computing general results with arbitrary integer exponents
n = sympy.symbols(
    "n", integer=True, positive=True)

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

$\displaystyle s^{- n - 1} \Gamma\left(n + 1\right)$

# laplace xform of composite expression
sympy.laplace_transform(
    (1 - a*t) * sympy.exp(-a*t), t, s, noconds=True)

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

Fourier Transforms

  • available functions: fourier_transform(), inverse_fourier_transform()
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}}$