Numerical integration of some awkward integrals

Q8.2.2Numerical integration of some awkward integrals

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$ ).

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

Show solution

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