document/GMC1035/diffintegra.tex

  \begin{encadre}{Objectif du chapitre}
   Le calcul dérivé et intégrale est présent dans toutes les modélisations des problèmes d'ingénierie. En pratique, nous avons souvent des modèles sous forme de fonction ce qui rend l'intégration et la dérivation difficile. Dans ce chapitre, l'objectif est d'effectuer des dérivées et des intégrales sans connaître les expressions à dériver ou à intégrer. 
   \end{encadre}
       \\[1cm]
     \begin{encadre}{Connaissances antérieures nécessaires}
      \begin{itemize}
       \item les mathématiques : calcul différentiel, calcul intégral.
       \item programmation de base en Basic ou C : les types, les opérations élémentaires et la structuration de programme en mode console
      \end{itemize}

    \end{encadre}
    
    
\section{Mise en situation}
Dans le chapitre précédent, on a recherché une courbe passant par des points. En pratique, ces points sont souvent des points de mesure. Ce dernier chapitre a donc permis de trouver une fonction $f(x)$ qui passe exactement par ces points.
Dans certaines applications, il peut être intéressant de connaître les fonctions dérivées ou les fonctions intégrales qui passent par ces points.
\\Par exemple dans un GPS, l'appareil effectue des relevés de position. La vitesse de déplacement est alors la dérivée de la fonction $f(x)$ et l'accélération est alors la dérivée seconde de $f(x)$ tandis que la distance parcourue est l'intégrale de $f(x)$.
Dans ce chapitre nous allons établir des formules qui permettent de connaître ces grandeurs sans forcement construire la fonction d'interpolation.
\\Même si ces formules sont destinées à être utilisées sans construire la fonction d'interpolation, elles sont déduites de la fonction d'interpolation.
\section{Différentiation numérique}
Au chapitre précédent on a montré que l'on pouvait interpoler des points d'une fonction $f(x)$ et on 
\begin{eqnarray*}
f(x)&=&p_n(x)+E_n(x)\\
E_n(x)&=&\frac{f^{(n+1)}(\zeta(x))}{(n+1)!}(x-x_0)...(x-x_n)
\end{eqnarray*}
$E_n(x)$ est le terme d'erreur d'ordre $n+1$ (il y a $n+1$ facteurs $x-x_i$).
\subsection{Dérivée d'ordre 1}
On peut calculer la dérivé d'ordre 1 $f'(x)$.
\begin{eqnarray*}
f'(x)&=&p'_n(x)+E'_n(x)\\
E'_n(x)&=&\frac{f^{(n+2)}(\zeta(x))\zeta'(x)}{(n+1)!}(x-x_0)...(x-x_n)+\frac{f^{(n+1)}(\zeta(x))}{(n+1)!}\left[(x-x_0)...(x-x_n)\right]'\\
&=&\frac{f^{(n+2)}(\zeta(x))\zeta'(x)}{(n+1)!}(x-x_0)...(x-x_n)+\frac{f^{(n+1)}(\zeta(x))}{(n+1)!}\left(\sum_{k=0}^{n}\prod_{j=0(j\ne k)}^n(x-x_j)\right)
\end{eqnarray*}
$E'_n(x)$ est le terme d'erreur de $f'(x)$.
\\Pour simplifier le terme d'erreur nous allons calculer l'erreur en un point d'interpolation $x_i$.
\begin{eqnarray*}
E'_n(x_i)&=&\frac{f^{(n+1)}(\zeta(x_i))}{(n+1)!}\left(\prod_{j=0(j\ne i)}^n(x_i-x_j)\right)
\end{eqnarray*}
$E'_n(x)$ est le terme d'erreur. Il est d'ordre $n$ aux points $x_i$ (il y a $n$ facteurs $x-x_i$).
\subsubsection{Approximation d'ordre 1 }
Pour obtenir la dérivée première avec une approximation d'ordre 1, on doit poser $n=1$. On doit trouver un polynôme d'ordre $n=1$. Un polynôme de degré 1 passe par 2 points $x_0$ et $x_1$:
\begin{eqnarray*}
p_1(x)&=&f(x_0)+f[x_0,x_1](x-x_0)\\
f(x)&=&p_1(x)+O(2)\\
f'(x)&=&p'_1(x)+O(1)\\
f'(x)&=&f[x_0,x_1]+O(1)\\
\end{eqnarray*}
On fait l'évaluation en $x_0$
\begin{eqnarray*}
f'(x_0)&=&f[x_0,x_1]=\frac{f(x_1)-f(x_0)}{x_1-x_0}\\
\end{eqnarray*}
Cette expression est la différence avant d'ordre 1.
\\On fait l'évaluation en $x_1$
\begin{eqnarray*}
f'(x_1)&=&f[x_0,x_1]=\frac{f(x_1)-f(x_0)}{x_1-x_0}\\
\end{eqnarray*}
Cette expression est la différence arrière d'ordre 1.


\subsubsection{Approximation d'ordre 2 }
Pour obtenir la dérivée première avec une approximation d'ordre 2, on doit poser $n=2$. On doit trouver un polynôme d'ordre $n=2$. Un polynôme de degré 2 passe par 3 points $x_0$, $x_1$ et $x_2$:
\begin{eqnarray*}
p_2(x)&=&f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)\\
f(x)&=&p_2(x)+O(3)\\
f'(x)&=&p'_2(x)+O(2)\\
f'(x)&=&f[x_0,x_1]+f[x_0,x_1,x_2](2x-x_0-x_1)+O(2)\\
\end{eqnarray*}
On fait l'évaluation en $x_0$ et on suppose que $h=x_1-x_0=x_2-x_1$
\begin{eqnarray*}
f'(x_0)&=&f[x_0,x_1]+f[x_0,x_1,x_2](x_0-x_1)\\
f'(x_0)&=&\frac{f(x_1)-f(x_0)}{h}+\frac{\frac{f(x_2)-f(x_1)}{h}-\frac{f(x_1)-f(x_0)}{h}}{2h}(-h)\\
f'(x_0)&=&\frac{-3f(x_0)+4f(x_1)-f(x_2)}{2h}\\
\end{eqnarray*}
Cette expression est la différence avant d'ordre 2.
\\
On fait l'évaluation en $x_1$ et on suppose que $h=x_1-x_0=x_2-x_1$
\begin{eqnarray*}
f'(x_1)&=&f[x_0,x_1]+f[x_0,x_1,x_2](x_1-x_0)\\
f'(x_1)&=&\frac{f(x_1)-f(x_0)}{h}+\frac{\frac{f(x_2)-f(x_1)}{h}-\frac{f(x_1)-f(x_0)}{h}}{2h}h\\
f'(x_1)&=&\frac{-f(x_0)+f(x_2)}{2h}\\
\end{eqnarray*}
Cette expression est la différence centrée d'ordre 2.
\\
On fait l'évaluation en $x_2$ et on suppose que $h=x_1-x_0=x_2-x_1$
\begin{eqnarray*}
f'(x_2)&=&f[x_0,x_1]+f[x_0,x_1,x_2](2x_2-x_0-x_1)\\
f'(x_2)&=&\frac{f(x_1)-f(x_0)}{h}+\frac{\frac{f(x_2)-f(x_1)}{h}-\frac{f(x_1)-f(x_0)}{h}}{2h}(3h)\\
f'(x_2)&=&\frac{f(x_0)-4f(x_1)+3f(x_2)}{2h}\\
\end{eqnarray*}
Cette expression est la différence arrière d'ordre 2.
\subsubsection{Formule de différences finies pour $f'(x)$ }
Les formules de différences finies qui permettent d'évaluer la dérivée d'une fonction en des points de mesures sont déduites des calculs précédents :
\\
\begin{tabular}{|l|}
  \hline
   \textbf{Formules de différences finies d'ordre 1 pour $f'(x)$}\\
  \hline
  \\
  $f'(x)=\frac{f(x+h)-f(x)}{h}$\\
  \\
  Différence avant d'ordre 1\\
  \hline
  \\
  $f'(x)=\frac{f(x)-f(x-h)}{h}$\\
  \\
  Différence arrière d'ordre 1\\
  \hline
  \textbf{Formules de différences finies d'ordre 2 pour $f'(x)$}\\
  \hline
  \\
  $f'(x)=\frac{-f(x+2h)+4f(x+h)-3f(x)}{2h}$\\
  \\
  Différence avant d'ordre 2\\
  \hline
  \\
  $f'(x)=\frac{f(x+h)-f(x-h)}{2h}$\\
  \\
  Différence centrée d'ordre 2\\
  \hline
  $f'(x)=\frac{3f(x)-4f(x-h)+f(x-2h)}{2h}$\\
  \\
  Différence arrière d'ordre 2\\
  \hline
  
\end{tabular}
\subsection{Dérivée d'ordre supérieur}
Les différences finies pour le calcul des dérivées d'ordre supérieur se fait de la même manière que pour les dérivées premières. Cependant l'ordre d'erreur est plus difficile à déterminer parce que le terme $E_n$ est difficile à dériver plusieurs fois.
A la place on établit les formules des différences finies puis nous déterminons l'ordre en développant la formule des différences finies en série de Taylor pour obtenir la table suivante  :\\ 
\\
\begin{tabular}{|l|}
  \hline
   \textbf{Formules de différences finies d'ordre 1 pour $f''(x)$}\\
  \hline
  \\
  $f''(x)=\frac{f(x+2h)-2f(x+h)+f(x)}{h^2}$\\
  \\
  Différence avant d'ordre 1\\
  \hline
  \\
  $f''(x)=\frac{f(x)-2f(x-h)+f(x-2h)}{h^2}$\\
  \\
  Différence arrière d'ordre 1\\
  \hline
   \textbf{Formules de différences finies d'ordre 2 pour $f''(x)$}\\
  \hline
  \\
  $f''(x)=\frac{f(x+h)-2f(x)+f(x-h)}{h^2}$\\
  \\
  Différence centrée d'ordre 2\\
  \hline
   \textbf{Formules de différences finies d'ordre 4 pour $f''(x)$}\\
  \hline
  \\
  $f''(x)=\frac{-f(x+2h)+16f(x+h)-30f(x)+16f(x-h)-f(x-2h)}{12h^2}$\\
  \\
  Différence centrée d'ordre 4\\
   \hline
   \hline
   \textbf{Formules de différences finies d'ordre 2 pour $f''''(x)$}\\
  \hline
  \\
  $f''''(x)=\frac{f(x+2h)-4f(x+h)+6f(x)-4ff(x-h)+f(x-2h)}{h^4}$\\
  \\
  Différence centrée d'ordre 2\\
   \hline
\end{tabular}


\subsection{Erreur}
L'ensemble des formules de différences finies fait intervenir l'écart entre les points d'interpolation $h$. On pourrait croire que plus $h$ diminue plus le calcul est précis.
En fait, ce n'est pas le cas parce que tous les dénominateurs sont en $h$ donc plus $h$ diminue plus la précision de calcul est affectée. De plus l'ensemble des différences finies s'approche de l'indéterminé $\frac{0}{0}$. Le processus de dérivation est un processus instable.
\subsection{Exemple}
Soit $f(x)=e^x$ que l'on veut dériver en $x=0$. Dans ce cas on sait que $f'(0)=e^0=1$. \\On veut vérifier la dérivation numérique et la comparer à la valeur théorique.
\\
\\À l'ordre 1  on a \\
\\
\begin{tabular}{|l|l|l|l|l|}
  \hline
   h & $f'(0)=\frac{e^h-1}{h}$ & erreur $E$ & $\frac{h_i}{h_{i-1}}$ & $\frac{E_i}{E_{i-1}}$ \\
  \hline
  $0.1$ & $1.051709$ & $0.051709$ & & \\
  $0.05$ & $1.0254219$ & $0.0254219$ & $0.5$ & $0.491633951$\\
  $0.025$ & $1.012604821$ & $0.012604821$ & $0.5$ & $0.495825292$\\
  \hline
\end{tabular}
\\ \\
Ceci confirme que l'on a une convergence d'ordre 1.\\
\\À l'ordre 2  on a \\
\\
\begin{tabular}{|l|l|l|l|l|}
  \hline
   h & $f'(0)=\frac{e^h-e^{-h}}{2h}$ & erreur $E$ & $\frac{h_i}{h_{i-1}}$ & $\frac{E_i}{E_{i-1}}$ \\
  \hline
  $0.1$ & $1.0016675$ & $0.0016675$ & & \\
  $0.05$ & $1.000416719$ & $0.000416719$ & $0.5$ & $0.249906446=0.499906438^2$\\
  $0.025$ & $1.00010417$ & $0.00010417$ & $0.5$ & $0.249976406=0.499976405^2$\\
  \hline
\end{tabular}
\\ \\
Ceci confirme que l'on a une convergence d'ordre 2.\\
  
\section{Extrapolation de Richardson}
\subsection{Construction de l'extrapolation de Richardson}
L'extrapolation de Richardson est une technique qui permet d'augmenter la précision d'une approximation. Soit une approximation $Q_{app}(h)$ d'ordre $n$ résultant d'une quelconque méthode d'approximation. La solution exacte $Q_{exa}$ peut s'écrire : 
\begin{eqnarray*}
Q_{exa}&=&Q_{app}(h)+O(h^n)\\
&=&Q_{app}(h)+C_nh^n+C_{n+1}h^{n+1}+C_{n+2}h^{n+2}+...
\end{eqnarray*}
où les différents coefficients $C$ sont dépendants de la méthode d'extrapolation.
\\Calculons maintenant la valeur de l'approximation pour un pas de $\frac{h}{2}$ : 
\begin{eqnarray*}
Q_{exa}&=&Q_{app}\left(\frac{h}{2}\right) +C_n \left( \frac{h}{2} \right)^n+C_{n+1}\left( \frac{h}{2} \right)^{n+1}+C_{n+2}\left( \frac{h}{2} \right)^{n+2}+...\\
2^nQ_{exa}&=&2^nQ_{app}\left(\frac{h}{2}\right) +C_nh^n +C_{n+1}\frac{h^{n+1}}{2} +C_{n+2}\frac{h^{n+2}}{2^2} +...\\
\end{eqnarray*}
En soustrayant les deux dernières équations il vient 
\begin{eqnarray*}
(2^n-1)Q_{exa}&=&2^nQ_{app}\left(\frac{h}{2}\right)-Q_{app}(h) -\frac{1}{2}C_{n+1}\frac{h^{n+1}}{2} -\frac{3}{4}C_{n+2}\frac{h^{n+2}}{2^2} +...\\
Q_{exa}&=&\frac{2^nQ_{app}\left(\frac{h}{2}\right)-Q_{app}(h)}{2^n-1}+\frac{ -\frac{1}{2}C_{n+1}\frac{h^{n+1}}{2} -\frac{3}{4}C_{n+2}\frac{h^{n+2}}{2^2} +...}{2^n-1}\\
&=&\frac{2^nQ_{app}\left(\frac{h}{2}\right)-Q_{app}(h)}{2^n-1}+O(h^{n+1})\\
\end{eqnarray*}
Le terme de droite est une approximation de la solution exacte $Q_{exa}$ d'ordre $n+1$. L'ordre de convergence a donc augmenté au moins de un. Selon les valeur de $C_{n+1}$ l'augmentation d'ordre pourrait être plus importante.
\subsection{Exemple}
On reprend l'exercice précédent et on effectue une extrapolation de Richardson.
\\
\\À l'ordre 1  on a \\
\\
\begin{tabular}{|l|l|l|l|l|l|}
  \hline
   h & $f'(0)=\frac{e^h-1}{h}$ &$\frac{2f'(h_i)-f'(h_{i-1})}{2-1}$  & erreur $E$ & $\frac{h_i}{h_{i-1}}$ & $\frac{E_i}{E_{i-1}}$ \\
  \hline
  $0.1$ & $1.051709$ & & & & \\
  $0.05$ & $1.0254219$ & $0.9991348$ & $-0.0008652$ & & \\
  $0.025$ & $1.012604821$ & $0.9999787742$& $-0.000212258$ & $0.5$ & $0.245328247=0.495306216^2$\\
  \hline
\end{tabular}
\\ \\
Ceci confirme que l'on a une convergence d'ordre 2.\\
\\À l'ordre 2  on a \\
\\
\begin{tabular}{|l|l|l|l|l|l|}
  \hline
   h & $f'(0)=\frac{e^h-e^{-h}}{2h}$&$\frac{2^2f'(h_i)-f'(h_{i-1})}{2^2-1}$  & erreur $E$ & $\frac{h_i}{h_{i-1}}$ & $\frac{E_i}{E_{i-1}}$ \\
  \hline
  $0.1$ & $1.0016675$ & & & & \\
  $0.05$ & $1.000416719$ & $0.999999792$ &$-0.000000208 $ &  & \\
  $0.025$ & $1.00010417$ & $0.999999987$ & $-0.000000013$& $0.5$ & $0.0625=0.5^4$\\
  \hline
\end{tabular}
\\ \\
Dans ce cas on a une convergence d'ordre 4.\\Ceci est explicable puisque l'erreur des différences finies centrées d'ordre 2 comporte des termes d'erreur en puissance paire uniquement. Dans l'erreur de l'extrapolation de Richardson le terme $C_{n+1}$ est nul, l'ordre de convergence est donc augmenté de 2.
\\
Dans le paragraphe précédent on a vu que la dérivation numérique est instable et que la diminution de $h$ ne permet pas d'augmenter forcement la précision de la méthode. En revanche une extrapolation de Richardson permet d'augmenter la précision sans diminuer $h$.
\section{Intégration numérique}
Au chapitre précédent on a montré que l'on pouvait interpoler des points d'une fonction $f(x)$ et on 
\begin{eqnarray*}
f(x)&=&p_n(x)+E_n(x)\\
E_n(x)&=&\frac{f^{(n+1)}(\zeta(x))}{(n+1)!}(x-x_0)...(x-x_n)
\end{eqnarray*}
$E_n(x)$ est le terme d'erreur d'ordre $n+1$ (il y a $n+1$ facteurs $x-x_i$).
Le problème de l'intégration revient à calculer
\begin{eqnarray*}
\int_{x_0}^{x_n} f(x)dx&=&\int_{x_0}^{x_n} p_n(x)dx+\int_{x_0}^{x_n} E_n(x)dx\\
\end{eqnarray*}

\subsection{Formule de Newton-Cotes}
\subsubsection{Méthode des trapèzes}
La fonction est connue en deux points donc $n=1$ et on a 
\begin{eqnarray*}
\int_{x_0}^{x_1} f(x)dx&=&\int_{x_0}^{x_1} p_1(x)dx+\int_{x_0}^{x_1} E_1(x)dx\\
E_1(x)&=&\frac{f''(\zeta(x))}{(2)!}(x-x_0)(x-x_1)\\
\int_{x_0}^{x_1} f(x)dx&=&\int_{x_0}^{x_1}\left( f(x_0)+f[x_0,x_1](x-x_0)  \right) dx+\int_{x_0}^{x_1}\frac{f''(\zeta(x))}{(2)!}(x-x_0)(x-x_1) dx\\
&=&\left[ f(x_0)x \right]_{x_0}^{x_1}+\left[ \frac{f(x_1)-f(x_0)}{x_1-x_0}\frac{(x-x_0)^2}{2} \right]_{x_0}^{x_1}+\int_{x_0}^{x_1}\frac{f''(\zeta(x))}{(2)!}(x-x_0)(x-x_1) dx\\
&=& f(x_0)(x_1-x_0) +\frac{(f(x_1)-f(x_0)(x_1-x_0)}{2} +\int_{x_0}^{x_1}\frac{f''(\zeta(x))}{(2)!}(x-x_0)(x-x_1) dx\\
&=& (x_1-x_0)\frac{f(x_0)+f(x_1)}{2} +\int_{x_0}^{x_1}\frac{f''(\zeta(x))}{(2)!}(x-x_0)(x-x_1) dx\\
\end{eqnarray*}
La partie du résultat de droite correspond à l'aire du trapèze qui est se situe sous le polynôme entre les points $x_0$ et $x_1$.  Il reste le terme d'erreur à calculer.
\begin{eqnarray*}
E_1(x)&=&\int_{x_0}^{x_1} \frac{f''(\zeta(x))}{(2)!}(x-x_0)(x-x_1)dx\\
\end{eqnarray*}
On effectue le changement de variable suivant 
\begin{eqnarray*}
h&=&x_1-x_0\\
s&=&\frac{x-x_0}{h}\quad x=sh+x_0\\
dx&=&hds\\
\end{eqnarray*}
et on obtient 
\begin{eqnarray*}
E_1(x)&=&\int_{0}^{1} \frac{f''(\zeta(x))}{(2)!}sh^2(s-1)hds\\
&=&h^3\int_{0}^{1} \frac{f''(\zeta(x))}{2}s(s-1)ds=O(h^3)\\
\end{eqnarray*}
La méthode des trapèzes a un ordre de convergence de 3.
\\Exemple :
On sait que $\int_0^\frac{\pi}{2}\sin x dx=1$.
Le calcul de l'intégrale numérique donne
\begin{eqnarray*}
\int_0^\frac{\pi}{2}\sin x dx&=&(\frac{\pi}{2}-0)\frac{\sin \frac{\pi}{2}-\sin 0}{2}=\frac{\pi}{4}
\end{eqnarray*}
Ce calcul donne un très mauvais résultat. Le problème ne vient pas de la méthode d'intégration elle-même mais du fait que l'on a interpolé une fonction sinus avec un polynôme de degré 1. Pour améliorer le tout on peut découper l'intervalle en $n$ intervalle de largeur $h=\frac{b-a}{n}$ :
\begin{eqnarray*}
\int_a^b f(x)dx&=&\sum_{i=0}^{n-1}\int_{x_i}^{x_{i+1}} f(x)dx\\
&=&\sum_{i=0}^{n-1}\frac{h}{2}\left( f(x_i)+f(x_{i+1})\right)\\
&=&\frac{h}{2}\left( f(x_0)+f(x_1)+f(x_1)+f(x_2)+...+f(x_n)\right)\\
&=&\frac{h}{2}\left( f(x_0)+2f(x_1)+2f(x_2)+...+2f(x_{n-1})+f(x_n)\right)\\
\end{eqnarray*}
Le terme d'erreur devient 
\begin{eqnarray*}
E_1(x)&=&\int_{a}^{b} \frac{f''(\zeta(x))}{(2)!}(x-x_0)(x-x_1)dx\\
&=&\sum_{i=0}^{n-1}\int_{x_i}^{x_{i+1}} \frac{f''(\zeta(x))}{(2)!}(x-x_0)(x-x_1)dx\\
&=&\sum_{i=0}^{n-1}h^3\int_{0}^{1} \frac{f''(\zeta(x))}{2}s(s-1)ds\\
&=&nh^3\int_{0}^{1} \frac{f''(\zeta(x))}{2}s(s-1)ds\\
&=&\frac{b-a}{h}h^3\int_{0}^{1} \frac{f''(\zeta(x))}{2}s(s-1)ds\\
&=&(b-a)h^2\int_{0}^{1} \frac{f''(\zeta(x))}{2}s(s-1)ds\\
&=&O(h^2)
\end{eqnarray*}
La méthode des trapèzes avec des intervalles a une convergence d'ordre 2.
\\On revient au calcul de $\int_0^\frac{\pi}{2}\sin x dx$
\begin{itemize}
 \item avec 4 intervalles
\begin{eqnarray*}
h&=&\frac{\frac{\pi}{2}-0}{4}=\frac{\pi}{8}\\
\int_0^\frac{\pi}{2}\sin x dx&=&\frac{\pi}{16}\left(\sin 0+2\sin\frac{\pi}{8}+2\sin\frac{2\pi}{8}+2\sin\frac{3\pi}{8}+\sin\frac{4\pi}{8} \right)\\&=&0.9871158\\
\end{eqnarray*}

 \item avec 8 intervalles
\begin{eqnarray*}
h&=&\frac{\frac{\pi}{2}-0}{8}=\frac{\pi}{16}\\
\int_0^\frac{\pi}{2}\sin x dx&=&\frac{\pi}{32}(\sin 0+2\sin\frac{\pi}{16}+2\sin\frac{2\pi}{16}+2\sin\frac{3\pi}{16}+2\sin\frac{4\pi}{16}+\\& & 2\sin\frac{4\pi}{16}+2\sin\frac{5\pi}{16}+2\sin\frac{7\pi}{16}+\sin\frac{8\pi}{16} )
\\&=&0.9967852\\
\end{eqnarray*}
\item L'erreur à 4 intervalles est $-0.0128842$ et l'erreur à 8 intervalles est $-0.0032148$ soit un rapport de $0.249514909=0.499514674^2$. Ce qui correspond à une convergence d'ordre 2.
\end{itemize}

\subsubsection{Méthode de Simpson $\frac{1}{3}$}
On refait le même raisonnement que la méthode des trapèzes mais on applique sur le polynôme d'interpolation de degré 2 $p_2(x)$ qui nécessite donc trois points d'interpolation et on obtient
\begin{eqnarray*}
\int_{x_0}^{x_2} f(x)dx&=&\frac{h}{3}\left( f(x_0)+4f(x_1)+f(x_2) \right)\\
\end{eqnarray*}
avec une convergence d'ordre 5.
\\La version avec division en $2n$ intervalles donne 
\begin{eqnarray*}
\int_{a}^{b} f(x)dx&=&\frac{h}{3}\left( f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+2f(x_4)+...+2f(x_{2n-2})+4f(x_{2n-1})+f(x_{2n}) \right)\\
\end{eqnarray*}
avec une convergence d'ordre 4.
\\Exemple : On revient au calcul de $\int_0^\frac{\pi}{2}\sin x dx$
\begin{eqnarray*}
\int_0^\frac{\pi}{2}\sin x dx=\frac{\frac{\frac{\pi}{2}-0}{2}}{3}(\sin 0 +4\sin \frac{\pi}{4}+\sin \frac{\pi}{2})=1.002279877
\end{eqnarray*}

\subsubsection{Méthode de Simpson $\frac{3}{8}$}
On refait le même raisonnement que la méthode des trapèzes mais on applique sur le polynôme d'interpolation de degré 3 $p_3(x)$ qui nécessite donc quatre points d'interpolation et on obtient
\begin{eqnarray*}
\int_{x_0}^{x_3} f(x)dx&=&\frac{3h}{8}\left( f(x_0)+3f(x_1)+3f(x_2)+f(x_3) \right)\\
\end{eqnarray*}
avec une convergence d'ordre 5.
\\La version avec division en $3n$ intervalles donne 
\begin{eqnarray*}
\int_{a}^{b} f(x)dx&=&\frac{3h}{8}\left( f(x_0)+3f(x_1)+3f(x_2)+2f(x_3)+3f(x_4)+...+2f(x_{3n-3})+3f(x_{3n-2})+3f(x_{3n-1})+f(x_{3n}) \right)\\
\end{eqnarray*}
avec une convergence d'ordre 4.
\\Exemple : On revient au calcul de $\int_0^\frac{\pi}{2}\sin x dx$
\begin{eqnarray*}
\int_0^\frac{\pi}{2}\sin x dx=\frac{3\frac{\frac{\pi}{2}-0}{3}}{8}(\sin 0 +3\sin \frac{\pi}{6}+3\sin \frac{2\pi}{6}+\sin \frac{3\pi}{6})=1.001004923
\end{eqnarray*}
\subsection{Formule de Gauss-Legendre}
Dans les formules de Newton-Cote, les points d'intégration sont fixés aux points d'interpolation. Dans la méthode de Gauss-Legendre on va chercher si d'autres points d'intégration permettent d'obtenir des résultats plus précis.
On voudrait trouver les $x_i$ pour que 
\begin{eqnarray*}
\int_a^bf(x)dx=\sum_{i=1}^r w_if(x_i)
\end{eqnarray*}
où $r$ est le nombre de points de Gauss. On appelle cette expression la quadrature de Gauss-Legendre à $r$ points. Les $x_i$ sont appelés les points d'intégrations et les $w_i$ sont les poids d'intégration.
Pour simplifier on s'intéresse tout d'abord à 
\begin{eqnarray*}
\int_{-1}^1f(x)dx=\sum_{i=1}^r w_if(x_i)
\end{eqnarray*}
\begin{itemize}
 \item si $r=1$, on a
\begin{eqnarray*}
\int_{-1}^1f(x)dx=w_1f(x_1)
\end{eqnarray*}
Soit deux inconnues $w_1$ et $x_1$. Il faut donc deux équations pour résoudre le problème. On prends alors un polynôme de degré 0 et un degré 1 pour établir les équations.
\begin{eqnarray*}
\int_{-1}^11dx&=&w_1*1=[x]_{-1}^1=2\\
\int_{-1}^1xdx&=&w_1*x_1=\left[\frac{x^2}{2}\right]_{-1}^1=0\\
\end{eqnarray*}
Il vient que $w_1=1$ et $x_1=0$ et qu'un polynôme de degré 1 est intégré de manière exact.
\item si $r=2$, on a
\begin{eqnarray*}
\int_{-1}^1f(x)dx=w_1f(x_1)+w_2f(x_2)
\end{eqnarray*}
Soit quatre inconnues $w_1$,$x_1$,$w_2$et$x_2$. Il faut donc quatre équations pour résoudre le problème. On prends alors des polynômes de degré 0,1,2 et 3 pour établir les équations.
\begin{eqnarray*}
\int_{-1}^11dx&=&w_1*1+w_2*1=[x]_{-1}^1=2\\
\int_{-1}^1xdx&=&w_1*x_1+w_2*x_2=\left[\frac{x^2}{2}\right]_{-1}^1=0\\
\int_{-1}^1x^2dx&=&w_1*x_1^2+w_2*x_2^2=\left[\frac{x^3}{3}\right]_{-1}^1=\frac{2}{3}\\
\int_{-1}^1x^3dx&=&w_1*x_1^3+w_2*x_2^3=\left[\frac{x^4}{4}\right]_{-1}^1=0\\
\end{eqnarray*}
on obtient un système d'équations non linéaires
\begin{eqnarray*}
\left\{
\begin{array}{l}
w_1+w_2=2\\
w_1x_1+w_2x_2=0\\
w_1x_1^2+w_2x_2^2=\frac{2}{3}\\
w_1x_1^3+w_2x_2^3=0\\
\end{array} \right.
\end{eqnarray*}
On multiplie la deuxième équation par $x_1^2$ et on la soustrait avec la quatrième, on obtient
\begin{eqnarray*}
w_2x_2x_1^2-w_2x_2^3=0\\ 
w_2x_2(x_1^2-x_2^2)=0\\ 
\end{eqnarray*}
Cette dernière équation s'annule pour $x_1=-x_2$. en effet tout les autres cas d'annulation ($w_2=0$,$x_2=0$ et $x_1=x_2$) conduisent à une formule de Gauss à un point. Le deuxième équation devient alors $x_1(w_1-w_2)$ ce qui conduit à $w_1=w_2=1$ (d'après la première équation).
\\
La troisième équation amène alors à $2x_1^2=\frac{2}{3}$ soit $x_1=\pm\sqrt{\frac{1}{3}}$.
\\On en déduit que $w_1=w_2=1$ et que $x_1=-\sqrt{\frac{1}{3}}$, $x_2=\sqrt{\frac{1}{3}}$
On a donc $\int_{-1}^1f(x)dx=f\left(-\sqrt{\frac{1}{3}}\right)+f\left(\sqrt{\frac{1}{3}}\right)$ et cette formule intègre exactement un polynôme de degré 3.
\item en continuant avec des valeurs de r plus élevé on trouve la table d'intégration de Gauss Legendre qui intègre exactement un polynôme de degré $2r-1$
\\
\begin{tabular}{|c|c|c|c|}
  \hline
   \multicolumn{4}{|c|}{\textbf{Quadrature de Gauss-Legendre}}\\
  \hline
 $r$ & $x_i$ & $w_i$ & Degré de précision\\
  \hline
$1$ & $0$ & $2$ & $1$\\
  \hline
$2$ & $-\sqrt{\frac{1}{3}}$ & $1$ & $3$\\
 & $\sqrt{\frac{1}{3}}$ & $1$ &\\
 \hline
$3$ & $0$ & $\frac{8}{9}$ & $5$\\
 & $-\frac{\sqrt{15}}{5}$ & $\frac{5}{9}$ & \\
 & $\frac{\sqrt{15}}{5}$ & $\frac{5}{9}$ & \\
 \hline
$4$ & $-\frac{\sqrt{525+70\sqrt{30}}}{35}$ & $\frac{18-\sqrt{30}}{36}$ & $7$\\
 & $-\frac{\sqrt{525-70\sqrt{30}}}{35}$ & $\frac{18+\sqrt{30}}{36}$ & \\
 & $\frac{\sqrt{525-70\sqrt{30}}}{35}$ & $\frac{18+\sqrt{30}}{36}$ & \\
 & $\frac{\sqrt{525+70\sqrt{30}}}{35}$ & $\frac{18-\sqrt{30}}{36}$ & \\
 \hline
$5$ & $0$ & $\frac{128}{225}$ & $9$\\
 & $-\frac{\sqrt{245+14\sqrt{70}}}{35}$ & $\frac{322-13\sqrt{70}}{900}$ & \\
 & $-\frac{\sqrt{245-14\sqrt{70}}}{35}$ & $\frac{322+13\sqrt{70}}{900}$ &\\
 & $\frac{\sqrt{245-14\sqrt{70}}}{35}$ & $\frac{322+13\sqrt{70}}{900}$ &\\
 & $\frac{\sqrt{245+14\sqrt{70}}}{35}$ & $\frac{322-13\sqrt{70}}{900}$ &\\
 \hline
  \end{tabular}
\end{itemize}
Maintenant nous pouvons revenir au problème initial qui est le calcul de  $\int_a^b f(x)dx$.
On effectue un changement de variable 
\begin{eqnarray*}
x&=&\frac{(b-a)t+(a+b)}{2}\\
dx&=&\frac{b-a}{2}dt
\end{eqnarray*}
et on obtient 
\begin{eqnarray*}
\int_a^b f(x)dx&=&\frac{b-a}{2}\int_{-1}^1 f\left( \frac{(b-a)t+(a+b)}{2} \right) dt
\end{eqnarray*}
et on peut appliquer la quadrature de Gauss à la fonction $f$ avec le changement de variable de $x$ en $t$.
\\Pour augmenter la précision on peut aussi diviser le calcul en plusieurs intervalles
\begin{eqnarray*}
\int_a^b f(x)dx&=&\sum_{i=0}^{n-1}   \int_{x_i}^{x_{i+1}} f(x) dx\\
&=&\sum_{i=0}^{n-1} \frac{x_{i+1}-x_i}{2}  \int_{-1}^{1} f\left(\frac{(x_{i+1}-x_i)t+(x_{i+1}+x_i)}{2}\right) dt\\
&=&\sum_{i=0}^{n-1} \frac{x_{i+1}-x_i}{2}  \sum_{j=1}^{r} w_jf\left(\frac{(x_{i+1}-x_i)t_j+(x_{i+1}+x_i)}{2}\right)
\end{eqnarray*}
\\
Exemple : On revient sur l'exemple de $\int_0^{\frac{\pi}{2}}\sin x dx$ et on utilise la quadrature à 2 points
\begin{eqnarray*}
\int_0^{\frac{\pi}{2}}\sin x dx&=&\frac{\frac{\pi}{2}-0}{2}\int_{-1}^1 \sin \frac{\pi t+\pi}{4} dt\\
&=&\frac{\pi}{4}\int_{-1}^1 \sin \frac{\pi t+\pi}{4} dt\\
&=&\frac{\pi}{4} \left( \sin \frac{-\frac{\pi}{\sqrt{3}}+\pi}{4}+\sin \frac{\frac{\pi}{\sqrt{3}}+\pi}{4}\right)\\
&=&0.998472
\end{eqnarray*}
On fait le même calcul avec deux intervalles.
\begin{eqnarray*}
\int_0^{\frac{\pi}{2}}\sin x dx&=&\sum_{i=0}^{1} \frac{x_{i+1}-x_i}{2}  \sum_{j=1}^{2} w_j\sin \frac{(x_{i+1}-x_i)t_j+(x_{i+1}+x_i)}{2}\\
&=&\frac{\frac{\pi}{4}-0}{2}  \sum_{j=1}^{2} w_j\sin \frac{(\frac{\pi}{4}-0)t_j+(\frac{\pi}{4}+0)}{2}+\frac{\frac{\pi}{2}-\frac{\pi}{4}}{2}  \sum_{j=1}^{2} w_j\sin \frac{(\frac{\pi}{2}-\frac{\pi}{4})t_j+(\frac{\pi}{2}+\frac{\pi}{4})}{2}\\
&=&\frac{\pi}{8}  \sum_{j=1}^{2} w_j\sin \frac{\pi t_j+\pi}{8}+\frac{\pi}{8} \sum_{j=1}^{2} w_j\sin \frac{\pi t_j+3\pi}{8}\\
&=&\frac{\pi}{8}  \sin \frac{-\pi \sqrt{\frac{1}{3}}+\pi}{8}+\frac{\pi}{8}  \sin \frac{\pi \sqrt{\frac{1}{3}}+\pi}{8}+\frac{\pi}{8} \sin \frac{-\pi \sqrt{\frac{1}{3}}+3\pi}{8}+\frac{\pi}{8} \sin \frac{\pi \sqrt{\frac{1}{3}}+3\pi}{8}\\
&=&0.99991017
\end{eqnarray*}
\begin{encadre}{À retenir}
Il est relativement aisé d'intégrer numériquement de manière précise. Au niveau de la dérivation cela est plus difficile.
\end{encadre}
Revision:	948
Committed:	Wed Aug 8 13:46:37 2018 UTC (6 years, 9 months ago) by francois
Content type:	application/x-tex
File size:	24482 byte(s)
Log Message:	mise a jour note de cours GMC1035