
\chapterimage{romberg}
\chaptertext{Il  est très  rare de  pouvoir  déterminer la
dérivée  ou  une  primitive  formelle d'une  fonction  donnée.  Il  faut
recourir  à  des  méthodes  numériques dont  l'utilisation  entraîne  de
nouveaux  problèmes:  il  faut   pouvoir  contrôler  l'approximation  de
manière théorique  et ne pas  oublier que cela  sera fait la  plupart du
temps par une machine qui elle aussi apporte son lot d'imprécision.

Werner \textsc{Romberg} (1909 - 2003) né à Berlin, s'installe en Norvège
en 1938.  Dix-sept ans  plus tard,  il publie un  article se  basant sur
l'extrapolation  de  Lewis Richardson  (1881  -  1953) et  qui  accélère
grandement l'efficacité de la méthode  dite des trapèzes utilisée depuis
longtemps  pour calculer  des approximations  d'intégrales. C'est  cette
méthode que nous étudierons aujourd'hui  en la comparant aux méthodes de
\textsc{Newton-Cotes},  \textsc{Gauss-Legendre},  \textsc{Kuncir}  après
une introduction sur les méthodes de dérivation numérique.
}

\chapter{Intégration numérique}



\section{Méthode de quadrature simplifiée déterministe}

On étudiera  des formules d'approximation  reposant sur des  formules du
type:

$$\displaystyle\int_{a}^{b}f(x)\;\mathrm{d}x=
\displaystyle\sum_{i=1}^{n}w_{i,n}f(x_i)+E_n(f)$$

où  les $w_i$  sont les  \textit{poids}, les  $x_i$ sont  les pivots  et
$E(f)$ l'erreur commise.

Une méthode de quadrature est associée à un ordre:

\begin{definition}
Une méthode de  quadrature est d'ordre $m$ si elle  est exacte pour tous
les polynômes de $\bbr_m[X]$.
\end{definition}


Nous admettrons  le théorème  suivant qui donne  une CNS  de convergence
d'une méthode de quadrature:

