Le but est d'obtenir une approximation d'une intégrale définie du type $$ J = \int_a^b f(x) \, \mathrm{d} x $$ pour une certaine fonction \( f:[a,b] \to \mathbb{R} \) trop compliquée pour a priori déterminer la valeur de \( J \) à la main. Des méthodes d'approximations déterministes et probabilistes seront introduites pour obtenir une approximation \( \tilde{J} \) de \( J \).
La méthode des rectangles est basée sur la définition de l'intégrale au sens de Riemann. La première étape est de découper l'intervalle \( [a,b] \) en \( N \) intervalles \( [x_n,x_{n+1}] \) de même taille \( \delta=\frac{b-a}{N} \), i.e. \( x_n = a+n\delta \) pour \( n\in\{0,1,\dots,N-1\} \). La seconde étape consiste à supposer que la fonction \( f \) est constante sur chaque intervalle \( [x_n,x_{n+1}] \), donc à faire l'approximation $$ J_n = \int_{x_n}^{x_{n+1}} f(x) \, \mathrm{d} x \approx \delta f(\tilde{x}_n) $$ pour \( \tilde{x}_n \) une certaine valeur à choisir dans l'intervalle \( [x_n,x_{n+1}] \). Le choix de \( \tilde{x}_n \) peut par exemple être fait par \( \tilde{x}_n = x_n + \alpha\delta \) avec \( \alpha\in[0,1] \). Finalement l'approximation de \( J \) est donnée par la somme des approximations de \( J_n \), $$ \tilde{J} = \sum_{n=0}^{N-1} \delta f(\tilde{x}_n)\,. $$
En supposant que \( f\in C^1([a,b]) \), alors il est possible de montrer que la méthode des rectangles converge et que sa vitesse de convergence est d'ordre un. Une méthode numérique est dite d'ordre \( k \) si l'erreur entre l'approximation numérique et le résultat exact est de l'ordre de \( N^{-k} \).
Pour chaque valeur de \( n \), et \( x\in[x_{n},x_{n+1}] \), par le théorème des accroissements finis, il existe \( c_{n} \) tel que: $$ f(x)-f(\tilde{x}_{n})=(x-\tilde{x}_{n})f^{\prime}(c_{n})\,. $$ Si \( f^{\prime} \) est continue sur \( [a,b] \), alors $$ \sup_{x\in[a,b]}|f^{\prime}(x)|\leq M\,, $$ et donc $$ |f(x)-f(\tilde{x}_{n})|\leq M|x-\tilde{x}_{n}|\leq M\delta\leq\frac{M(b-a)}{N}\,. $$ Ainsi $$ \begin{align*} E_{N} & = \left|\tilde{J}-J\right| = \left|\sum_{n=0}^{N-1}\left(\int_{x_{n}}^{x_{n+1}}f(x)\,\mathrm{d} x-\delta f(\tilde{x}_{n})\right)\right|\\ & \leq\sum_{n=0}^{N-1}\left|\int_{x_{n}}^{x_{n+1}}\left(f(x)-f(\tilde{x}_{n})\right)\,\mathrm{d} x\right|\leq\sum_{n=0}^{N-1}\int_{x_{n}}^{x_{n+1}}|f(x)-f(\tilde{x}_{n})|\,\mathrm{d} x\\ & \leq\sum_{n=0}^{N-1}\int_{x_{n}}^{x_{n+1}}\frac{M(b-a)}{N}\,\mathrm{d} x\leq\frac{M(b-a)^{2}}{N}\,. \end{align*} $$ Par conséquent si \( f\in C^1([a,b]) \), alors la méthode des rectangles est d'ordre un.
a) Choisir une fonction continue \( f:[a,b]\to\mathbb{R} \) et définir la fonction Python f(x) correspondante. Pour tester le code, il est judicieux de choisir une fonction \( f \) dont l'intégrale peut être facilement calculable à la main.
La liste des fonctions mathématiques de base disponibles en Python dans le module math est disponible ici. À noter que Numpy définit également des fonctions mathématiques, voir la documentation ici.
b) Écrire une fonction rectangles(f,a,b,N) qui retourne l'approximation de l'intégrale \( J \) par la méthode des rectangles par exemple en choisissant \( \tilde{x}_n=x_n \), i.e. le bord gauche de l'intervalle \( [x_n,x_{n+1}] \).
Il n'est pas nécessaire de stoker toutes les valeurs des approximations de \( J_n \), mais il est possible d'incrémenter une variable pour chaque approximation de \( J_n \).
c) Modifier la fonction précédente pour que celle-ci prenne un paramètre optionnel alpha déterminant le choix du paramètre \( \alpha\in[0,1] \).
d) Écrire une fonction plot_rectangles(f,a,b,N,alpha=0.5) qui représente graphiquement l'approximation par la méthode des rectangles.
e) Déterminer empiriquement la vitesse de convergence de la méthode des rectangles en fonction de \( N \).
La méthode des trapèzes est basée sur une approximation linéaire sur chaque intervalle \( [x_n,x_{n+1}] \), plus spécifiquement: $$ 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) Écrire une fonction Python trapezes(f,a,b,N) qui retourne l'approximation de l'intégrale \( J \) par la méthode des trapèzes. Tester la fonction trapezes(f,a,b,N) pour différentes fonctions \( f \).
b) L'implémentation de votre fonction trapezes(f,a,b,N) est-elle optimale quant au nombre d'évaluations de \( f \) effectuées par rapport au nombre d'évaluations nécessaires ? Une implémentation optimale de la fonction trapezes(f,a,b,N) devrait effectuer \( N+1 \) évaluations de \( f \).
c) Déterminer empiriquement la vitesse de convergence de la méthode des trapèzes en fonction de \( N \).
d) ! Déterminer analytiquement la convergence de la méthode des trapèzes. Quelles sont les hypothèses nécessaires sur \( f \) ?
La méthode de Monte-Carlo (du nom des casinos, pas d'une personne) est une approche probabiliste permettant d'approximer la valeur d'une intégrale. L'idée de base est que l'intégrale \( J \) peut être vue comme l'espérance d'une variable aléatoire uniforme \( X \) sur l'intervalle \( [a,b] \): $$ J = \int_a^b f(x) \, \mathrm{d} x = (b-a)\mathbb{E}(f(X)) \,. $$ Par la loi des grands nombres cette espérance peut être approximée par la moyenne empirique: $$ \tilde{J} = \frac{b-a}{N}\sum_{i=0}^{N-1} f(x_i) \,, $$ où les \( x_i \) sont tirés aléatoirement dans l'intervalle \( [a,b] \) avec une loi de probabilité uniforme.
Il est possible de montrer que \( \tilde{J} \) converge vers \( J \) comme \( N^{-1/2} \) et cela indépendamment de la dimension et de la régularité de \( f \).
Selon le théorème central limite, si \( Y_i \) est une suite de variables aléatoires indépendantes d'espérance \( \mu \) et de variance \( \sigma^2 \), alors la variable aléatoire: $$ S_N = \frac{1}{N}\sum_{i=0}^{N-1} Y_i \,, $$ a une espérance \( \mu \) et une variance: $$ \mathrm{Var}(S_N) = \frac{\sigma^2}{N} \,. $$ En prenant \( Y_i=f(X_i) \) avec \( X_i \) une suite de variables aléatoires indépendantes uniformément distribuées sur \( [a,b] \), alors l'espérance de \( Y_i \) est la moyenne de \( f \) et donc l'espérance de \( S_N \) est donnée par: $$ \mathbb{E}(S_N) = \frac{1}{b-a} \int_a^b f(x) \, \mathrm{d} x \,. $$ La variance de \( Y_i \) est également celle de \( f \), \( \sigma^2 = \mathrm{Var}(f(X)) \) et donc la variance de \( S_N \) est: $$ \mathrm{Var}(S_N) = \frac{\mathrm{Var}(f(X))}{N} \,. $$ Par conséquent, cela montre que \( \tilde{J} \) converge vers \( J \) comme \( N^{-1/2} \) vu que la variance est proportionnelle à \( N^{-1} \).
A remarquer que pour établir ce résultat, aucune condition de régularité sur \( f \) n'est nécessaire, de l'intégrabilité suffit.
En pratique la variance de \( f \) peut être estimée par la variance empirique.
a) Écrire une fonction montecarlo(f,a,b,N) qui détermine une approximation \( \tilde{J} \) de \( J \) par la méthode de Monte-Carlo.
Pour générer un vecteur de nombres aléatoires, le sous-module random de Numpy peut être utile, voir la documentation ici.
b) Modifier la fonction précédente, pour qu'elle retourne en plus de la moyenne \( \tilde{J} \) également la variance empirique: $$ \tilde{V} = \frac{(b-a)^2}{N}\sum_{i=0}^{N-1} \left(f(x_i)-\frac{\tilde{J}}{b-a}\right)^2 \,. $$
c) Étudier empiriquement la convergence de la méthode de Monte-Carlo en fonction de \( N \) en faisant pour chaque valeur de \( N \) une statistique sur \( k \) exécutions. Plus précisément cela consiste à faire \( k \) évaluations de \( \tilde{J} \) par la méthode de Monte-Carlo et de calculer la moyenne et la variance des \( k \) résultats obtenus.
La méthode de Simpson consiste à approximer la fonction \( f \) sur chaque intervalle \( [x_n,x_{n+1}] \) par un polynôme de degré deux. Le choix le plus naturel est le polynôme \( p_n \) de degré deux passant par les points \( (x_n,f(x_n)) \), \( (\frac{x_n+x_{n+1}}{2},f(\frac{x_n+x_{n+1}}{2})) \) et \( (x_{n+1},f(x_{n+1})) \).
a) Déterminer la forme explicite du polynôme \( p_n \).
Le polynôme \( L(x) = \frac{(x-c)(x-b)}{(a-c)(a-b)} \) prend la valeur 1 en \( x=a \) et la valeur 0 en \( x=b \) et \( x=c \). Faire une combinaison linéaire de trois polynômes de ce type.
b) Calculer l'approximation donnée par \( J_n \approx \int_{x_n}^{x_{n+1}} p_n(x)\,\mathrm{d} x \,. \)
Il est possible de calculer cette intégrale à la main ou bien de le faire avec le module Sympy, voir la documentation ici.
c) Simplifier à la main la somme \( \tilde{J} \) des approximations de \( J_n \).
Le résultat est: $$ \tilde{J} = \frac{\delta}{3}\left[\frac{f(b)-f(a)}{2}+\sum_{n=0}^{N-1}\left(f(x_{n})+2f\left(\frac{x_{n}+x_{n+1}}{2}\right)\right)\right] \,. $$
d) Écrire une fonction simpson(f,a,b,N) permettant de calculer une approximation de \( J \) avec la méthode de Simpson.
e) Comparer la précision des méthodes des rectangles, des trapèzes et de Simpson en fonction de \( N \) pour une fonction lisse et la fonction \( f(x)=\sqrt{1-x^2} \) sur \( [0,1] \) (dont l'intégrale vaut \( \frac{\pi}{4} \)).
f) Si ce n'est pas déjà fait proposer une implémentation parallèle de la méthode de Simpson en utilisant l'indexage Numpy.
Les méthodes d'intégrations précédentes et d'autres sont définies dans le module integrate de Scipy. Ce module permet en particulier de traiter des cas plus compliqués: intégrales singulières, généralisées ou en plusieurs dimensions.
a) Définir une fonction E(n,x) calculant numériquement l'intégrale suivante: $$ E_n(x) = \int_1^\infty \frac{e^{-xt}}{t^n} \mathrm{d} t \,. $$
Lire la documentation du sous-module integrate de Scipy disponible ici.
b) Déterminer une approximation de l'intégrale double: $$ I = \int_{0}^{\pi} \left(\int_{0}^{y} x \sin(xy) \,\mathrm{d} x \right) \mathrm{d} y \,. $$