5  Integration

Exercise 5.1 - Rectangle rule

The rectangle rule is based on the definition of the integral in the Riemann sense. The first step is to split the interval \([a,b]\) into \(N\) intervals \([x_n,x_{n+1}]\) of the same size \(\delta=\frac{b-a}{N}\), i.e., \(x_n = a+n\delta\) for \(n\in\{0,1,\dots,N-1\}\). The second step consists in assuming that the function \(f\) is constant on each interval \([x_n,x_{n+1}]\), thus to make the approximation:

\[ J_n = \int_{x_n}^{x_{n+1}} f(x) \, \mathrm{d} x \approx \delta f(\tilde{x}_n) \,, \]

for \(\tilde{x}_n\) a certain value to choose in the interval \([x_n,x_{n+1}]\). The choice of \(\tilde{x}_n\) can, for example, be done by \(\tilde{x}_n = x_n + \alpha\delta\) with \(\alpha\in[0,1]\). Finally, the approximation of \(J\) is given by the sum of the approximations of \(J_n\),

\[ \tilde{J} = \sum_{n=0}^{N-1} \delta f(\tilde{x}_n)\,. \]

Assuming that \(f\in C^1([a,b])\), it is then possible to show that the rectangle rule converges and that its speed of convergence is of order one. A numerical method is of order \(k\) if the error between the numerical approximation and the exact result is of order \(N^{-k}\).

a. Choose a continuous function \(f:[a,b]\to\mathbb{R}\) and define the corresponding Python function f(x). To test the code, it is wise to choose a function \(f\) whose integral can be easily computed by hand.

b. Write a function rectangles(f,a,b,N) that returns the approximation of the integral \(J\) by the rectangle rule, for example, by choosing \(\tilde{x}_n=x_n\), i.e., the left edge of the interval \([x_n,x_{n+1}]\).

c. Modify the previous function so that it takes an optional parameter alpha determining the choice of parameter \(\alpha\in[0,1]\).

d. Write a function plot_rectangles(f,a,b,N,alpha=0.5) that graphically represents the approximation by the rectangle rule.

e. Determine empirically the speed of convergence of the rectangle rule as a function of \(N\).

f. Determine analytically the convergence of the rectangle rule. What are the necessary assumptions on \(f\)?

Exercise 5.2 - Trapezoidal rule

The trapezoidal rule is based on a linear approximation on each interval \([x_n,x_{n+1}]\), more specifically:

\[ J_n = \int_{x_n}^{x_{n+1}} f(x) \, \mathrm{d} x \approx \delta \frac{f(x_n) + f(x_{n+1})}{2} \,. \]

a. Write a Python function trapezes(f,a,b,N) that returns the approximation of the integral \(J\) by the trapezoidal rule. Test the function trapezes(f,a,b,N) for different functions \(f\).

b. Is your implementation of the function trapezes(f,a,b,N) optimal in terms of the number of evaluations of \(f\) performed compared to the number of evaluations needed? An optimal implementation of the function trapezes(f,a,b,N) should perform \(N+1\) evaluations of \(f\).

c. Determine empirically the speed of convergence of the trapezoidal rule as a function of \(N\).

d. ! Determine analytically the convergence of the trapezoidal rule. What are the necessary assumptions on \(f\)?

Exercise 5.3 - Monte Carlo method

The Monte Carlo method (named after casinos, not a person) is a probabilistic approach to approximate the value of an integral. The basic idea is that the integral \(J\) can be seen as the expectation of a uniform random variable \(X\) on the interval \([a,b]\):

\[ J = \int_a^b f(x) \, \mathrm{d} x = (b-a)\mathbb{E}(f(X)) \,. \]

By the law of large numbers, this expectation can be approximated by the empirical mean:

\[ \tilde{J} = \frac{b-a}{N}\sum_{i=0}^{N-1} f(x_i) \,, \]

where \(x_i\) are drawn randomly in the interval \([a,b]\) with a uniform probability distribution.

a. Write a function montecarlo(f,a,b,N) that determines an approximation \(\tilde{J}\) of \(J\) by the Monte Carlo method.

b. Modify the previous function, so that it returns in addition to the mean \(\tilde{J}\) also the empirical variance:

\[ \tilde{V} = \frac{(b-a)^2}{N}\sum_{i=0}^{N-1} \left(f(x_i)-\frac{\tilde{J}}{b-a}\right)^2 \,. \]

c. Study empirically the convergence of the Monte Carlo method as a function of \(N\) by making for each value of \(N\) a statistic on \(k\) executions. More precisely, this consists in making \(k\) evaluations of \(\tilde{J}\) by the Monte Carlo method and to calculate the mean and the variance of the \(k\) results obtained.

d. Determine analytically the convergence of the Monte Carlo method. What are the necessary assumptions on \(f\)?

Exercise 5.4 - Simpson’s rule !

Simpson’s rule consists in approximating the function \(f\) on each interval \([x_n,x_{n+1}]\) by a polynomial of degree 2. The most natural choice is the polynomial \(p_n\) of degree 2 passing through the points \((x_n,f(x_n))\), \((\frac{x_n+x_{n+1}}{2},f(\frac{x_n+x_{n+1}}{2}))\), and \((x_{n+1},f(x_{n+1}))\).

a. Determine the explicit form of the polynomial \(p_n\).

b. Compute the approximation given by \(J_n \approx \int_{x_n}^{x_{n+1}} p_n(x)\,\mathrm{d} x \,.\)

c. Simplify by hand the sum \(\tilde{J}\) of the approximations of \(J_n\).

d. Write a function simpson(f,a,b,N) to approximate \(J\) with Simpson’s rule.

e. Compare the accuracy, i.e., the convergence speed of the rectangle, trapezoid, and Simpson methods as a function of \(N\) for a smooth function and the function \(f(x)=\sqrt{1-x^2}\) on \([0,1]\) (whose integral is \(\frac{\pi}{4}\)).

f. If not already done, propose a parallel implementation of Simpson’s rule using NumPy indexing.

Exercise 5.5 - Integration with SciPy !!

The above and other integration methods are defined in the integrate module of SciPy. This module allows in particular to handle more complicated cases: singular, generalized, or multidimensional integrals.

a. Define a function E(n,x) computing numerically the following integral:

\[ E_n(x) = \int_1^\infty \frac{e^{-xt}}{t^n} \mathrm{d} t \,. \]

b. Determine an approximation of the double integral:

\[ I = \int_{0}^{\pi} \left(\int_{0}^{y} x \sin(xy) \,\mathrm{d} x \right) \mathrm{d} y \,. \]