\begin{theoreme}
  Pour toute fonction continue sur $]a,b[$, $
  \displaystyle\lim_{n\to +\infty}E_n(f)=0$ si, et seulement si, 
  \begin{itemize}
  \item                $(\exists               M\in\bbr^\star_+)(\forall
        n\in\bbn)(\displaystyle\sum_{i=0}^{n}|w_{i,n}|\leqslant M)$
      \item $(\forall k\in \bbn)(
        \displaystyle\lim_{n\to +\infty}E_n(x\mapsto x^k)=0$
)
  \end{itemize}
\end{theoreme}


On  peut se  contenter,  pour la  condition  suffisante, d'avoir  chaque
$w_{i,n}$ positif et $f$ continue par morceaux.

Ce qui nous  intéressera au plus haut point est  le théorème suivant qui
donne une estimation de l'erreur commise:



\begin{theoreme}
  On    suppose     que    la    méthode    est     d'ordre    $m$    et
  $f\in\CR^{(m+1)}([a,b],\bbr)$. Alors


$$E(f)=\frac{1}{m!}\displaystyle\int_{a}^{b}K_m(x)f^{(m+1)}\;\mathrm{d}x$$

avec  $K_n(x)=E\left(x\mapsto  \max  \left((x-t)^n,0\right)\right)$.  La
famille   $(K_n)_{n\in\bbn}$  s'appelle   la  famille   des  noyaux   de
\textsc{Peano} associée à la méthode considérée.

\end{theoreme}


On utilise la  formule de \textsc{Taylor} intégrale puis le  fait que la
méthode est d'ordre $m$.

On en déduit le corollaire suivant:

\begin{theoreme}
  Si $K_m$ est de signe constant sur $[a,b]$:

$$(\exists           \xi           \in           ]a,b[)           (E(f)=
\frac{1}{(m+1)!}f^{(m+1)}(\xi)E(x\mapsto x^{m+1}))$$
\end{theoreme}


On utilise la deuxième formule de la moyenne.






\section{Méthodes de \textsc{Newton-Cotes} composées}



Pour donner une valeur approchée de l'intégrale $I$, une première idée est de remplacer la
fonction $f$ par un polynôme $P$ qui interpole $f$ en plusieurs points.
Cependant, dès qu'on augmente le nombre de pivots, on risque fort
de se retrouver confronté au phénomène de \textsc{Runge}.

On modifie alors l'idée de départ ainsi : on commence par subdiviser régulièrement
l'intervalle d'intégration en $n$ sous-intervalles $[x_j,x_{j+1}]$ où
\[
x_0 = a,\quad x_n = b,\quad \text{et} \quad
∀j∈⟦0, n⟧,\quad x_j = a + j h,\quad
\text{ avec } h = \frac{b-a}{n}
\]
Ensuite, on remplace $f$ par un polynôme $P_j$ qui interpole $f$ sur chacun des
petits segments $[x_j,x_{j+1}]$.
\begin{itemize}
    \item Si $P_j$ interpole $f$ au point $x_j$, alors le
        graphe de $P_j$ est une droite horizontale : c'est la
        méthode des rectangles.
    \item Si $P_j$ interpole $f$ aux points $x_j$ et $x_{j+1}$, alors le
        graphe de $P_j$ est une droite affine : c'est la méthode des trapèzes.
    \item Si $P_j$ interpole $f$ aux points $x_j$, $\frac{x_j+x_{j+1}}{2}$ et $x_{j+1}$, alors le
        graphe de $P_j$ est une parabole: c'est la méthode de \textsc{Simpson} $1/3$.
\end{itemize}
Enfin, on somme les intégrales de chaque polynôme $P_j$ sur $[x_j,x_{j+1}]$, et
on obtient une valeur approchée de $I$ :
\[
I = ∑_{j=0}^{n-1} ∫_{x_j}^{x_{j+1}} f(x) \dd x
= ∑_{j=0}^{n-1} \left(∫_{x_j}^{x_{j+1}} P_j(x) \dd x + E_j\right)
≈ ∑_{j=0}^{n-1} ∫_{x_j}^{x_{j+1}} P_j(x) \dd x
\]


On se place sur l'intervalle $[x_j,x_{j+1}]$ et on note
$ξ_j=\tfrac{x_j+x_{j+1}}{2}$ ; alors d'après les formules des différences
divisées (cf. TP précédent...), le polynôme d'interpolation de $f$
aux points $(x_j, ξ_j, x_{j+1})$ est
donné par
\[
P_j(x) = f[x_j]  + f[x_j, ξ_j] (x-x_j) + f[x_j, ξ_j, x_{j+1}] (x-x_j)(x-ξ_j)
\]
On en déduit, après calculs,
\[
∫_{x_j}^{x_{j+1}} P_j(x) \dd x
= \frac{h}{6}\left(f(x_j) + 4 f(ξ_j) + f(x_{j+1})\right)
\]
Par sommation on obtient
\begin{align*}
I
&≈ \frac{b-a}{6n} \left(
f(x_0) + 4 f\left(\tfrac{x_0+x_1}{2}\right) + 2f(x_1)
+  4 f\left(\tfrac{x_1+x_2}{2}\right) +
⋯
+ 2f(x_{n-1}) +  4 f\left(\tfrac{x_{n-1}+x_n}{2}\right) + f(x_n)
\right)  \\
&≈ \frac{b-a}{6n} \left(f(x_0) + 2∑_{j=1}^{n-1} f(x_j)
+ 4∑_{j=0}^{n-1} f\left(\tfrac{x_{j}+x_{j+1}}{2}\right)
+ f(x_n) \right) \\
&≈ \frac{b-a}{6n} \sum_{k=0}^{2n} w_k f(a_k)
\quad \text{ avec }\quad
a_k = a + k⋅\frac{b-a}{2n}
\text{ et }
w_k = \begin{cases}
1 \text{ si }k=0 \text{ ou }2n\\
4 \text{ si } 0<k<2n \text{ et $k$ impair}\\
2 \text{ si } 0<k<2n \text{ et $k$ pair}\\
\end{cases}
\end{align*}
On peut montrer que l'erreur commise
est majorée par
\[
|E| \leqslant  \frac{1}{2880}⋅h^4⋅\big\|{f^{(4)}}\big|_\infty⋅(b-a)
\tag{∗∗}
\]
On dit alors que la méthode est d'ordre 4; noter que la méthode est
exacte pour tout polynôme de degré au plus 3.



\begin{recherche}
  Programmez   les   trois  méthodes   mentionnées:   \verb+Rectangles+,
  \verb+Trapezes+ et \verb+Simpson+.


On introduira la fonction $x\mapsto \frac{4}{1+x^2}$ et on comparera son
intégrale sur [0,1] avec $\pi$ en utilisant les différentes méthodes.


\end{recherche}

En  pratique,  il  est  difficile  d'avoir un  majorant  de  la  dérivée
quatrième d'une fonction. On préfère donc en pratique  appliquer la méthode avec deux pas
différents ($h$ et $2h$ pour minimiser les évaluations de $f$) et on utilise la
différence des deux approximations numériques comme estimation de l'erreur du
moins bon résultat.




En prenant  5 points  d'interpolation, on construit  la méthode  dite de
\textsc{Boole-Villarceau} et  avec 7 celle de  \textsc{Weddle-Hardy}. On
ne  va  pas plus  loin  car  sinon on  obtient  des  poids négatifs  qui
amplifient les erreurs d'arrondi.




%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\section{La méthode de \textsc{Romberg}}
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

La méthode de \textsc{Romberg} utilise le procédé d'extrapolation de
\textsc{Richardson} pour
accélérer la convergence de la méthode des trapèzes composée.

Soit $f$ une fonction de classe $\mathcal{C}^{∞}$ sur le segment
$[a,b]$. On décompose le segment $[a,b]$ en $n$ sous-intervalles égaux.
Alors l'approximation de l'intégrale $I$ de $f$ sur ce segment par la méthode des
trapèzes vaut
\[
T(h) = h\left(\frac12 f(x_0) + ∑_{j=1}^{n-1} f(x_j) + \frac12 f(x_n)\right)
\quad \text{où}\quad
∀j∈⟦0, n⟧,\quad x_j = a + j h,\quad
\text{et}\enskip h = \frac{b-a}{n}
\]

On peut alors montrer à l'aide de la formule
d'\textsc{Euler-Maclaurin} que
l'erreur commise dans la méthode des trapèzes possède un développement limité à
tout ordre de la forme
\[
E(h) = a_1h^2 + a_2h^4 + ⋯+ a_{k-1}h^{2k-2} + O({h^{2k}})
\quad \text{où}\quad
E(h) = T(h) - ∫_{a}^{b} f(t) \dd t
\]
De plus, les $a_i$ ne dépendent ni de $n$, ni de $h$.

L'idée est alors de faire une combinaison linéaire de $T(h)$ et $T(h/2)$ pour
annuler le premier terme du développement de $E$ :
\[
\renewcommand*\arraystretch{1.3}
\renewcommand*\arraycolsep{1pt}
\begin{array}{rcll}
    T(h) & = &  ∫_{a}^{b} f(t) \dd t
    + a_1h^2 &+ a_2h^4 + ⋯+ a_{k-1}h^{2k-2} + O({h^{2k}}) \\
    T\left(\frac h2\right) & = &  ∫_{a}^{b} f(t) \dd t
    + a_1\frac{h^2}{2^2} &+ a_2\frac{h^4}{2^4} + ⋯+
    a_{k-1}\frac{h^{2k-2}}{2^{2k-2}} + O({h^{2k}}) \\[3pt]
    \hline
    \frac13 \left( 4T\left(\frac h2\right) - T(h)\right)
    & = &  ∫_{a}^{b} f(t) \dd t
    &+ \underbrace{b_2h^4 + ⋯+ b_{k-1}h^{2k-2} + O({h^{2k}})}_{O({h^4}) }
\end{array}
\]
On remarque qu'ainsi, on a gagné un ordre dans le développement de l'erreur.
C'est un bon exercice de vérifier que la nouvelle expression obtenue correspond
à la méthode de \textsc{Simpson}.


Le procédé précédent se généralise aisément : on utilise des dichotomies
successives avec des pas de la forme $h=\frac{b-a}{2^m}$ et on remplit un
tableau triangulaire comprenant les nombres :
\[
A_{m, 0} = T\left(\frac{b-a}{2^m}\right)
\quad \text{et}\quad
∀n∈⟦1, m⟧,\enskip
A_{m, n} = \frac{4^n A_{m, n-1} - A_{m-1,n-1}}{4^n-1}
\]
\begin{center}
\begin{tikzpicture}
    \matrix [column sep=4ex, row sep=1.5ex]
    {
    \node (b) {étape};      & \node {$t=0$};    & \node {$t=1$};  & \node {$t=2$};
    & \node{$⋯$}; & \node (n) {$t=m$}; \\
    \node (t) {tableau[0]};
    & \node [inner sep=1pt, ellipse, fill=gray!30] (x0)   {\small $A_{0,0}$}; \\
    \node {tableau[1]};
    & \node [inner sep=1pt, ellipse, fill=gray!10](x1)   {\small $A_{1,0}$};
    & \node (x11) [inner sep=1pt, ellipse, fill=gray!30]  {\small $A_{1, 1}$}; \\
    \node {tableau[2]};
    & \node [inner sep=1pt, ellipse, fill=gray!10](x2)   {\small $A_{2,0}$};
    & \node (x21) [inner sep=1pt, ellipse, fill=gray!10] {\small $A_{2, 1}$};
    & \node [inner sep=1pt, ellipse, fill=gray!30](xk1)   {\small $A_{2, 2}$};\\
    \node {\small $⋮$}; &\node {\small $⋮$}; &\node (x3) {\small $⋮$}; &\node
    (x31) {\small $⋮$}; & \node {\small $⋱$};\\
    \node [outer xsep=4pt] (c) {tableau[n]};
    & \node (xn)  [inner sep=1pt, ellipse, fill=gray!10] {\small $A_{m,0}$};
    & \node (xn1) [inner sep=1pt, ellipse, fill=gray!10] {\small $A_{m, 1}$};
    & \node [inner sep=1pt, ellipse, fill=gray!10] (xn2) {\small $A_{m, 2}$};
    & \node {\small $⋯$}; & \node [inner sep=1pt, ellipse, fill=gray!30] (xnn)
    {\small $A_{m, m}$};\\
    };
    \draw [thin] (c.south east) -- (c.north east |- b) --++(0,5pt);
    \draw [thin] (n.south west -| b) --++(-25pt, 0) -- (n.south east) --++(15pt,0);
    \draw [->,gray,thick] (x0) .. controls  + (-50:0.5cm) and +(170:0.5cm) .. (x11);
    \draw [->,gray,thick] (x1) .. controls  + (40:0.6cm)  and +(170:0.5cm) .. (x11);
    \draw [->,gray,thick] (x1) .. controls  + (-50:0.5cm) and +(170:0.5cm) .. (x21);
    \draw [->,gray,thick] (x2) .. controls  + (40:0.6cm)  and +(170:0.5cm) .. (x21);
    \draw [->,gray,thick] (x11) .. controls + (-50:0.5cm) and +(170:0.5cm) .. (xk1);
    \draw [->,gray,thick] (x21) .. controls + (40:0.6cm)  and +(170:0.5cm) .. (xk1);
\end{tikzpicture}
\end{center}
Le dernier élément de la diagonale principale $A_{m, m}$ est alors la meilleure
approximation de l'intégrale $I$ que l'on puisse obtenir en évaluant $f$ en
$2^m$ points.
%%% \[
%%% A_{m, m} =   ∫_{a}^{b} f(t) \dd t + O{4^{-m(m+1)}}
%%% \]

En pratique, on commence par remplir la première colonne en se servant de la
relation de récurrence :
\[
A_{k, 0} = \frac12 A_{k-1, 0} + A'_{k, 0}
\quad \text{ où }\quad
A'_{k, 0} = \frac{b-a}{2^k}  ∑_{j=1}^{2^{k-1}} f\left(a+(2j+1)\frac{b-a}{2^k} \right)
\]
En effet, en notant $h=\frac{b-a}{2^k}$, on a
\begin{align*}
    A_{k, 0}
    & = h \left( \frac12 f(a)
    + ∑_{j=1}^{2^k-1} f\left(a+jh \right)
    + \frac12 f(b)\right) \\
    & = h  \left( \frac12 f(a)
    + ∑_{p=1}^{2^{k-1}-1} f\left(a+2ph \right)
    + ∑_{p=0}^{2^{k-1}-1} f\left(a+(2p+1)h \right)
    + \frac12 f(b)\right) \\
    & = \frac12 A_{k-1, 0} +
    \underbrace{h   ∑_{j=0}^{2^{k-1}-1} f\left(a+(2j+1)h \right)}_{A'_{k, 0}}
\end{align*}
Ceci permet de calculer une fois et une fois seulement les valeurs de $f$
nécessaires.

Ensuite, on remplit les autres colonnes en utilisant la formule
\[
∀n∈⟦1, m⟧,\enskip
A_{m, n} = \frac{4^n A_{m, n-1} - A_{m-1,n-1}}{4^n-1}
\]


On montre alors par récurrence que:

$$A_{m,k}=\displaystyle\int_{a}^{b}f(x)\;\mathrm{d}x+O \left(\frac{1}{2^{m(2k+2)}}\right)$$


\begin{recherche}
Programmez Python pour obtenir $A_{m,m}$.

On estime l'erreur comme d'habitude à $\left|A_{m,m}-A_{m-1,m-1}\right|$.
\end{recherche}


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "PolyIPT2014_15"
%%% End: 
