Use scipy.integrate.quad
to evaluate the following definite integrals (which can also be expressed in closed form over the range given but are awkward).
(a) $$ \int_0^1 \frac{x^4(1-x)^4}{1+x^2}\;\mathrm{d}x. $$ (Compare with $22/7 - \pi$).
(b) The following integral appears in the Debye theory of the heat capacity of crystals at low temperature $$ \int_0^\infty \frac{x^3}{e^x-1}\;\mathrm{d}x $$ (Compare with $\pi^4/15$).
(c) The integral sometimes known as the Sophomore's dream: $$ \int_0^1 x^{-x}\;\mathrm{d}x $$ (Compare the value you obtain from the summation $\sum_{n=1}^\infty n^{-n}$).
(d) $$ \int_0^1 [\ln(1/x)]^p\;\mathrm{d}x $$ (Compare with $p!$ for integer $0 \le p \le 10$).
(e) $$ \int_0^{2\pi} e^{z\cos\theta}\;\mathrm{d}\theta $$ (Compare with $2\pi I_0(z)$, where $I_0(z)$ is a modified Bessel function of the first kind, for $0 \le z \le 2$).
In the following we assume the following imports:
In [x]: import numpy as np
In [x]: from scipy.integrate import quad
(a)
In [x]: f1 = lambda x: x**4 * (1 - x)**4/(1 + x**2)
In [x]: quad(f1, 0, 1)
Out[x]: (0.0012644892673496185, 1.1126990906558069e-14)
In [x]: 22/7 - np.pi
Out[x]: 0.0012644892673496777
(b)
In [x]: f2 = lambda x: x**3/(np.exp(x) - 1)
In [x]: quad(f2, 0, np.inf)
Out[x]: (6.49393940226683, 2.628470028924825e-09)
In [x]: np.pi**4 / 15
Out[x]: 6.493939402266828
(c)
In [x]: f3 = lambda x: x**-x
In [x]: quad(f3, 0, 1)
Out[x]: (1.2912859970626633, 3.668398917966442e-11)
In [x]: n, I, TOL = 0, 0, 1.e-16
In [x]: while True:
...: Iold = I
...: n += 1
...: I += n**-n
...: if I-Iold < TOL:
...: break
...:
In [x]: I
Out[x]: 1.2912859970626636
(d)
In [x]: from scipy.misc import factorial
In [x]: f4 = lambda x, p: np.log(1/x)**p
In [x]: for p in range(10):
...: print(quad(f4, 0, 1, args=(p,))[0], factorial(p))
...:
1.0 1.0
0.9999999999999999 1.0
1.9999999999999991 2.0
6.000000000000064 6.0
24.000000000000014 24.0
119.9999999999327 120.0
719.9999999989705 720.0
5039.99999945767 5040.0
40320.00000363255 40320.0
362880.00027390465 362880.0
(e) Here we find the maximum deviation from the Bessel function $I_0(z)$ for 100 values of $z$ between 0 and 2:
In [x]: from scipy.special import i0
In [x]: z = np.linspace(0,2,100)
In [x]: y1 = i0(z)
In [x]: f5 = lambda theta, z: np.exp(z*np.cos(theta))
In [x]: y2 = np.array([quad(f5, 0, 2*np.pi, args=(zz,))[0] for zz in z])
In [x]: y2 /= 2 * np.pi
In [x]: np.max(abs(y2-y1))
Out[x]: 3.4796610037801656e-12