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()
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
phi = [sympy.Lambda(X, X**n) for n in range(len(x))]
phi
eqs = [q_rule.subs(f, phi[n]) - sympy.integrate(phi[n](X), (X, a, b)) for n in range(len(phi))]
eqs
w_sol = sympy.solve(eqs, w)
w_sol
q_rule.subs(w_sol).simplify()
integrate
¶def f(x):
return np.exp(-x**2)
val, err = integrate.quad(f, -1, 1)
val
err
val, err = integrate.quadrature(f, -1, 1)
val
err
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
from scipy.special import jv
val, err = integrate.quad(lambda x: jv(0, x), 0, 5)
val
err
f = lambda x: np.exp(-x**2)
val, err = integrate.quad(f, -np.inf, np.inf)
val
err
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.
integrate.quad(f, a, b, points=[0])
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")
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
val_simps = integrate.simps(y, x)
val_simps
val_exact = 2.0/3.0 * (b-a)**(3.0/2.0)
val_exact
val_exact - val_trapz
val_exact - val_simps
x = np.linspace(a, b, 1 + 2**6)
len(x)
y = f(x)
val_exact - integrate.romb(y, dx=(x[1]-x[0]))
val_exact - integrate.simps(y, dx=x[1]-x[0])
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
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)
integrate.dblquad(lambda x, y: np.exp(-x**2-y**2), 0, 1, lambda x: 0, lambda x: 1)
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)
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)
integrate.nquad(f, [(0, 1), (0, 1), (0, 1)])
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
%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
%time integrate.nquad(f, [(0,1)] * 3)
CPU times: user 169 ms, sys: 3.38 ms, total: 172 ms Wall time: 171 ms
%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
%time integrate.nquad(f, [(0,1)] * 5)
CPU times: user 1min 9s, sys: 237 ms, total: 1min 9s Wall time: 1min 9s
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
%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
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
mpmath.mp.dps = 75
f_mpmath = sympy.lambdify(x, f, 'mpmath')
val = mpmath.quad(f_mpmath, (a, b))
sympy.sympify(val)
sympy.N(val_sym, mpmath.mp.dps+1) - val
%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)
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)
integrate.tplquad(f3, 0, 1, lambda x : 0, lambda x : 1, lambda x, y : 0, lambda x, y : 1)
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
f2_numpy = sympy.lambdify((x, y), f2, 'numpy')
integrate.dblquad(f2_numpy, 0, 1, lambda x: 0, lambda x: 1)
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)
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)
%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
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])
sympy.line_integrate(x**2 * y**2, C, [x, y])
s = sympy.symbols("s")
a, t = sympy.symbols("a, t", positive=True)
f = sympy.sin(a*t)
sympy.laplace_transform(f, t, s)
F = sympy.laplace_transform(f, t, s, noconds=True)
F
sympy.inverse_laplace_transform(F, s, t, noconds=True)
[sympy.laplace_transform(f, t, s, noconds=True) for f in [t, t**2, t**3, t**4]]
n = sympy.symbols("n", integer=True, positive=True)
sympy.laplace_transform(t**n, t, s, noconds=True)
sympy.laplace_transform((1 - a*t) * sympy.exp(-a*t), t, s, noconds=True)
w = sympy.symbols("omega")
f = sympy.exp(-a*t**2)
F = sympy.fourier_transform(f, t, w)
F
sympy.inverse_fourier_transform(F, w, t)
sympy.fourier_transform(sympy.cos(t), t, w) # not good
%reload_ext version_information
%version_information numpy, matplotlib, scipy, sympy, mpmath, skmonaco
Software | Version |
---|---|
Python | 3.6.8 64bit [GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)] |
IPython | 7.5.0 |
OS | Darwin 18.2.0 x86_64 i386 64bit |
numpy | 1.16.3 |
matplotlib | 3.0.3 |
scipy | 1.2.1 |
sympy | 1.4 |
mpmath | 1.1.0 |
skmonaco | 0.2.1 |
Mon May 06 14:55:36 2019 JST |