% Jean-Pierre Demailly
% Universit\'e de Grenoble I, Institut Fourier
% \'Equations Diff\'erentielles et Analyse Num\'erique
% Collection Grenoble-Sciences, Presses Universitaires de Grenoble
% Chapitre IX

\input equa-mac.tex

\blankline\blankline 
\hfill\chtitle Chapitre IX|%
\blankline\blankline 
\hfill\title M\'ethodes \`a pas multiples|
\chaptitlerunning IX -- M\'ethodes \`a pas multiples|
\vskip4pc

\noindent Comme dans le chapitre pr\'ec\'edent, on s'int\'eresse \`a la 
r\'esolution num\'erique du probl\`eme de Cauchy relatif \`a une
\'equation diff\'erentielle
$$
y'=f(t,y),\qquad (t,y)\in[t_0,t_0+T]\times\bR.\leqno {\rm (E)}
$$
Si $(t_n)_{0\leq n\leq N}$ est une subdivision de $[t_0,t_0+T]$ de
pas successifs $h_n=t_{n+1}-t_n$, on appelle {\it m\'ethode
num\'erique \`a $r+1$ pas} toute m\'ethode num\'erique de la forme
$$
y_{n+1}=\Psi(t_n,y_n,h_n\,;\dots;\,t_{n-r},y_{n-r},h_{n-r}).
$$
L'int\'er\^et de ces m\'ethodes vient du fait qu'on peut obtenir un
ordre \'elev\'e pour une complexit\'e de calcul nettement inf\'erieure
\`a celle des m\'ethodes de Runge-Kutta. L'un des probl\`emes
essentiels, n\'eanmoins, est de s'assurer que la stabilit\'e
num\'erique reste suffisamment bonne.

\section 1. Une classe de m\'ethodes avec pas constant|

\noindent On suppose ici que le pas $h_n=h$ est constant. On s'int\'eresse
aux m\'ethodes \`a $r+1$ pas permettant un calcul r\'ecurrent des
points $(t_n,y_n)$ et des pentes $f_n=f(t_n,y_n)$ sous la forme
$$
{\rm (M)}\quad\cases{
\displaystyle y_{n+1}=\sum_{0\leq i\leq r}\alpha_iy_{n-i}+h\sum_{0\leq i\leq
r}\beta_if_{n-i}\hfill\cr 
\noalign{\vskip3pt}
t_{n+1}=t_n+h\cr
\noalign{\vskip3pt}
f_{n+1}=f(t_{n+1},y_{n+1})\cr
}
$$
o\`u les $\alpha_i$, $\beta_i$, $0\leq i\leq r$ sont des constantes r\'eelles.

\bigskip
\state D\'emarrage de l'algorithme -- |  Le point initial $(t_0,y_0)$ \'etant
donn\'e, l'algo\-rithme ne peut d\'emarrer que si les valeurs
$(y_1,f_1),\dots,(y_r,f_r)$ ont d\'ej\`a \'et\'e calcul\'ees.

Ce calcul ne peut \^etre fait que par une m\'ethode \`a un pas pour
$(y_1,f_1)$, \`a au plus 2 pas pour $(y_2,f_2)$, $\dots$ au plus $r$
pas pour $(y_r,f_r)$. L'initialisation des $r$ premi\`eres valeurs
$(y_i,f_i)$, $1\leq i\leq r$, sera g\'en\'eralement faite \`a l'aide
d'une m\'ethode de Runge-Kutta d'ordre sup\'erieur ou \'egal \`a celui
de la m\'ethode (M), ou \`a la rigueur un de moins (voir le d\'ebut du
\S\ 1.2 sur ce point).

\subsection 1.1. Erreur de consistance et ordre|

\noindent La d\'efinition g\'en\'erale de l'erreur de consistance pour une
m\'ethode \`a $r+1$ pas est la suivante (on ne suppose pas
n\'ecessairement dans cette d\'efinition que le pas est constant).

\bigskip
\majdeb
\state D\'efinition -- |  {\it Soit $z$ une solution exacte de l'\'equation 
{\rm (E)}. {\it L'erreur de consistance} $e_n$ relative \`a $z$ est l'\'ecart
$$
e_n=z(t_{n+1})-y_{n+1},\quad r\leq n<N,
$$
obtenu en calculant $y_{n+1}$ \`a partir  des $r+1$ valeurs pr\'ec\'edentes
suppos\'ees exactes
\hfill\break $y_n=z(t_n),\dots,y_{n-r}=z(t_{n-r})$.\hfill\break
La m\'ethode est dite {\rm d'ordre $p$} si pour toute solution $z$ il existe
une constante $C$ telle que 
$$
|e_n|\leq Ch_nh^p_{\rm max}.
$$}
\vskip-\baselineskip
\majfin
\vskip\baselineskip

D\'eterminons $e_n$ dans le cas de la m\'ethode (M) ci-dessus. On a
$$
z(t_{n+1})=z(t_n+h)=\sum_{0\leq k\leq p}\frac{h^k}{k!}\ z^{(k)}(t_n)+O(h^{p+1})
$$
d\`es que $f$ est de classe $C^p$ ($z$ est alors de classe $C^{p+1}$). Par ailleurs
$$
\eqalign{
y_{n-i}&=z(t_{n-i})=z(t_n-ih)=\sum_{0\leq k\leq p}\frac{(-ih)^k}{k!}\ z^{(k)}(t_n)+O(h^{p+1}),\cr
f_{n-i}&=f(t_{n-i},z(t_{n-i}))=z'(t_{n-i})=z'(t_n-ih)\cr
&=\sum_{0\leq k\leq p-1}\frac{(-ih)^k}{k!}\ z^{(k+1)}(t_n)+O(h^p)\cr
&=\sum_{0\leq k\leq p}k\,\frac{(-ih)^{k-1}}{k!}\ z^{(k)}(t_n)+O(h^p).\cr
}
$$
Il vient par cons\'equent
$$
\eqalign{
e_n&=z(t_{n+1})-y_{n+1}=z(t_{n+1})-\sum_{0\leq i\leq r}(\alpha_iy_{n-i}+h\beta_i f_{n-i})\cr
&=\sum_{0\leq k\leq p}\frac{h^k}{k!}\ z^{(k)}(t_n)\Big[1-\sum_{0\leq i\leq
r}(\alpha_i(-i)^k+k\beta_i(-i)^{k-1})\Big]+O(h^{p+1})\cr
&=\sum_{0\leq k\leq p}\frac{h^k}{k!}\ z^{(k)}(t_n)\Big[1-(-1)^k\sum_{0\leq i\leq
r}i^k\alpha_i-ki^{k-1}\beta_i\Big]+O(h^{p+1}).\cr  
}
$$
La m\'ethode (M) est donc d'ordre $\geq p$ si et seulement si elle v\'erifie les conditions
$$
\sum_{0\leq i\leq r}i^k\alpha_i-ki^{k-1}\beta_i=(-1)^k,\qquad 0\leq k\leq p.
$$
En particulier, elle est d'ordre $\geq 1$ (ou consistante) si et seulement si
$$
\cases{
\alpha_0+\alpha_1+\dots+\alpha_r=1\hfill\cr
\noalign{\vskip3pt}
\qquad\alpha_1+\dots+r\alpha_r-(\beta_0+\dots+\beta_r)=-1.\cr
}
$$
Pour qu'elle soit d'ordre $\geq p$, la $p$-i\`eme condition s'explicite de la mani\`ere suivante~:
$$
\alpha_1+2^p\alpha_2+\dots+r^p\alpha_r-p(\beta_1+2^{p-1}\beta_2+\dots+r^{p-1}\beta_r)=(-1)^p.
$$

\subsection 1.2. Stabilit\'e|

\noindent On dit qu'une m\'ethode \`a pas multiples est stable si une petite
perturbation des valeurs initiales $y_0,\dots,y_r$ et de petites erreurs
$\varepsilon_n$ dans le calcul r\'ecurrent des valeurs $y_{n+1}$, 
$r\leq n<N$, provoquent une erreur finale contr\^olable. De fa\c{c}on 
pr\'ecise~:

\bigskip
\majdeb
\state D\'efinition -- |  {\it On dit qu'une m\'ethode \`a $r+1$ pas est
stable de constante de stabilit\'e $S$ si pour toutes suites $y_n$, 
$\widetilde y_n$ avec
$$
\eqalign{
y_{n+1}&=\Psi (t_{n-i},y_{n-i},h_{n-i}),\qquad\quad\kern1pt r\leq n<N,\cr
\widetilde y_{n+1}&=\Psi(t_{n-i},\widetilde y_{n-i},h_{n-i})+\varepsilon_n,
\quad r\leq n<N,\cr}
$$
alors
$$
\max_{0\leq n\leq N}|\widetilde y_n-y_n|\leq S\Big(\max_{0\leq n\leq r}|\widetilde
y_n-y_n|+\sum_{r\leq n<N}|\varepsilon_n|\Big). 
$$}
\vskip-\baselineskip
\majfin
\vskip\baselineskip

En appliquant cette d\'efinition \`a $\widetilde y_n=z(t_n)$, on voit que l'erreur globale de
la suite $y_n$ par rapport \`a la solution exacte $z(t_n)$ admet la majoration

\bigskip
\majdeb
\hfil$\displaystyle\max_{0\leq n\leq N}|y_n-z(t_n)|\leq S\Big(\max_{0\leq n\leq
r}|y_n-z(t_n)|+\sum_{r\leq n<N}|e_n|\Big).$\hfill 
\par\majfin\par

\bigskip
Si la m\'ethode est d'ordre $p$ avec $|e_n|\leq Ch_nh^p_{\rm max}$, alors on a
$\displaystyle\sum_{r\leq n<N}|e_n|\leq CTh^p_{\rm max}$ car $\displaystyle\sum h_n=T$. Pour la
phase d'initialisation, il convient donc de choisir une m\'ethode conduisant \`a une erreur
d'initialisation $\displaystyle\max_{0\leq n\leq r}|y_n-z(t_n)|$ de l'ordre de $h^p_{\rm max}$
au plus. Ceci conduit \`a choisir une m\'ethode d'initialisation d'ordre $\geq p-1$. Il est
toutefois pr\'ef\'erable de choisir une m\'ethode d'ordre $\geq p$, car l'erreur
d'initialisation est alors born\'ee par $C'h^{p+1}_{\rm max}$ et donc n\'egligeable.

\bigskip
\state Condition n\'ecessaire de stabilit\'e -- |  On cherche \`a 
d\'eterminer les conditions n\'ecessaires \`a la stabilit\'e de la m\'ethode 
(M) d\'ecrite au \S 1.1. Pour cela, on consid\`ere l'\'equation 
diff\'erentielle la plus simple qui soit~:
$$
y'=0.\leqno{\rm (E)}
$$
La suite $y_n$ est alors d\'efinie par
$$
y_{n+1}=\sum_{0\leq i\leq r}\alpha_i y_{n-i},\qquad n\geq r.
$$
L'ensemble des suites v\'erifiant cette relation de r\'ecurrence
est un espace vectoriel de dimension $r+1$, car chaque suite est
d\'efinie de mani\`ere unique par la donn\'ee de
$(y_0,y_1,\dots,y_r)$. On a de mani\`ere \'evidente des solutions
particuli\`eres
$$
y_n=\lambda^n,
$$
o\`u $\lambda$ est racine du polyn\^ome caract\'eristique
$$
\lambda^{r+1}-\alpha_0\lambda^r-\alpha_1\lambda^{r-1}-\dots-\alpha_r=0.
$$
Soient $\lambda_j$ les racines complexes de ce polyn\^ome, et $m_j$
les multiplicit\'es correspondantes. On sait (par une th\'eorie en
tout point analogue \`a celle des \'equations diff\'erentielles
lin\'eaires \`a coefficients constants), qu'une base de l'espace
vectoriel des suites consid\'er\'ees est form\'ee des suites
$$
n\mapsto n^q\lambda_j^n,\qquad 0\leq q<m_j.
$$
Consid\'erons maintenant la suite $y_n\equiv 0$ et la suite
$\widetilde y_n=\varepsilon\lambda_j^n$, $0\leq n\leq N$ avec
$\varepsilon>0$ petit (on a ici $\varepsilon_n=0$, seule l'erreur
d'initialisation intervient). Si la m\'ethode (M) est stable, on doit
avoir
$$
|\widetilde y_N-y_N|=\varepsilon|\lambda_j|^N\leq S\max_{0\leq n\leq r}
\varepsilon|\lambda_j|^n,
$$
ce qui \'equivaut \`a
$$
|\lambda_j|^N\leq S\ \max(1,|\lambda_j|^r).
$$
Si l'on fait tendre $h$ vers 0 et $N=\dfrac{T}{h}$ vers $+\infty$,
ceci n'est possible que pour $|\lambda_j|\leq 1$.

Supposons maintenant que le polyn\^ome caract\'eristique admette une
racine $\lambda_j$ de module $|\lambda_j|=1$ et de multiplicit\'e
$m_j\geq 2$. En regardant la suite $\widetilde y_n=n\lambda_j^n$\break
on trouve
$$
|\widetilde y_N-y_N|=\varepsilon N,\qquad \max_{0\leq n\leq r}
|\widetilde y_n-y_n|=\varepsilon r
$$
ce qui contredit la stabilit\'e. On obtient donc la condition
n\'ecessaire suivante~: $|\lambda_j|\leq 1$ pour tout $j$, et si
$|\lambda_j|=1$ alors cette racine doit \^etre simple.

\bigskip
\state Condition suffisante de stabilit\'e* -- | On va voir que la
condition n\'ecessaire qui vient d'\^etre trouv\'ee est en fait suffisante.

\bigskip
\majdeb
\state Th\'eor\`eme -- |  {\it On suppose que $f(t,y)$ est $k$-lipschitzienne
en $y$. Alors la m\'ethode {\rm (M)} est stable si et seulement si
$$
\lambda^{r+1}-\alpha_0\lambda^r-\dots-\alpha_r=0
$$
a toutes ses racines de module $\leq 1$ et si les racines de module
1 sont simples.}
\par\majfin\par

\bigskip
\noindent{\it Remarque.} On  observera que $\lambda=1$ est toujours racine
d\`es que la m\'ethode est consistante, puisqu'alors $\alpha_0+\dots+
\alpha_r=1$.

\bigskip
\noindent{\prem D\'emonstration.} Soient deux suites $y_n,\widetilde y_n$
telles que
$$
\left.\matrix{
\displaystyle y_{n+1}=\sum_{0\leq i\leq r}\alpha_iy_{n-i}+h\beta_if_{n-i}
\hfill\cr
\noalign{\vskip3pt}
\displaystyle\widetilde y_{n+1}=\sum_{0\leq i\leq r}
\alpha_i\widetilde y_{n-i}+h\beta_i\widetilde f_{n-i}+\varepsilon_n\hfill\cr}
\right\}\quad r\leq n<N
$$
avec $\widetilde f_n=f(t_n,\widetilde y_n)$. Posons
$\theta_n=\widetilde y_n-y_n$. Il vient $|\widetilde f_n-f_n|\leq
k|\theta_n|$ et
$$
\theta_{n+1}-\sum_{0\leq i\leq r}\alpha_i\theta_{n-i}=h\sum_{0\leq i\leq r}\beta_i(\widetilde
f_{n-i}-f_{n-i})+\varepsilon_n. 
$$
Posons $\sigma_n=\theta_n-\alpha_0\theta_{n-1}-\dots-\alpha_r\theta_{n-r-1}$. On a donc
$$
|\sigma_{n+1}|\leq kh\sum_{0\leq i \leq r}|\beta_i||\theta_{n-i}|+|\varepsilon_n|,\quad r\leq
n\leq N.\eqno(*) 
$$
Pour exploiter cette in\'egalit\'e, on cherche \`a majorer
$|\theta_{n+1}|$ en fonction des $|\sigma_i|$. Pour cela, on observe
que la relation de d\'efinition de $\sigma_n$ \'equivaut \`a
l'\'egalit\'e formelle
$$
\sum\sigma_nX^n=(\sum\theta_nX^n)(1-\alpha_0X-\dots-\alpha_rX^{r+1})
$$
avec la convention $\sigma_n=0$, $\theta_n=0$ pour $n<0$. On a donc
inversement
$$
\sum\theta_n X^n=\frac{1}{1-\alpha_0 X-\dots-\alpha_rX^{r+1}}\ \sum\sigma_nX^n.
$$
Consid\'erons le d\'eveloppement en s\'erie en $X=0$~:
$$
\frac{1}{1-\alpha_0X-\dots-X^{r+1}}=\sum\gamma_nX^n.
$$

\bigskip
\majdeb
\state Lemme -- |  {\it Sous les hypoth\`eses du th\'eor\`eme, les 
coefficients $\gamma_n$ sont born\'es.} 
\par\majfin\par

\bigskip
En effet, si les racines du polyn\^ome caract\'eristique sont les complexes $\lambda_j$ de
multiplicit\'e $m_j$, on a
$$
1-\alpha_0 X-\dots-\alpha_rX^{r+1}=\prod_j (1-\lambda_j X)^{m_j}.
$$
Il existe par cons\'equent une d\'ecomposition en \'el\'ements simples
$$
\frac{1}{1-\alpha_0 X-\dots-\alpha_rX^{r+1}}=\sum_{j,q\leq m_j}\frac{c_{jq}}{(1-\lambda_jX)^q}.
$$
Par r\'ecurrence sur $q$ et par d\'erivation, on v\'erifie facilement que 
$$
\frac{1}{(1-\lambda_jX)^q}=\sum_{n=0}^{+\infty}\frac{(n+1)(n+2)\dots(n+q-1)}{(q-1)!}\
\lambda_j^n X^n. 
$$
Les coefficients de cette derni\`ere s\'erie admettent l'\'equivalent
$n^{q-1}\lambda^n_j/(q-1)!$ et sont donc born\'es si et seulement si $|\lambda_j|<1$ ou bien
$|\lambda_j|=1$ et $q=1$.\hfill\blacksquare

\medskip
Notons $\displaystyle\Gamma=\sup_{n\in\bN}|\gamma_n|<+\infty$. Comme $\gamma_0=1$, on a toujours
$\Gamma\geq 1$. La relation
$$
\sum\theta_nX^n=\sum\gamma_nX^n\sum\sigma_nX^n
$$
\'equivaut \`a $\theta_n=\gamma_0\sigma_n+\gamma_1\sigma_{n-1}+\dots+\gamma_n\sigma_0$, d'o\`u
$$
|\theta_n|\leq\Gamma(|\sigma_0|+|\sigma_1|+\dots+|\sigma_{n}|).\eqno(**)
$$
En combinant $(*)$ et $(**)$ il vient
$$
\eqalign{
|\theta_{n+1}|&\leq \Gamma\sum_{0\leq j\leq n+1}|\sigma_j|\leq\Big[\sum_{r\leq j\leq
n}|\sigma_{j+1}|+\sum_{0\leq j\leq r}|\sigma_j|\Big]\cr 
|\theta_{n+1}|&\leq\Gamma\Big[\sum_{r\leq j\leq n}\Big(kh\sum_{0\leq i\leq
r}|\beta_i||\theta_{j-i}|+|\varepsilon_j|\Big)+\sum_{0\leq j\leq r}|\sigma_j|\Big].\cr 
}
$$
Or $\displaystyle\sum_{r\leq j\leq n}\ \sum_{0\leq i\leq
r}|\beta_i||\theta_{j-i}|\leq\Big(\sum_{0\leq i\leq r}|\beta_i|\Big)\Big(\sum_{0\leq j\leq
n}|\theta_j|\Big)$

\noindent et la relation de d\'efinition
$\sigma_j=\theta_j-\alpha_0\theta_{j-1}-\dots-\alpha_r\theta_{j-r-1}$ donne par ailleurs
$$
\sum_{0\leq j\leq r}|\sigma_j|\leq\Big(1+\sum_{0\leq i\leq
r}|\alpha_i|\Big)\Big(|\theta_0|+\dots+|\theta_r|\Big). 
$$
On obtient donc
$$
\eqalign{
|\theta_{n+1}|&\leq \Gamma kh\sum_{0\leq i\leq r}|\beta_i|(|\theta_0|+\dots+|\theta_n|)\cr
&+\Gamma\left[\sum_{r\leq j\leq n}|\varepsilon_j|+\Big(1+\sum_{0\leq i\leq
r}|\alpha_i|\Big)\sum_{0\leq i\leq r}|\theta_i|\right].\cr 
}
$$
Posons $\delta_n=|\theta_0|+\dots+|\theta_n|$ et
$$
\eqalign{
\Lambda&=\Gamma k\sum_{0\leq i\leq r}|\beta_i|,\cr
\eta_n&=\Gamma\left[\sum_{r\leq j\leq n}|\varepsilon_j|+\Big(1+\sum_{0\leq i\leq
r}|\alpha_i|\Big)\sum_{0\leq i\leq r}|\theta_i|\right]. 
}
$$
La derni\`ere in\'egalit\'e s'\'ecrit maintenant
$$
\delta_{n+1}-\delta_n\leq\Lambda h\delta_n+\eta_n\Leftrightarrow\delta_{n+1} \leq (1+\Lambda
h)\delta_n+\eta_n 
$$
et le lemme de Gronwall discret (cf.\ VIII 2.3) donne
$$
\delta_n\leq e^{\Lambda nh}\delta_0+\sum_{0\leq j\leq n-1}e^{\Lambda(n-1-j)h}\eta_j.
$$
Or pour $j\leq n-1$ on a $\eta_j\leq \eta_n$ et $\delta_0=|\theta_0|\leq \eta_n$. On en d\'eduit
$$
\delta_n\leq \eta_n(1+e^{\Lambda h}+\dots+e^{\Lambda
nh})=\frac{e^{\Lambda(n+1)h}-1}{e^{\Lambda h}-1}\ \eta_n.
$$
Cette in\'egalit\'e entra\^\i ne aussit\^ot une majoration de $\theta_n$~:
$$
\eqalign{
|\theta_n|&=\delta_n-\delta_{n-1}\leq\Lambda h\delta_{n-1}+\eta_{n-1}\cr
&=\Big(\Lambda h\ \frac{e^{\Lambda nh}-1}{e^{\Lambda h}-1}+1\Big)\ \eta_{n-1}.\cr
}
$$
Comme $e^{\Lambda h}-1\geq \Lambda h$, il vient $|\theta_n|\leq e^{\Lambda nh}\eta_{n-1}$, soit
$$
|\theta_n|\leq \Gamma e^{\Lambda nh}\Big[\Big(1+\sum_{0\leq i\leq
r}|\alpha_i|\Big)\sum_{0\leq i\leq r}|\theta_i|+\sum_{r\leq j\leq n-1}|\varepsilon_j|\Big],
$$

\bigskip
\majdeb
\hfil$\max_{0\leq n\leq N}|\theta_n|\leq S'\Big[\Big(1+\sum_{0\leq n\leq
r}|\alpha_i|\Big)\sum_{0\leq i\leq r}|\theta_n|+\sum_{r\leq n\leq N}|\varepsilon_n|\Big]$\hfill
\par\majfin\par

\bigskip
\noindent avec $S'=\Gamma e^{\Lambda T}$, c'est-\`a-dire

\bigskip
\majdeb
\hfil$S'=\Gamma e^{\Gamma kT\sum|\beta_i|}$\quad o\`u\quad $\Gamma=\sup\,|\gamma_n|$.\hfill
\par\majfin\par

\bigskip
Si l'erreur initiale $\displaystyle\max_{0\leq n\leq r}|\theta_n|$ est n\'egligeable (comme
c'est souvent le cas) on prendra $S=S'$, sinon on peut prendre
$$
S=(1+r)\Big(1+\sum_{0\leq i\leq r}|\alpha_i|\Big)S'.
$$

\bigskip
\state Remarque -- | 
Pour une fonction $f(t,y)$ de constante de Lipschitz $k$ et pour une
dur\'ee d'int\'egration $T$ fix\'ees, la stabilit\'e de la m\'ethode
d\'epend essentiellement de la grandeur de la constante
$\Gamma\sum|\beta_i|$. On a donc int\'er\^et \`a choisir une m\'ethode
pour laquelle cette constante soit la plus petite possible.

\subsection 1.3. Exemples|

\noindent{\bf $\bullet$ M\'ethode de Nystr\"om}
\bigskip

C'est la m\'ethode \`a 2 pas d\'efinie par
$$
y_{n+1}=y_{n-1}+2hf_n,\quad n\geq 1.
$$
On a ici $\alpha_0=0$, $\alpha_1=1$, $\beta_0=2$, $\beta_1=0$. Le principe de cette m\'ethode est
analogue \`a celui de la m\'ethode du point milieu~:

\InsertFig 18.000 67.000 {
1 mm unit
  0.000  10.000 moveto 
 90.000   0.000 vector 
 10.000   5.000 moveto 
 53.000  90.000 vector 
  0.340 setlinewidth 
 10.000  36.000 moveto 
  0.700   0.000 indent 
 10.000  25.000 moveto 
  0.700   0.000 indent 
 30.000  10.000 moveto 
  0.700  90.000 indent 
 45.000  10.000 moveto 
  0.700  90.000 indent 
 60.000  10.000 moveto 
  0.700  90.000 indent 
 20.000  30.000 moveto 
[  20.000  30.000   30.000  25.000   45.000  26.000   60.000  36.000  
 70.000  54.000  ] curve 
stroke 
  0.113 setlinewidth 
 30.000  25.000 moveto 
  0.700 disk
 60.000  36.000 moveto 
 45.000  26.000 moveto 
  0.700 disk 
 60.000  36.000 moveto 
  0.700 disk 
 25.000  20.000 moveto 
 65.000  32.000 lineto 
stroke 
 25.000  23.500 moveto 
 65.000  35.500 lineto 
stroke 
 60.000  34.000 moveto 
  0.700 disk 
 45.000   3.000 moveto 
 15.000   0.000 vector 
 15.000 180.000 vector 
 10.000  25.000 moveto 
[ 2.000  1.000 ] 0 setdash
 30.000  25.000 moveto 
 60.000  36.000 lineto 
stroke 
[   0.400   1.000 ] 0 setdash 
 11.500  25.000 moveto 
 30.000  25.000 lineto 
stroke 
 30.000  11.500 moveto 
 30.000  25.000 lineto 
stroke 
 45.000  11.500 moveto 
 45.000  26.000 lineto 
stroke 
 60.000  11.500 moveto 
 60.000  36.000 lineto 
stroke 
 11.500  36.000 moveto 
 60.000  36.000 lineto 
stroke}
\LabelTeX   28.000   6.200 $t_{n-1}$\ELTX
\LabelTeX   44.200   6.200 $t_n$\ELTX
\LabelTeX   58.000   6.200 $t_{n+1}$\ELTX
\LabelTeX   88.000   6.200 $t_n$\ELTX
\LabelTeX    6.500  56.000 $y$\ELTX
\LabelTeX    1.500  35.000 $y_{n+1}$\ELTX
\LabelTeX    1.500  24.000 $y_{n-1}$\ELTX
\LabelTeX   69.000  48.000 $f$\ELTX
\LabelTeX   66.500  31.000 $\hbox{\rm pente}~f_n$\ELTX
\LabelTeX   43.700  -1.000 $2h$\ELTX
\EndFig
\vskip8mm

\noindent on approxime la pente $\dfrac{1}{2h}\ (y_{n+1}-y_{n-1})$ de la corde par la pente de
la tangente au point milieu $t_n$. Les calculs du \S\ 1.1 donnent
$$
e_n=\frac{h^3}{3!}\ z^{(3)}(t_n)\cdot 2+O(h^4)=\frac{h^3}{3}\ f^{[2]}(t_n,y_n)+O(h^4).
$$
La m\'ethode de Nystr\"om est donc d'ordre 2. Pour la phase d'initialisation, on choisira une
m\'ethode \`a 1 pas d'ordre 2, par exemple la m\'ethode du point milieu~:
$$
\eqalign{
&f_0=f(t_0,y_0)~;\cr
&y_{1/2}=y_0+\frac{h}{2}\ f_0~;\quad t_{1/2}=t_0+\frac{h}{2}~;\quad
f_{1/2}=f(t_{1/2},y_{1/2})~;\cr
&y_1=y_0+hf_{1/2}~;\quad t_1=t_0+h~;\quad f_1=f(t_1,y_1).\cr 
}
$$
Le polyn\^ome caract\'eristique est $\lambda^2-1$. D'apr\`es le \S\ 1.2 la m\'ethode est stable
et
$$
\frac{1}{1-\alpha_0X-\alpha_1X^2}=\frac{1}{1-X^2}=\sum X^{2n}.
$$
On a donc $\Gamma=1$, $|\beta_0|+|\beta_1|=2$, d'o\`u la constante de stabilit\'e
$$
S'=e^{2kt}.
$$
On observe n\'eanmoins que le polyn\^ome caract\'eristique admet la racine $\lambda=-1$ situ\'ee
sur le cercle limite de stabilit\'e $|\lambda|=1$, en plus de la racine oblig\'ee $\lambda=1$. 
Ceci laisse suspecter que la stabilit\'e n'est peut-\^etre pas tr\`es bonne. Pour mettre en
\'evidence ce ph\'enom\`ene, regardons le cas de l'\'equation diff\'erentielle 
$$
y'=-y\leqno{\rm (E)}
$$
avec donn\'ee initiale $t_0=0$, $y_0=1$. La solution exacte est donn\'ee par
$$
z(t)=e^{-t},\quad z(t_n)=e^{-nh}.
$$
La suite $y_n$ sera donn\'ee par la relation de r\'ecurrence
$$
y_{n+1}=y_{n-1}-2hy_n,\qquad n\geq 1\eqno(*)
$$
(car ici $f_n=-y_n)$, avec valeurs initiales~:
$$
\eqalign{
&y_0=1,\qquad\qquad f_0=-1,\cr
&y_{1/2}=1-\frac{h}{2},\quad f_{1/2}=-1+\frac{h}{2},\cr
&y_1=1-h+\frac{h^2}{2}.\cr
}
$$
La solution g\'en\'erale de $(*)$ peut s'\'ecrire
$$
y_n=c_1\lambda_1^n+c_2\lambda_2^n
$$
o\`u $\lambda_1,\lambda_2$ sont les racines de l'\'equation
$$
\lambda^2+2hk-1=0,
$$
\`a savoir $\lambda_1=-h+\sqrt{1+h^2}$, $\lambda_2=-h-\sqrt{1+h^2}$.
Les constantes $c_1,c_2$ sont d\'etermin\'ees par
$$
\cases{
y_0=c_1+c_2=1\hfill\cr
\noalign{\vskip3pt}
y_1=c_1\lambda_1+c_2\lambda_2=1-h+\dfrac{h^2}{2}.\hfill\cr
}
$$
On obtient donc
$$
\eqalign{
c_1&=\frac{1-h+\dfrac{h^2}{2}-\lambda_2}{\lambda_1-\lambda_2}=\frac{1+\dfrac{h^2}{2}+
\sqrt{1+h^2}}{2\sqrt{1+h^2}},\cr
c_2&=\frac{\lambda_1-\Big(1-h+\dfrac{h^2}{2}\Big)}{\lambda_1-\lambda_2}=\frac{\sqrt{1+h^2}
-\Big(1+\dfrac{h^2}{2}\Big)}{2\sqrt{1+h^2}}.\cr
}
$$
Comme $\sqrt{1+h^2}=1+\dfrac{h^2}{2}-\dfrac{h^4}{8}+O(h^6)$, on voit que
$$
\eqalign{
&c_1=1+O(h^4)\cr
&c_2=-\frac{1}{16}\ h^4+O(h^6)\cr
}
$$
Par ailleurs $\lambda_1=1-h+\dfrac{h^2}{2}+O(h^4)=e^{-h}+O(h^3)$, 

\noindent tandis que\kern8pt $\lambda_2=-\Big(1+h+\dfrac{h^2}{2}+O(h^4)\Big)=-e^h+O(h^3)$.

\noindent On voit donc que $y_n$ est somme de deux termes
$$
c_1\lambda_1^n\simeq e^{-nh},\qquad c_2\lambda_2^n\simeq -\dfrac{1}{16}\ h^4(-1)^ne^{nh}.
$$
Le premier de ces termes approxime bien la solution exacte, mais le second est un terme de
perturbation qui diverge quand $n\rightarrow+\infty$, bien qu'il soit n\'egligeable quand le
temps $t_n=nh$ n'est pas trop grand. Si la dur\'ee d'int\'egration est trop longue (plus
pr\'ecis\'ement, si $h^4e^T$ n'est plus n\'egligeable),  on va obtenir un trac\'e de la forme
ci-dessous.

\InsertFig -30.000 118.000 {
1 mm unit
 40.000  65.000 moveto 
110.000   0.000 vector 
 50.000  50.000 moveto 
 55.000  90.000 vector
 50.000  95.000 moveto 
[  50.000  95.000   58.000  89.600   66.000  84.000   74.000  81.600  
 82.000  76.900   90.000  76.300   98.000  71.800  106.000  73.300  
114.000  67.300  122.000  75.300  130.000  57.300  ] polygon 
stroke 
 50.000  95.000 moveto 
  0.700 disk 
 58.000  89.600 moveto 
  0.700 disk 
 66.000  84.000 moveto 
  0.700 disk 
 74.000  81.600 moveto 
  0.700 disk 
 82.000  76.900 moveto 
  0.700 disk 
 90.000  76.300 moveto 
  0.700 disk 
 98.000  71.800 moveto 
  0.700 disk 
106.000  73.300 moveto 
  0.700 disk 
114.000  67.300 moveto 
  0.700 disk 
122.000  75.300 moveto 
  0.700 disk 
130.000  57.300 moveto 
  0.700 disk 
  0.340 setlinewidth
 50.000  95.000 moveto 
[  50.000  95.000    50.444  94.647    50.891  94.295    51.340  93.947   
 51.792  93.600    52.246  93.256    52.702  92.914    53.161  92.575   
 53.622  92.238    54.085  91.904    54.551  91.572    55.019  91.242   
 55.490  90.915    55.963  90.590    56.438  90.267    56.916  89.947   
 57.397  89.630    57.880  89.315    58.365  89.002    58.853  88.692   
 59.343  88.384    59.836  88.079    60.331  87.776    60.829  87.476   
 61.329  87.178    61.831  86.883    62.337  86.590    62.844  86.299   
 63.355  86.012    63.867  85.726    64.383  85.443    64.900  85.163   
 65.421  84.885    65.943  84.610    66.469  84.337    66.996  84.067   
 67.527  83.799    68.060  83.534    68.595  83.271    69.133  83.011   
 69.673  82.753    70.216  82.498    70.762  82.245    71.310  81.995   
 71.860  81.748    72.413  81.502    72.969  81.260    73.527  81.020   
 74.087  80.782    74.651  80.547    75.216  80.314    75.784  80.084   
 76.355  79.856    76.928  79.631    77.503  79.408    78.081  79.188   
 78.661  78.970    79.244  78.754    79.830  78.541    80.417  78.331   
 81.008  78.122    81.600  77.917    82.195  77.713    82.793  77.512   
 83.393  77.314    83.995  77.118    84.599  76.924    85.206  76.732   
 85.816  76.543    86.428  76.357    87.042  76.172    87.658  75.990   
 88.277  75.810    88.898  75.633    89.521  75.457    90.146  75.284   
 90.774  75.114    91.404  74.945    92.037  74.779    92.671  74.615   
 93.308  74.453    93.947  74.293    94.588  74.136    95.231  73.980   
 95.877  73.827    96.524  73.676    97.174  73.527    97.826  73.380   
 98.480  73.235    99.136  73.092    99.794  72.952   100.454  72.813   
101.116  72.676   101.780  72.541   102.446  72.409   103.114  72.278   
103.784  72.149   104.456  72.022   105.130  71.897   105.805  71.774   
106.483  71.652   107.162  71.533   107.844  71.415   108.527  71.300   
109.212  71.186   109.898  71.073   110.587  70.963   111.277  70.854   
111.969  70.747   112.663  70.642   113.358  70.538   114.055  70.436   
114.754  70.336   115.454  70.237   116.156  70.140   116.860  70.044   
117.565  69.950   118.272  69.858   118.980  69.767   119.690  69.678   
120.401  69.590   121.114  69.503   121.828  69.418   122.544  69.335   
123.261  69.253   123.979  69.172   124.699  69.093   125.421  69.015   
126.143  68.938   126.867  68.863   127.592  68.789   128.319  68.716   
129.047  68.645   129.776  68.575   130.506  68.506   131.238  68.438   
131.971  68.371   132.705  68.306   133.440  68.242   134.176  68.179   
134.914  68.117   135.652  68.056   136.392  67.996   137.133  67.938   
137.875  67.880   138.618  67.824   139.362  67.768   140.107  67.714   
140.853  67.660   141.600  67.608   142.348  67.556   143.097  67.506   
143.847  67.456   144.597  67.408   145.349  67.360   ] scurve 
106.000  65.000 moveto 
  0.700  90.000 indent 
130.000  65.000 moveto 
  0.700  90.000 indent}
\LabelTeX   46.500  61.200 $0$\ELTX
\LabelTeX   46.500  94.000 $1$\ELTX
\LabelTeX   46.500 103.000 $y$\ELTX
\LabelTeX  104.800  61.200 $t_n$\ELTX
\LabelTeX  130.000  61.200 $T$\ELTX
\LabelTeX  148.500  61.200 $t$\ELTX
\LabelTeX  105.000  76.500 $y_n$\ELTX
\LabelTeX  135.000  70.500 $y=e^{-t}$\ELTX
\EndFig
\vskip-44mm

\noindent{\bf $\bullet$ M\'ethode de Milne}
\medskip

C'est la m\'ethode \`a 4 pas d\'efinie par
$$
y_{n+1}=y_{n-3}+h\Big(\frac{8}{3}\ f_n-\frac{4}{3}\ f_{n-1}+\frac{8}{3}\ f_{n-2}\Big).
$$
On v\'erifie qu'elle est stable, d'ordre 4, mais comme pour la m\'ethode de Nystr\"om la
stabilit\'e n'est pas tr\`es bonne \`a cause des racines de module 1 du polyn\^ome
caract\'eristique $\lambda^4-1=0$.

\section 2. M\'ethodes d'Adams-Bashforth|

\subsection 2.1. Description|

\noindent On ne suppose plus ici que le pas $h_n$ soit n\'ecessairement constant. Si $z$ est
une solution exacte de l'\'equation, on \'ecrit
$$
z(t_{n+1})=z(t_n)+\int^{t_{n+1}}_{t_n}f(t,z(t))dt.
$$
Supposons que pour $0\leq i\leq r$ on ait d\'ej\`a calcul\'e les points $z(t_{n-i})$ et les pentes
$f_{n-i}=f(t_{n-i}$, $z(t_{n-i}))$. 

L'id\'ee de la m\'ethode est d'approximer la fonction $f(t,z(t))$ sur $[t_n,t_{n+1}]$ par son
polyn\^ome d'interpolation aux points $t_n,t_{n-1},\dots,t_{n-r}$. Consid\'erons donc le polyn\^ome
$p_{n,r}(t)$ qui interpole les points $(t_{n-i},f_{n-i})$ pour $0\leq i\leq r$~:
$$
p_{n,r}(t)=\sum_{0\leq i\leq r}f_{n-i}L_{n,i,r}(t),\quad \deg\,(p_{n,r})=r,
$$
o\`u $\displaystyle L_{n,i,r}(t)=\prod_{\scriptstyle 0\leq j\leq r\atop \scriptstyle j\not=i}\
\dfrac{t-t_{n-j}}{t_{n-i}-t_{n-j}}$. On \'ecrit maintenant~:
$$
\eqalign{
z(t_{n+1})&=z(t_n)+\int^{t_{n+1}}_{t_n}f(t,z(t))dt\cr
&\simeq z(t_n)+\int^{t_{n+1}}_{t_n}p_{n,r}(t)dt\cr
&=z(t_n)+h_n\sum_{0\leq i\leq r}b_{n,i,r}f_{n,i}\cr
}
$$
avec
$$
b_{n,i,r}=\frac{1}{h_n}\int^{t_{n+1}}_{t_n}L_{n,i,r}(t)dt.
$$
L'algorithme de la m\'ethode d'Adams-Bashforth \`a $r+1$ pas (en abr\'eg\'e $AB_{r+1}$) va donc
s'\'ecrire~:
$$
\cases{
\displaystyle y_{n+1}=y_n+h_n\sum_{0\leq i\leq r}b_{n,i,r}f_{n-i},\quad n\geq r,\hfill\cr
\noalign{\vskip3pt}
t_{n+1}=t_n+h_n\hfill\cr
\noalign{\vskip3pt}
f_{n+1}=f(t_{n+1},y_{n+1}).\hfill\cr
}
$$
L'int\'er\^et de cette m\'ethode provient de sa relative simplicit\'e et du fait qu'une seule
\'evaluation de la fonction $f$ est n\'ecessaire \`a chaque \'etape (contrairement aux
m\'ethodes de Runge-Kutta qui en r\'eclamaient plusieurs). Il va en r\'esulter un gain assez
important sur le temps de calcul.


\bigskip
\state Exemples|

\medskip
\noindent$\bullet$ $r=0$~:\quad on a $p_{n,0}(t)=$ constante $=f_n$, 
d'o\`u $AB_1$~:\quad $y_{n+1}=y_n+h_nf_n$. Il s'agit de la m\'ethode d'Euler.

\medskip
\noindent$\bullet$ $r=1$~:  le polyn\^ome $p_{n,1}$ est la fonction affine qui interpole
$(t_n,f_n)$ et $(t_{n-1},f_{n-1})$, d'o\`u les formules 
$$
\eqalign{
p_{n,1}(t)=f_n&+\frac{f_n-f_{n-1}}{t_n-t_{n-1}}\ (t-t_n)\cr
\int^{t_{n+1}}_{t_n}p_{n,1}(t)dt&=f_nh_n+\frac{f_n-f_{n-1}}{h_{n-1}}\ \Big[\frac{1}{2}\
(t-t_n)^2 \Big]^{t_{n+1}}_{t_n}\cr
&=b_n\Big(f_n+\frac{h_n}{2h_{n-1}}\ (f_n-f_{n-1})\Big).\cr
}
$$
L'algorithme s'\'ecrit donc
$$
AB_2\quad\cases{
y_{n+1}=y_n+h_n\Big(f_n+\displaystyle
\dfrac{h_n}{2h_{n-1}}\ (f_n-f_{n-1})\Big)\hfill\cr
\noalign{\vskip3pt}
t_{n+1}=t_n+h_n\hfill\cr
\noalign{\vskip3pt}
f_{n+1}=f(t_n,y_n).\hfill\cr
}
$$
Dans le cas o\`u le pas $h_n=h$ est constant, la formule de r\'ecurrence se r\'eduit \`a
$$
y_{n+1}=y_n+h\Big(\frac{3}{2}\ f_n-\frac{1}{2}\ f_{n-1}\Big).
$$
$\bullet$ De mani\`ere g\'en\'erale, lorsque le pas est constant
les coefficients $b_{n,i,r}$ ne d\'ependent pas de $n$ car la
m\'ethode est invariante par translation. Pour les petites valeurs
de~$r$, les coefficients $b_{i,r}$ correspondants sont donn\'es par le
tableau~:

$$
\vbox{\offinterlineskip\halign{
\hfill#\hfill\quad&
\tv\hfill\quad#\hfill\quad&
\hfill\quad#\hfill\quad&
\hfill\quad#\hfill\quad&
\hfill\quad#\hfill\quad&
\tv\hfill\quad#\hfill\quad\cr
$r$&$b_{0,r}$&$b_{1,r}$&$b_{2,r}$&$b_{3,r}$&$\displaystyle\beta_r=\sum_i|
b_{i,r}|$\cr
\noalign{\hrule}
0&1&&&&1\phantom{\bigg(}\cr
1&$\dfrac{3}{2}$&$-\dfrac{1}{2}$&&&2\phantom{\Bigg(}\cr
2&$\dfrac{23}{12}$&$-\dfrac{16}{12}$&$\dfrac{5}{12}$&&\quad 3,66$\ldots$
\phantom{\Bigg(}\cr
3&$\dfrac{55}{24}$&$-\dfrac{59}{24}$&$\dfrac{37}{24}$&$-\dfrac{9}{24}$&
\quad6,6$\ldots$\phantom{0\Bigg(}
\cr
}}
$$

\bigskip
\state Remarque -- | On a toujours $\sum_{0\leq i\leq r}
b_{n,i,r}=1$, car pour $f_n=\dots=f_{n-r}=1$ on a $p_{n,r}(t)\equiv 1$,
par cons\'equent
$$
\int^{t_{n+1}}_{t_n}p_{n,r}(t)dt=h_n=h_n\sum_{0\leq i\leq r}b_{n,i,r}\cdot 1.
$$
Comme on le verra plus loin, la quantit\'e $\beta_r$ intervient
dans le calcul de la constante de stabilit\'e $S$.

\subsection 2.2.\ Erreur de consistance et ordre de la 
m\'ethode~$AB_{r+1}$\kern-2pt|

\noindent Soit $z$ une solution exacte du probl\`eme de Cauchy. L'erreur de
consistance est donn\'ee par
$$
\eqalign{
e_n&=z(t_{n+1})-y_{n+1}\cr
&=z(t_{n+1})-\Big(z(t_n)+\int^{t_{n+1}}_{t_n}p_{n,r}(t)dt\Big),\cr
e_n&=\int^{t_{n+1}}_{t_n}\Big(z'(t)-p_{n,r}(t)\Big)dt,\cr
}
$$
o\`u $p_{n,r}$ est pr\'ecis\'ement le polyn\^ome d'interpolation de
la fonction $z'(t)=f(t,z(t))$ aux points $t_{n-i}$, $0\leq i\leq r$.
D'apr\`es le th\'eor\`eme de la moyenne, il existe un point
$\theta\in\ ]t_n,t_{n+1}[$ tel que
$$
e_n=h_n\Big(z'(\theta)-p_{n,r}(\theta)\Big).
$$
La formule donnant l'erreur d'interpolation (voir chapitre II, 
\S\ 1.2) implique
$$
z'(\theta)-p_{n,r}(\theta)=\frac{1}{(r+1)!}\ z^{(r+2)}(\xi)\pi_{n,r}(\xi),
$$
o\`u $\xi\in\ ]t_{n-r},t_{n+1}[$ est un point interm\'ediaire entre
$\theta$ et les points $t_{n-i}$, et o\`u
$$
\pi_{n,r}(t)=\prod_{0\leq i\leq r}(t-t_{n-i}).
$$
Si $\xi\in\ ]t_{n-j},t_{n-j+1}[$, $0\leq j\leq r$, on a
l'in\'egalit\'e $|\xi-t_{n-i}|\leq (1+|j-i|)h_{\rm max}$, d'o\`u
$$
\eqalign{
|\pi_{n,r}(\xi)|&\leq h^{r+1}_{\rm max}(1+j)\dots(1+1)1(1+1)\dots(1+r-j)\cr
&=h^{r+1}_{\rm max}(j+1)!(r-j+1)!\leq h^{r+1}_{\rm max}(r+1)!,\cr
}
$$
en majorant 2 par $j+2,\dots,(r-j+1)$ par $(r+1)$. On en d\'eduit
par cons\'equent
$$
|z'(\theta)-p_{n,r}(\theta)|\leq|z^{(r+2)}(\xi)|\,h^{r+1}_{\rm max},
$$
ce qui donne la majoration cherch\'ee de l'erreur de consistance~:

\bigskip
\majdeb
\hfil$|e_n|\leq|z^{(r+2)}(\xi)|h_nh^{r+1}_{\rm max}\leq 
Ch_nh^{r+1}_{\rm max}$\hfill
\par\majfin\par

\bigskip
\noindent avec $\displaystyle C=\max_{t\in[t_0,t_0+T]}|z^{(r+2)}(t)|$.

\bigskip
{\it La m\'ethode d'Adams-Bashforth \`a $r+1$ pas est donc d'ordre $r+1$} 
(le lecteur pourra v\'erifier \`a titre d'exercice que l'ordre n'est pas 
$\geq r+2$ en consid\'erant le cas de la fonction $f(t,y)=t^{r+1}$).

\bigskip
\state Phase d'initialisation -- |  De ce qui pr\'ec\`ede, il r\'esulte 
qu'on choisira une m\'ethode de Runge-Kutta d'ordre $r+1$ (ou $r$ \`a la 
rigueur) pour initialiser les premi\`eres valeurs 
$y_1,\dots,y_r,f_0,f_1,\dots,f_r$.

\subsection 2.3. Stabilit\'e de la m\'ethode $AB_{r+1}$|

\noindent Nous allons d\'emontrer le r\'esultat suivant.

\bigskip
\majdeb
\state Th\'eor\`eme -- |  {\it On suppose que $f(t,y)$ est
$k$-lipschitzienne en $y$ et que les sommes
$\sum_{0\leq i\leq r}|b_{n,i,r}|$ sont major\'ees
ind\'ependamment de $n$ par une constante $\beta_r$.\hfill\break
Alors la m\'ethode d'Adams-Bashforth \`a $r+1$ pas est stable, avec
constante de stabilit\'e
$$
S=\exp\,(\beta_rkT).
$$} 
\vskip-\baselineskip
\majfin
\vskip\baselineskip

\noindent{\prem D\'emonstration.} Soit $\widetilde y_n$ la suite r\'ecurrente
perturb\'ee telle que 
$$
\cases{
\displaystyle\widetilde y_{n+1}=\widetilde y_n+h_n\sum_{0\leq i\leq r}b_{n,i,r}\widetilde
f_{n-i}+\varepsilon_n,\quad r\leq n<N,\hfill\cr 
\noalign{\vskip3pt}
\widetilde f_{n-i}=f(t_{n-i},\widetilde y_{n-i}).\hfill\cr}
$$
Posons $\displaystyle\theta_n=\max_{0\leq i\leq n}|\widetilde y_i-y_i|$. On a
$$
\eqalign{
|\widetilde f_{n-i}-f_{n-i}|&\leq k|\widetilde y_{n-i}-y_{n-i}|\leq k\theta_n,
\cr
|\widetilde y_{n+1}-y_{n+1}|&\leq \theta_n+h_n\sum_{0\leq i\leq r}|b_{n,i,r}|
\cdot k\theta_n+|\varepsilon_n|\cr
&\leq (1+\beta_rkh_n)\theta_n+|\varepsilon_n|.\cr}
$$
Comme $\theta_{n+1}=\max\,(|\widetilde y_{n+1}-y_{n+1}|,\theta_n)$, on en
d\'eduit
$$
\theta_{n+1}\leq (1+\beta_rkh_n)\theta_n+|\varepsilon_n|.
$$
Le lemme de Gronwall implique alors
$$
\theta_N\leq {\rm exp}\,\Big(\beta_rk(t_N-t_r)\Big)\Big(\theta_r+\sum_{r\leq
n<N}|\varepsilon_n|\Big), 
$$
ce qui entra\^\i ne bien la stabilit\'e, avec constante $S=\exp\ (\beta_rkT)$.
D'apr\`es le tableau donn\'e au \S\ 2.1, on voit que la constante $\beta_r$ 
cro\^{\i}t assez vite quand $r$ augmente. La stabilit\'e devient donc de
moins en moins bonne quand le nombre de pas augmente. Cette stabilit\'e
m\'ediocre est un des inconv\'enients les plus s\'erieux des m\'ethodes
d'Adams-Bashforth lorsque $r$ est grand. En pratique, on se limitera le
plus souvent aux cas $r=1$ ou $r=2$.

\bigskip
\state Remarque -- |  L'exemple $AB_2$ donn\'e au \S\ 2.1 montre que les
coefficients $b_{n,i,r}$ ne sont en g\'en\'eral born\'es que si le rapport
$h_n/h_{n-1}$ de 2 pas cons\'ecutifs reste born\'e. Dans la pratique,
il est raisonnable de supposer que
$$
\frac{h_n}{h_{n-1}}\leq\delta
$$
avec disons $\delta\leq 2$. Les formules du \S\ 2.1 donnent dans ce cas~:
$$
|b_{n,i,r}|\leq\max_{t\in[t_n,t_{n+1}]}|L_{n,i,r}(t)|=\prod_{1\leq
j\leq n}\frac{t_{n+1}-t_{n-j}}{|t_{n-i} -t_{n-j}|}
$$
Comme $t_{n+1}-t_{n-j}=h_{n-j}+\dots+h_n\leq
(1+\delta+\dots+\delta^j)h_{n-j}$ et
$$
h_{n-j}\leq\cases{
|t_{n-i}-t_{n-j}|\hfill&\quad si\quad $j>i$\hfill\cr
\noalign{\vskip3pt}
\delta|t_{n-i}-t_{n-j}|\hfill&\quad si\quad $j<i$\hfill\cr
}
$$
Il vient
$$
|b_{n,i,r}|\leq \delta^i\prod_{j\not=i}(1+\delta+\dots+\delta^j)\leq
\prod_{0\leq j\leq r}(1+\delta+\dots+\delta^j).
$$
On a donc l'estimation assez grossi\`ere
$$
\beta_r=\max_n\sum_{0\leq i\leq r}|b_{n,i,r}|\leq (r+1)\prod_{0\leq j\leq r}
(1+\delta+\dots+\delta^j). 
$$

\section 3. M\'ethodes d'Adams-Moulton|

\subsection 3.1. Description|

\noindent L'id\'ee en est la m\^eme que celle des m\'ethodes d'Adams-Bashforth, mais on
approxime ici $f(t,z(t))$ par son polyn\^ome d'interpolation aux points
$t_{n+1},t_n,\dots,t_{n-r}$~; le point $t_{n+1}$ est donc pris en plus. On consid\`ere  le
polyn\^ome $p^*_{n,r}(t)$ de degr\'e $r+1$ qui interpole les points $(t_{n-i},f_{n-i})$ pour
$-1\leq i\leq r$~: 
$$
\leqalignno{
p^*_{n,r}(t)&=\sum_{-1\leq i\leq r}f_{n-i}L^*_{n,i,r}(t),\cr
L^*_{n,i,r}(t)&=\prod_{\scriptstyle -1\leq j\leq r\atop\scriptstyle
j\not=i}\frac{t-t_{n,j}}{t-t_{n,i}}.&\hbox{\rm d'o\`u}\cr 
}
$$
On obtient donc comme au \S 2.2
$$
z(t_{n+1})\simeq z(t_n)+h_n\sum_{-1\leq i\leq r}b^*_{n,i,r}f_{n-i}
$$
avec\hskip7pc $\displaystyle b^*_{n,i,r}=\frac{1}{h_n}\int^{t_{n+1}}_{t_n}L^*_{n,i,r}(t)dt$.
\bigskip
\noindent L'algorithme correspondant $AM_{r+1}$ s'\'ecrit

\bigskip
\majdeb
\hfil$\displaystyle y_{n+1}-h_nb^*_{n,-1,r}f(t_{n+1},y_{n+1})=y_n+h_n\sum_{0\leq i\leq
r}b^*_{n,i,r}f_{n-i}$.\hfill 
\par\majfin\par

\bigskip
On observera ici que $y_{n+1}$ n'est pas donn\'e explicitement en fonction des quantit\'es
$y_n$, $f_{n-i}$ ant\'erieurement calcul\'ees, mais seulement comme solution d'une \'equation
dont la r\'esolution n'est pas {\it a priori} imm\'ediate. Pour cette raison, on dit que la
m\'ethode d'Adams-Moulton est une {\it m\'ethode implicite} (la m\'ethode d'Adams-Bashforth est
dite par opposition explicite). Pour r\'esoudre l'\'equation ci-dessus, on aura recours en
g\'en\'eral \`a une m\'ethode it\'erative. Notons $u_n$ la quantit\'e (explicite)
$$
u_n=y_n+h_n\sum_{0\leq i\leq r}b^*_{n,i,r} f_{n-i}.
$$
Le point $y_{n+1}$ cherch\'e est la solution $x$ de l'\'equation
$$
x=u_n+h_nb^*_{n,-1,r}f(t_{n+1},x).
$$
On va donc calculer la suite it\'er\'ee $x_{p+1}=\varphi(x_p)$ o\`u 
$$
\varphi(x)=u_n+h_nb^*_{n,-1,r}f(t_{n+1},x).
$$
Comme $\varphi'(x)=h_nb^*_{n,-1,r}f'_y(t_{n+1},x)$, l'application $\varphi$ va \^etre
contractante (avec une petite constante de Lipschitz) lorsque $h_n$ est assez petit. Si
$f(t,y)$ est $k$-lipschitzienne en $y$, il suffit que $h_n<\dfrac{1}{|b^*_{n,-1,r}|k}$ pour avoir
convergence. La solution $y_{n+1}$ est alors unique d'apr\`es le th\'eor\`eme du point fixe, et
l'algorithme it\'eratif s'\'ecrit
$$
\cases{
F_p=f(t_{n+1},x_p),\hfill\cr
\noalign{\vskip3pt}
x_{p+1}=u_n+h_nb_{n,-1,r}F_p.\hfill\cr
}
$$
On choisira une valeur initiale $x_0$ qui soit une approximation de
$y_{n+1}$ (la meilleure possible~!), par exemple la valeur donn\'ee par
la m\'ethode d'Adams-Bashforth~:
$$
x_0=y_n+h_n\sum_{0\leq i\leq r}b_{n,i,r}f_{n-i}.
$$
On arr\^ete l'it\'eration pour $|x_{p+1}-x_p|\leq 10^{-10}$ (par exemple) et on prend
$$
\cases{
y_{n+1}=&derni\`ere valeur $x_{p+1}$ calcul\'ee\hfill\cr
\noalign{\vskip3pt}
f_{n+1}=&$f(t_{n+1},y_{n+1})$\quad (ou par \'economie $=F_p$)\hfill\cr
}
$$


\bigskip
\state Exemples|

\smallskip
\noindent$\bullet$ $r=0$~:\quad le polyn\^ome $p^*_{n,0}$ est le polyn\^ome
de degr\'e 1 qui interpole $(t_{n+1},f_{n+1})$ et $(t_n,f_n)$, soit
$$
\displaylines{
p^*_{n,0}(t)=f_n+\frac{f_{n+1}-f_n}{h_n}\ (t-t_n)~;\cr
\int^{t_{n+1}}_{t_n}p^*_{n,0}(t)dt=h_n\Big(\frac{1}{2}\ f_{n+1}+\frac{1}{2}\ 
f_n\Big).\cr}
$$
On obtient ainsi la m\'ethode dite des trap\`ezes (ou m\'ethode de
Crank-Nicolson)~:
$$
y_{n+1}=y_n+h_n\Big(\frac{1}{2}\ f_{n+1}+\frac{1}{2}\ f_n)
$$
ou encore
$$
y_{n+1}-\frac{1}{2}\ h_nf(t_{n+1},y_{n+1})=y_n+\frac{1}{2}\ h_nf_n.
$$

\smallskip
\noindent$\bullet$ $r=1$~:\quad le polyn\^ome $p^*_{n,1}$ interpole les points
$(t_{n+1},f_{n+1})$, $(t_n,f_n)$, $(t_{n-1},f_{n-1})$, d'o\`u les formules
$$
p^*_{n,1}(t)=f_{n+1}\,\frac{(t-t_n)(t-t_{n-1})}{h_n(h_n+h_{n-1})}-f_n\,
\frac{(t-t_{n-1})(t-t_{n-1})}{h_nh_{n-1}}+f_{n-1}\,
\frac{(t-t_n)(t-t_{n+1})}{h_{n-1}(h_n+h_{n-1})},  
$$
$$
\eqalign{
y_{n+1}&=y_n\int^{t_{n+1}}_{t_n}p_{n,1}^*(t)dt,\cr
y_{n+1}&=y_n+h_n\left[\frac{2h_n+3h_{n-1}}{6(h_n+h_{n-1})}\,
f_{n+1}+\frac{3h_{n-1}+h_n}{6h_{n-1}}\,f_n-\frac{h^2_n}{6h_{n-1}(h_n+h_{n-1})
}\,f_{n-1}\right]\kern-1pt.\cr 
}
$$
Dans le cas o\`u le pas  $h_n=h$ est constant, cette formule se r\'eduit \`a
$$
y_{n+1}=y_n+h\Big(\frac{5}{12}\ f_{n+1}+\frac{8}{12}\ f_n-\frac{1}{12}\ f_{n-1}\Big).
$$

\smallskip \noindent$\bullet$ De mani\`ere g\'en\'erale, les
coefficients $b^*_{n,i,r}$ sont des nombres $b^*_{i,r}$ ind\'ependants
de $n$ lorsque le pas est constant. On a la table suivante~:
$$
\vbox{\offinterlineskip\halign{
\hfill#\hfill\quad&
\tv\hfill\quad#\hfill\quad&
\hfill\quad#\hfill\quad&
\hfill\quad#\hfill\quad&
\hfill\quad#\hfill\quad&
\hfill\quad#\hfill\quad&
\tv\quad#\hfill\quad&
\tv\quad#\hfill\cr
$r$&$b^*_{\rm-1,r}$&$b^*_{0,r}$&$b_{1,r}^*$&$b_{2,r}^*$&$b_{3,r}^*$&$
\displaystyle
\beta^*_r=\sum_i|b^*_{i,r}|$&$\beta_{r+1}$\cr
\noalign{\hrule}
$0$&$\dfrac{1}{2}$&$\dfrac{1}{2}$&&&&1&2\phantom{\Bigg(}\cr
$1$&$\dfrac{5}{12}$&$\dfrac{8}{12}$&$-\dfrac{1}{12}$&&&1,16$\ldots$
&3,66$\ldots$\phantom{\Bigg(}\cr
$2$&$\dfrac{9}{24}$&$\dfrac{19}{24}$&$-\dfrac{5}{24}$&$\dfrac{1}{24}$&
&1,41$\ldots$&6,66$\ldots$\phantom{\Bigg(}\cr
$3$&$\dfrac{251}{720}$&$\dfrac{646}{720}$&$-\dfrac{264}{720}$&$\dfrac{106}{720}$&
$-\dfrac{19}{720}$&1,78$\ldots$&12,64$\ldots$\phantom{\Bigg(}\cr
}}
$$
\noindent Les coefficients $b^*_{n,i,r}$ v\'erifient toujours 
$\displaystyle\sum_{-1\leq i\leq r}b^*_{n,i,r}=1$.

\subsection 3.2.\ Erreur de consistance et ordre de la 
m\'ethode~$AM_{r+1}$\kern-5.1pt|

\noindent Soit $z$ une solution exacte du probl\`eme de Cauchy. Sous l'hypoth\`ese
$y_{n-i}=z(t_{n-i})$, $0\leq i\leq n$, on a
$$
\eqalign{
e_n&=z(t_{n+1})-y_{n+1}\cr
&=z(t_{n+1})-\Big[z(t_n)+h_n\sum_{0\leq i\leq
r}b^*_{n,i,r}f(t_{n-i},z(t_{n-i}))+h_nb^*_{n,-1,r}f(t_{n+1},y_{n+1})\Big]\cr
&=z(t_{n+1})-\Big[z(t_n)+h_n\sum_{-1\leq i\leq r}b^*_{n,i,r}f(t_{n-i},z(t_{n-i}))\Big]\cr
&\qquad\qquad\kern6pt+h_n b^*_{n,-1,r}[f(t_{n+1},z(t_{n+1}))-f(t_{n+1},y_{n+1})],\cr 
e_n&=\int^{t_{n+1}}_{t_n}(z'(t)-p^*_{n,r}(t))dt+h_nb^*_{n,-1,r}\Big[f(t_{n+1},z(t_{n+1}))-f(t_{n+1},
y_{n+1})\Big].\cr
}
$$
Supposons que $f(t,y)$ soit lipschitzienne de rapport $k$ en $y$. Alors il vient
$$
\eqalign{
|e_n|&\leq\Big|\int^{t_{n+1}}_{t_n}(z'(t)-p^*_{n,r}(t))dt\Big|+h_nb^*_{n,-1,r}k|e_n|,\cr
|e_n|&\leq\frac{1}{1-h_nb^*_{n,-1,r}k}\ \Big|\int^{t_{n+1}}_{t_n}(z'(t)-p^*_{n,r}(t))dt\Big|.\cr
}
$$
Quand le pas $h_n$ est suffisamment petit, on a donc
$$
|e_n|=\Big|\int^{t_{n+1}}_{t_n}(z'(t)-p^*_{n,r}(t))dt\Big|(1+O(h_n))
$$
comme dans la m\'ethode d'Adams-Bashforth. Par ailleurs la formule de la moyenne donne
$$
\eqalign{
&\int^{t_{n+1}}_{t_n}(z'(t)-p^*_{n,r}(t))dt=h_n(z'(\theta)-p^*_{n,r}(\theta)),\quad
\theta\in\ ]t_n,t_{n+1}[,\cr 
&z'(\theta)-p^*_{n,r}(\theta)=\frac{1}{(r+2)!}\ z^{(r+3)}(\xi)\pi^*_{n,r}(\xi),\quad
\xi\in\ ]t_{n,r},t_{n+1}[.\cr 
}
$$
o\`u $\displaystyle\pi^*_{n,r}(t)=\prod_{-1\leq i\leq r}(t-t_{n-i})$. Il r\'esulte du \S\
2.2 que
$$
\eqalign{
|\pi^*_{n,r}(\xi)|&=|\xi-t_{n+1}|\kern2pt|\pi_{n,r}(\xi)|\cr
&\leq (r+1)h_{\rm max}\cdot (r+1)!\,h^{r+1}_{\rm max}\leq (r+2)!\,h^{r+2}_{\rm max},\cr
\Big|\int^{t_{n+1}}_{t_n}(z'(t)-p^*_{n,r}(t))dt\Big|&\leq|z^{(r+3)}(\xi)|h_nh^{r+2}_{\rm max},\cr
}
$$
par cons\'equent l'erreur de consistance admet la majoration

\bigskip
\majdeb
\hfil$|e_n|\leq Ch_nh_{\rm max}^{r+2}(1+O(h_n))$,\hfill
\par\majfin\par

\bigskip
\noindent avec\quad $\displaystyle C=\max_{t\in[t_0,t_0+T]}|z^{(r+3)}|$.

\bigskip
\noindent{\it La m\'ethode $AM_{r+1}$ est donc d'ordre $r+2$.}
On initialisera les $r$ premi\`eres valeurs $y_1,\dots,y_r$ \`a l'aide d'une m\'ethode de
Runge-Kutta d'ordre $r+2$ (ou \`a la rigueur $r+1$).

\subsection 3.3.$^*$ Stabilit\'e de la m\'ethode $AM_{r+1}$|

\noindent On suppose que les rapports $h_n/h_{n-1}$ restent born\'es, de sorte que les quantit\'es
$$
\beta^*_r=\max_n\sum_{\leq i\leq r}|b^*_{n,i,r}|,\qquad \gamma^*_r=\max_n|b^*_{n,-1,r}|
$$
sont contr\^ol\'ees. Supposons \'egalement que $f(t,y)$ soit $k$-lipschitzienne en $y$. La m\'ethode
de r\'esolution it\'erative pour $y_{n+1}$ fonctionne d\`es que $h_n<\dfrac{1}{|b^*_{n,-1,r}|k}$, en
particulier d\`es que 
$$
h_{\rm max}<\frac{1}{\gamma^*_rk},
$$ 
ce que nous supposons d\'esormais. Soit $\widetilde y_n$ une suite
perturb\'ee telle que
$$
\cases{
\displaystyle\widetilde y_{n+1}=\widetilde y_n+h_n\Big(b^*_{n,-1,r}\widetilde f_{n+1}+\sum_{0\leq
i\leq r}b^*_{n,i,r}\widetilde f_{n-i}\Big)+\varepsilon_n\hfill\cr 
\noalign{\vskip3pt}
\widetilde f_{n-i}\kern1pt=f(t_{n-i},y_{n-i}),\qquad r\leq n<N,\hfill\cr
}
$$
et posons $\displaystyle\theta_n=\max_{0\leq i\leq n}|\widetilde y_i-y_i|$. Comme on a
$\theta_{n+1}=\max\,(|\widetilde y_{n+1}-y_{n+1}|,\theta_n)$, il vient 
$$
\eqalign{
\theta_{n+1}&\leq \theta_n+kh_n\Big(|b^*_{n,-1,r}|\theta_{n+1}+\sum_{0\leq i\leq
r}|b^*_{n,i,r}|\theta_n\Big)+|\varepsilon_n|,\cr 
\theta_{n+1}&(1-|b^*_{n,-1,r}|kh_n)\leq\theta_n\Big(1+\sum_{0\leq i\leq
r}|b^*_{n,i,r}|kh_n\Big)+|\varepsilon_n|\cr 
&\qquad\qquad\qquad\quad\kern6pt\leq\Big(1+\sum_{0\leq i\leq
r}|b^*_{n,i,r}|kh_n\Big)(\theta_n+|\varepsilon_n|).\cr 
}
$$
Or $1-|b^*_{n,-1,r}|kh_n\geq 1-\gamma^*_rkh_{\rm max}>0$, par suite
$$
\eqalign{
\theta_{n+1}&\leq\frac{\displaystyle 1+\sum_{0\leq i\leq r}|b^*_{n,i,r}|kh_n}{1-|b^*_{n,-1,r}|kh_n}\
(\theta_n+|\varepsilon_n|),\cr 
\noalign{\vskip4pt}
\theta_{n+1}&\leq\left(1+\frac{\displaystyle\sum_{-1\leq i\leq
r}|b^*_{n,i,r}|kh_n}{1-|b^*_{n,-1,r}|kh_n}\right)(\theta_n+|\varepsilon_n|),
\cr 
\noalign{\vskip4pt}
\theta_{n+1}&\leq (1+\Lambda h_n)(\theta_n+|\varepsilon_n|),\cr
}
$$
avec $\Lambda=\beta^*_rk/(1-\gamma^*_rkh_{\rm max})$. Un
raisonnement analogue \`a la d\'emonstration du lemme de Gronwall
discret donne par r\'ecurrence sur $n$~:
$$
\theta_n\leq e^{\Lambda(t_n-t_r)}\Big(\theta_r+\sum_{r\leq i\leq n}|\varepsilon_i|\Big),
$$
d'o\`u la constante de stabilit\'e

\bigskip
\majdeb
\hfil$S=e^{\Lambda T}=\exp\ \left(
\displaystyle\dfrac{\beta^*_rkT}{1-\gamma^*_rkh_{\rm max}}\right)$\hfill
\par\majfin\par

\bigskip
\noindent Lorsque $h_{\rm max}$ est assez petit devant $1/\gamma_r^*k$, on a donc sensiblement 
$$
S\simeq\exp\,(\beta^*_rkT).
$$
Le tableau du \S\ 3.1 montre qu'\`a ordre $r+2$ \'egal, la m\'ethode $AM_{r+1}$ est beaucoup plus
stable que $AB_{r+2}$. Il n'en reste pas moins que malgr\'e cet avantage important, la m\'ethode
d'Adams-Moulton est d'utilisation d\'elicate \`a cause de son caract\`ere implicite. Les m\'ethodes
de pr\'ediction-correction que nous allons d\'ecrire permettent d'obtenir une stabilit\'e
\'equivalente, tout en fournissant un sch\'ema explicite de r\'esolution.
\vfill\eject

\section 4. M\'ethodes de pr\'ediction-correction|

\subsection 4.1. Principe g\'en\'eral|

\noindent On se donne une m\'ethode dite de {\it pr\'ediction} (ou {\it
pr\'edicteur}), fournissant (explicitement) une premi\`ere valeur
approch\'ee $py_{n+1}$ du point $y_{n+1}$ \`a atteindre~:
$$
\eqalign{
py_{n+1}&=\hbox{pr\'ediction de $y_{n+1}$},\cr
\pf_{n+1}&=f(t_{n+1},py_{n+1})=\hbox{pr\'ediction de $f_{n+1}$}.\cr
}
$$
En substituant la valeur $\pf_{n+1}$ ainsi trouv\'ee \`a $f_{n+1}$ dans la formule d'Adams-Moulton,
on obtient alors une nouvelle valeur {\it corrig\'ee} $y_{n+1}$ qui est retenue en vue des calculs
ult\'erieurs.

De fa\c con pr\'ecise, une m\'ethode PECE (pr\'ediction, \'evaluation, 
correction, \'evaluation)\kern-0.54pt\break
\`a $r+1$ pas va s'\'ecrire de la mani\`ere suivante~: 
$y_{n-r},f_{n-r},\dots,y_n,f_n$ \'etant d\'ej\`a calcul\'es, on pose

$$
\cases{
\hbox{Pr\'ediction~:}\hfill&$py_{n+1}=\dots~~$\hbox{(\`a partir des 
$y_{n-i},f_{n-i},0\leq i\leq
r$)}\hfill\cr 
&$t_{n+1}=t_n+h_n$\hfill\cr
\noalign{\vskip3pt}
\hbox{Evaluation~:}\hfill&$\pf_{n+1}=f(t_{n+1},py_{n+1})$\hfill\cr
\noalign{\vskip3pt}
\displaystyle\hbox{Correction~:}\hfill&$y_{n+1}=y_n+h_n\Big(b^*_{n,-1,r}\pf_{n+1}+\sum_{0\leq i\leq
r}b^*_{n,i,r}f_{n-i}\Big)$\hfill\cr
\noalign{\vskip3pt}
\hbox{Evaluation~:}\hfill&$f_{n+1}=f(t_{n+1},y_{n+1})$\hfill\cr 
}
$$

\noindent Nous avons utilis\'e ici la m\'ethode d'Adams-Moulton comme {\it correcteur}, mais cela
pourrait \^etre {\it a priori} n'importe quelle autre m\'ethode implicite.

Ici encore, le d\'emarrage de l'algorithme PECE n\'ecessite le calcul pr\'ealable des points
$y_1,\dots,y_r$ et des pentes $f_0,\dots,f_r$ \`a l'aide d'une m\'ethode \`a 1 pas. On peut estimer
que le co\^ut en temps de calcul d'une m\'ethode PECE est environ le double de celui d'une
m\'ethode d'Adams-Bashforth d'ordre \'egal (mais on verra que la stabilit\'e est beaucoup
meilleure). Ce temps de calcul est en g\'en\'eral inf\'erieur \`a celui des m\'ethodes de
Runge-Kutta sophistiqu\'ees (d'ordre $\geq 3$).

\subsection 4.2. Erreur de consistance dans PECE|

\noindent Soit $z$ une solution exacte du probl\`eme de Cauchy. L'erreur de consistance est
$$
e_n=z(t_{n+1})-y_{n+1}\quad{\rm avec}\quad y_{n-i}=z(t_{n-i}),\quad 0\leq i\leq r.
$$
Soit $y^*_{n+1}$ la valeur qui serait obtenue \`a l'aide du seul correcteur (Adams-Moulton), de
sorte que
$$
\cases{
\displaystyle y^*_{n+1}=y_n+h_n\Big(b^*_{n,-1,r}f^*_{n+1}+\sum_{0\leq i\leq
r}b^*_{n,i,r}f_{n,i}\Big)\hfill\cr \noalign{\vskip3pt}
f^*_{n+1}=f(t_{n+1},y^*_{n+1}).\hfill\cr
}
$$
L'erreur de consistance correspondante est
$$
e^*_n=z(t_{n+1})-y^*_{n+1}.
$$
Le pr\'edicteur introduit lui aussi une erreur de consistance
$$
pe_n=z(t_{n+1})-py_{n+1}.
$$
\'Ecrivons
$$
\eqalign{
e_n&=(z(t_{n+1})-y^*_{n+1})+(y^*_{n+1}-y_{n+1})\cr
e_n&=e^*_n+(y^*_{n+1}-y_{n+1}).\cr
}
$$
On a par ailleurs
$$
y^*_{n+1}-y_{n+1}=h_nb^*_{n,-1,r}(f^*_{n+1}-\pf_{n+1}).
$$
Si $f(t,y)$ est $k$-lipschitzienne en $y$, on en d\'eduit
$$
|y^*_{n+1}-y_{n+1}|\leq h_n|b^*_{n,-1,r}|\,k\,|y^*_{n+1}-py_{n+1}|
$$
et $y^*_{n+1}-py_{n+1}=(z(t_{n+1})-py_{n+1})-(z(t_{n+1})-y^*_{n+1})$, d'o\`u
$$
|y^*_{n+1}-py_{n+1}|\leq |pe_n|+|e^*_n|.
$$
Il en r\'esulte finalement
$$
\eqalign{
&|y^*_{n+1}-y_{n+1}|\leq |b^*_{n,-1,r}|\,kh_n(|pe_n|+|e^*_n|)\cr
&|e_n|\leq |e^*_n|+|y^*_{n+1}-y_{n+1}|\cr
&|e_n|\leq \big(1+|b^*_{n,-1,r}|kh_n\big)|e^*_n|+|b^*_{n,-1,r}|\,kh_n\,|pe_n|.\cr
}
$$
On voit que l'influence du pr\'edicteur est nettement moindre que celle du correcteur puisque son
erreur de consistance est en facteur d'un terme $O(h_n)$. Le {\it correcteur $AM_{r+1}$ \'etant
d'ordre $r+2$} (c'est-\`a-dire $|e^*_n|\leq Ch_nh^{r+2}_{\rm max}$), on voit qu'il convient de {\it
choisir un pr\'edicteur d'ordre $r+1$.} Les contributions de $|e^*_n|$ et $|pe_n|$ dans $|e_n|$
seront alors toutes deux $\leq Ch_nh^{r+2}_{\rm max}$, {\it l'ordre global de PECE est donc $r+2$
dans ce cas.}

\subsection 4.3. Exemples|

\noindent{\tenpointbf $\bullet$ Pr\'edicteur~: Euler (ordre 1),~
Correcteur~: $AM_1$ (ordre 2).}
$$
\cases{
P~:\quad py_{n+1}=y_n+h_nf_n\hfill\cr
\noalign{\vskip3pt}
E~:\quad \pf_{n+1}=f(t_{n+1},py_{n+1})\hfill\cr
\noalign{\vskip3pt}
C~:\quad y_{n+1}=y_n+h_n\big(\dfrac{1}{2}\ \pf_{n+1}+\dfrac{1}{2}\ f_n
\big)\hfill\cr
\noalign{\vskip3pt}
E~:\quad f_{n+1}=f(t_{n+1},y_{n+1})\hfill\cr}
$$
Cet algorithme co\"\i ncide avec la m\'ethode de Heun, qui n'est  autre que la m\'ethode de
Runge-Kutta d\'efinie par  $$\matrix{
0&\big|&0&0\cr
\noalign{\vskip-2pt}
1&\big|&0&1\cr
\noalign{\vskip-8pt}
\noalign{\hbox to 5pc{\hrulefill}}
\noalign{\vskip-3pt}
&\Bigg|&\dfrac{1}{2}&\dfrac{1}{2}\cr
}.
$$

{\parindent=1pc\tenpointbf 
\item{$\bullet$} Pr\'edicteur~: Nystr\"om (ordre 2) 
avec pas constant $h_n=h$,\hfil\break Correcteur~: $AM_2$ (ordre 3).\par}

$$
\cases{
P~:\quad py_{n+1}=y_{n-1}+2hf_n\hfill\cr
\noalign{\vskip3pt}
E~:\quad \pf_{n+1}=f(t_{n+1},py_{n+1})\hfill\cr
\noalign{\vskip3pt}
C~:\quad
y_{n+1}=y_n+h\Big(\dfrac{5}{12}\,\pf_{n+1}+\dfrac{8}{12}\,f_n-\dfrac{1}{12}\,f_{n-1}\Big)\hfill\cr 
\noalign{\vskip3pt}
E~:\quad f_{n+1}=f(t_{n+1},y_{n+1})\hfill\cr
}
$$
\medskip

\noindent{\tenpointbf $\bullet$ Pr\'edicteur~: $AB_{r+1}$ (ordre $r+1$),
Correcteur~: $AM_{r+1}$ (ordre $r+2$).}
$$
\cases{
\displaystyle P~:\quad py_{n+1}=y_n+h_n\sum_{0\leq i\leq r}b_{n,i,r}f_{n-i}\hfill\cr
\noalign{\vskip3pt}
E~:\quad \pf_{n+1}=f(t_{n+1},py_{n+1})\hfill\cr
\noalign{\vskip3pt}
\displaystyle C~:\quad y_{n+1}=y_n+h_n\Big(b^*_{n,-1,r}\pf_{n+1}+\sum_{0\leq i\leq
r}b^*_{n,i,r}f_{n-i}\Big)\hfill\cr
\noalign{\vskip3pt} 
E~:\quad f_{n+1}=f(t_{n+1},y_{n+1}).\hfill\cr
}
$$


\bigskip
\majedeb
\state Exercice -- |  {\it V\'erifier que ce dernier algorithme PECE 
\'equivaut \`a la m\'ethode d'Adams-Moulton dans laquelle
l'algorithme it\'eratif est arr\^et\'e \`a la premi\`ere \'etape,
soit $y_{n+1}=x_1$, calcul\'e \`a partir de la valeur $x_0$ fourni
par la m\'ethode d'Adams-Bashforth.}
\par\majefin\par


\subsection 4.4. Stabilit\'e de la m\'ethode PECE|

\noindent Supposons que le pr\'edicteur soit de la forme
$$
py_n=\sum_{0\leq i\leq r}\alpha_{n,i}y_{n-i}+h_n\sum_{0\leq i\leq r}
\beta_{n,i}f_{n-i}
$$
et notons
$$
A=\max_n\sum_i|\alpha_{n,i}|,\qquad B=\max_n\sum_i|\beta_{n,i}|.
$$
Soit $\widetilde y_n$ une suite perturb\'ee telle que
$$
\cases{
\displaystyle p\widetilde y_{n+1}=\sum_{0\leq i\leq r}
\alpha_{n,i}\widetilde y_{n-i}+h_n\sum_{0\leq i\leq r}
\beta_{n,i}\widetilde f_{n-i}\hfill\cr
\noalign{\vskip3pt}
\displaystyle \widetilde y_{n+1}=\widetilde y_n+h_n\Big(b^*_{n,-1,r}p
\widetilde f_{n+1}+\sum_{0\leq i\leq r}
b^*_{n,i,r}\widetilde f_{n-i}\Big)+\varepsilon_n.\hfill\cr}
$$

\smallskip
\state Remarque -- |  Dans la r\'ealit\'e il s'introduit \'egalement une
erreur d'arrondi $p\varepsilon_n$ au niveau de la pr\'ediction, mais pour
notre calcul il sera plus simple de comptabiliser cette erreur dans
$\varepsilon_n$ (ceci n'est visiblement pas restrictif).
\smallskip

\noindent Posons $\displaystyle\theta_n=\max_{0\leq i\leq n}|\widetilde
y_i-y_i|$ et $p\theta_n=|p\widetilde y_n-py_n|$. Comme $f(t,y)$ est
suppos\'ee $k$-lipschitzienne en $y$, il vient~:
$$
\cases{
p\theta_{n+1}\leq A\theta_n+h_nBk\theta_n\hfill\cr
\noalign{\vskip3pt}
\displaystyle\theta_{n+1}\leq \theta_n+kh_n\Big(|b^*_{n,-1,r}|p\theta_{n+1}
+\sum_{0\leq i\leq r}|b^*_{n,i,r}|\theta_n\Big)+|\varepsilon_n|.\hfill\cr}
$$
En substituant $p\theta_{n+1}\leq \theta_n+\theta_n(A-1+Bkh_{n})$ dans la
deuxi\`eme ligne il vient
$$
\eqalign{
\theta_{n+1}&\leq \theta_n\left(1+\Big[\sum_{-1\leq i\leq r}|b^*_{n,i,r}|
+|b^*_{n,-1,r}|(A-1+Bkh_n)\Big]kh_n\right)+|\varepsilon_n|,\cr
\theta_{n+1}&\leq \theta_n(1+\Lambda h_n)+|\varepsilon_n|,\cr}
$$
avec $\Lambda=(\beta^*_r+\gamma^*_r(A-1+Bkh_{\rm max}))k$. Le lemme de
Gronwall montre donc que la m\'ethode PECE est stable, avec constante
de stabilit\'e

\bigskip
\majdeb
\hfil$S=e^{\Lambda T}=\exp\Big((\beta^*_r+\gamma^*_r
(A-1+Bkh_{\rm max})kT\Big)$.\hfill
\par\majfin\par

\bigskip
\noindent Si $h_{\rm max}$ est petit, on va avoir
$$
S\simeq \exp\Big((\beta^*_r+\gamma^*_r(A-1))kT\Big).
$$
On voit que la stabilit\'e du pr\'edicteur n'a pas d'incidence sur
la stabilit\'e de la m\'ethode PECE, seule la valeur de la constante
$A$ peut influer sur cette stabilit\'e~; on pourrait en th\'eorie
utiliser un pr\'edicteur instable~! La consistance du pr\'edicteur
implique $\displaystyle\sum_i\alpha_{n,i}=1$. Si les coefficients
$\alpha_{n,i}$ sont $\geq 0$, alors on a $A=1$ (c'est le cas des
m\'ethodes de Nystr\"om, Milne, ou $AB_{r+1}$), par cons\'equent la
constante de stabilit\'e $$
S\simeq\exp\,(\beta^*_rkT)
$$
sera peu diff\'erente de celle de la m\'ethode d'Adams-Moulton seule. On obtient donc des
m\'ethodes assez stables, d'un co\^ut mod\'er\'e en temps de calcul et d'ordre aussi \'elev\'e que
l'on veut. Par comparaison avec les m\'ethodes de Runge-Kutta, elles sont un peu plus rapides mais
un peu moins stables \`a ordre \'egal.

\subsection 4.5. M\'ethodes PEC|

\noindent Comme leur nom l'indique, il s'agit de m\'ethodes de
pr\'ediction-correction dans lesquelles la derni\`ere \'etape
d'\'evaluation est omise (en vue bien s\^ur de gagner du temps).  Ceci
signifie que les pentes corrig\'ees $f_{n+1}$ ne sont pas calcul\'ees,
il faudra donc se contenter de faire intervenir les pentes pr\'edites
$\pf_{n-i}$. On obtient alors l'algorithme suivant~:
$$
\cases{
\hbox{Pr\'ediction}~:\hfill&$\displaystyle py_{n+1}=\sum_{0\leq i\leq
r}\alpha_{n,i}y_{n-i}+h_n\sum_{0\leq i\leq r}\beta_{n,i}\pf_{n-i}$\hfill\cr 
&$t_{n+1}=t_n+h_n$\hfill\cr
\noalign{\vskip3pt}
\hbox{Evaluation}~:\hfill&$\pf_{n+1}=f(t_{n+1},py_{n+1})$\hfill\cr
\noalign{\vskip3pt}
\hbox{Correction}~:\hfill&$\displaystyle y_{n+1}=y_n+h_n\sum_{-1\leq i\leq
r}b^*_{n,i,r}\pf_{n-i}$.\hfill\cr }
$$
Le d\'emarrage de l'algorithme n\'ecessite le calcul pr\'ealable des quantit\'es $y_1,\dots,y_r$,
$\pf_0,\dots,\pf_r$.

\medskip
\state Erreur de consistance* -- | L'algorithme PEC n'entre pas tout 
\`a fait dans le cadre g\'en\'eral des m\'ethodes que nous avons
consid\'er\'ees jusqu'\`a pr\'esent. Il convient de red\'efinir $e_n$
comme suit. Si $z$ est une solution exacte, on pose
$$
e_n=z(t_{n+1})-y_{n+1}
$$
o\`u $y_{n+1}$ est calcul\'ee \`a partir des valeurs ant\'erieures
$py_{n-i}=y_{n-i}=z(t_{n-i})$, $0\leq i\leq r$. Avec cette
d\'efinition, il est facile de voir que l'erreur de consistance est
identique \`a celle de la m\'ethode PECE, d'o\`u
$$
|e_n|\leq \Big(1+|b^*_{n,-1,r}|kh_n\Big)|e^*_n|+|b^*_{n,-1,r}|kh_n|pe_n|.
$$
Le correcteur \'etant d'ordre $r+2$, on choisira ici encore le pr\'edicteur
d'ordre $r+1$.

\bigskip
\state Stabilit\'e de la m\'ethode PEC* -- | Avec les notations et
hypoth\`eses du \S~4.4, consid\'erons une suite perturb\'ee $\widetilde y_n$
telle que
$$
\widetilde y_{n+1}=\widetilde y_n+h_n\sum_{-1\leq i\leq r}b^*_{n,i,r}p\widetilde
f_{n-i}+\varepsilon_n 
$$
et posons $\displaystyle\theta_n=\max_{0\leq i\leq n}|\widetilde y_i-y_i|$,\ \ 
$\displaystyle p\theta_n=\max_{0\leq i\leq n}|p\widetilde y_i-py_i|$. Il vient~:
$$
\cases{
p\theta_{n+1}\leq A\theta_n+Bkh_np\theta_n\hfill\cr
\noalign{\vskip3pt}
\theta_{n+1}\leq\theta_n+\beta^*_rkh_np\theta_{n+1}+|\varepsilon_n|.\hfill\cr
}
$$
La premi\`ere ligne entra\^\i ne
$$
p\theta_{n+1}\leq A\theta_n+Bkh_{\rm max}p\theta_{n+1},
$$
d'o\`u $p\theta_{n+1}\leq\dfrac{A}{1-Bkh_{\rm max}}\ \theta_n$ si $Bkh_{\rm max}<1$. En substituant
dans la deuxi\`eme ligne il vient~: $$
\theta_{n+1}\leq\left(1+\frac{\beta^*_rAk}{1-Bkh_{\rm max}}\ h_n\right)\theta_n+|\varepsilon_n|.
$$
Le lemme de Gronwall donne la constante de stabilit\'e

\bigskip
\majdeb
\hfil$S=\exp\,\left(\displaystyle
\dfrac{\beta^*_rAkT}{1-Bkh_{\rm max}}\right)$.\hfill
\par\majfin\par

\bigskip
\noindent Si $h_{\rm max}$ est assez petit, on aura
$$
S\simeq\exp\,(\beta_{r}^*AkT).
$$
Ceci est un peu moins bon que dans le cas de la m\'ethode PECE, car $\gamma^*_r<\beta^*_r$.
N\'eanmoins, pour $A=1$ la constante de stabilit\'e est la m\^eme~:
$$
S\simeq \exp\,(\beta^*_r kT),
$$
c'est-\`a-dire pr\'ecis\'ement la constante de stabilit\'e de la
m\'ethode d'Adams-Moulton seule.  Par rapport \`a PECE, on \'economise
donc un peu de temps de calcul, mais on perd un peu en stabilit\'e et
en pr\'ecision.

\section 5. Probl\`emes|

{\parindent1.4pc
  
\noindent{\bf 5.1.} On consid\`ere le probl\`eme de Cauchy
$y'(t)=f(t,y(t))$, $y(t_0)=y_0$, o\`u
$f:[t_0,t_0+T]\times\bR\rightarrow\bR$ est une fonction de classe
$C^5$.  Pour r\'esoudre num\'eriquement ce probl\`eme, on se donne un
entier $N\geq 2$ et on consid\`ere la subdivision $t_n=t_0+nh$, $0\leq
n\leq N$, de pas constant $h=\dfrac{T}{N}$.

\noindent On \'etudie les m\'ethodes \`a 2 pas de la forme
$$
y_{n+1}=\alpha y_{n-1}+\alpha'y_n+h(\beta f_{n-1}+\beta'f_n+\beta''f_{n+1}),
\leqno{\rm (M)}
$$
avec $f_n=f(t_n,y_n)$. La m\'ethode est donc explicite si
$\beta''=0$ et implicite si $\beta''\not=0$.  \medskip

\item{(a)} Soit $g$ une fonction de classe $C^5$ au voisinage de 0.

\item{} Calculer un d\'eveloppement limit\'e \`a l'ordre 4 en $h=0$ de la 
quantit\'e
$$
\Delta(h)=g(h)-[\alpha g(-h)+\alpha'g(0)+h(\beta g'(-h)+\beta'g'(0)+\beta'' g'(h))].
$$
\item{} \'Ecrire la condition n\'ecessaire et suffisante pour que la m\'ethode (M) soit d'ordre
$\geq 1$ (respectivement $\geq 2$, $\geq 3$, $\geq 4$). Montrer que la seule m\'ethode (M) qui soit
d'ordre $\geq 4$ est 
$$
y_{n+1}=y_{n-1}+h\Big(\dfrac{1}{3}\ f_{n-1}+\dfrac{4}{3}\ f_n+\dfrac{1}{3}\
f_{n+1}\Big).\leqno({\rm M}^4) 
$$
\item{} Quelle interpr\'etation peut-on donner de cette m\'ethode~?

\medskip
\noindent NB~: Dans les questions qui suivent, on supposera sans le repr\'eciser chaque fois que les
m\'ethodes (M) \'etudi\'ees sont d'ordre 1 au moins.

\medskip
\item{(b)} On cherche \`a tester la stabilit\'e de la m\'ethode (M) en consid\'erant l'\'equation
diff\'erentielle triviale $y'=0$.

\item{} Les r\'eels $y_0$ et $y_1=y_0+\varepsilon$ \'etant donn\'es {\it a priori}, exprimer $y_n$
en fonction de $y_0$, $\varepsilon$, $\alpha$, $n$. En d\'eduire qu'une condition n\'ecessaire pour
que la m\'ethode (M) soit stable est que $-1<\alpha\leq 1$.

\medskip
\item{(c)} On se propose ici de montrer inversement que la m\'ethode (M) est stable, si $0\leq
\alpha\leq 1$ et si $h$ est assez petit. On suppose que pour tout $t\in[t_0,t_0+T]$ la fonction
$y\mapsto f(t,y)$ est $k$-lipschitzienne. Soient deux suite $(y_n)$ et $(z_n)$ telles que pour
$n\geq 1$ on ait
$$
\eqalign{
y_{n+1}&=\alpha y_{n-1}+\alpha'y_n+h(\beta
f(t_{n-1},y_{n-1})+\beta'f(t_n,y_n)+\beta''f(t_{n+1},y_{n+1})),\cr 
z_{n+1}&=\alpha z_{n-1}+\alpha' z_n+h(\beta
f(t_{n-1},z_{n-1})+\beta'f(t_n,z_n)+\beta''f(t_{n+1},z_{n+1}))
+\varepsilon_n.\kern-8pt\cr 
}
$$
On pose
$$
\theta_n=\max_{0\leq i\leq n}|z_i-y_i|.
$$

\itemitem{$(\alpha)$} Majorer $|z_{n+1}-y_{n+1}|$ en fonction de 
$\theta_n$, $\theta_{n+1}$ et $\varepsilon_n$.

\itemitem{} En d\'eduire que si $|\beta''|kh<1$ on a
$$
\theta_{n+1}\leq\left(1+\frac{(|\beta|+|\beta'|+|\beta''|)kh}{1-
|\beta''|kh}\right)\ (\theta_n+|\varepsilon_n|).
$$

\itemitem{$(\beta)$} En d\'eduire l'existence d'une constante de 
stabilit\'e $S(h)$, qui reste born\'ee quand $h$ tend vers 0, telle que
$$
\theta_N\leq S(h)\Big(\theta_1+\sum_{k=1}^{N-1}|\varepsilon_n|\Big).
$$

\item{(e)} D\'eterminer en fonction de $\alpha$ les m\'ethodes (M) qui
sont d'ordre 3~; on notera celles-ci (M$^3_{\alpha})$. A quoi
correspond (M$^3_0)$~? Existe-t-il une m\'ethode (M$^3_{\alpha})$ qui
soit explicite et stable~?

\item{} Montrer qu'il existe une unique m\'ethode (M$^3_{\alpha_1})$
pour laquelle $\beta=0$.

\medskip
\item{(f)}\vskip-\parskip\vskip-\baselineskip
\itemitem{$(\alpha)$} Expliciter l'algorithme PECE dont le pr\'edicteur 
est la m\'ethode de Nystr\"om et dont le correcteur est la m\'ethode 
(M$^3_{\alpha_1}$).

\itemitem{} L'initialisation sera faite au moyen de la m\'ethode de
Runge-Kutta d'ordre 4 usuelle.

\smallskip
\itemitem{$(\beta)$} \'Ecrire un programme Pascal mettant en {\oe}uvre 
l'algorithme pr\'ec\'edent dans le cas de la fonction $f(t,y)=\sin(ty-y^2)$.

\itemitem{} L'utilisateur fournit la donn\'ee initiale $(t_0,y_0)$, le 
pas $h$, et le nombre $N$ d'it\'erations. L'ordinateur affichera alors 
les valeurs $(t_n,y_n)$ successives pour $0\leq n\leq N$.

\bigskip
\noindent{\bf 5.2.} L'objet de ce probl\`eme est d'\'etudier les m\'ethodes 
d'Adams-Bashforth et d'Adams-Moulton avec pas constant.
\medskip

\item{(a)} Montrer que la m\'ethode d'Adams-Bashforth \`a $(r+1)$ pas, de 
pas constant $h$, s'\'ecrit
$$
y_{n+1}=y_n+h\sum^r_{i=0}b_{i,r}f(t_{n-i},y_{n-i})
$$
avec
$$
b_{i,r}=(-1)^i\int^1_0\frac{s(s+1)\dots\widehat {(s+i)}\dots(s+r)}{i!(r-i)!},\quad 0\leq i\leq r.
$$

\item{(b)} On pose
$$
\gamma_r=\int^1_0\frac{s(s+1)\dots(s+r-1)}{r!}\ ds.
$$
\item{} D\'emontrer les formules
$$
\eqalign{
&b_{i,r}-b_{i,r-1}=(-1)^iC^i_{r}\gamma_r,\quad 0\leq i\leq r-1,\cr
&b_{r,r}=(-1)^r\gamma_r.\cr
}
$$

\item{(c)} Montrer que pour $|t|<1$ on a
$$
\int^1_0(1-t)^{-s}ds=\sum^{+\infty}_{r=0}\gamma_r t^r.
$$
\item{} En d\'eduire la valeur de l'expression
$\displaystyle\log\,(1-t)\sum^{+\infty}_{r=0}\gamma_rt^r$, puis la valeur des sommes
$$
\dfrac{\gamma_0}{r+1}+\frac{\gamma_1}{r}+\dots+\frac{\gamma_{r-1}}{2}+\gamma_r.
$$

\item{(d)} \'Ecrire un programme Pascal mettant en {\oe}uvre la m\'ethode d'Adams-Bashforth \`a un
nombre arbitraire de pas. On utilisera  les formules de r\'ecurrence ci-dessus pour \'evaluer
$\gamma_r$ et $b_{i,r}$. L'initialisation sera faite au moyen de la m\'ethode de Runge-Kutta
d'ordre 4.

\medskip
\item{(e)} D\'emontrer des formules analogues \`a celles de $(a)$, $(b)$, $(c)$ pour la m\'ethode
d'Adams-Moulton.

\bigskip
\noindent{\bf 5.3.} Soit \`a r\'esoudre num\'eriquement un probl\`eme de Cauchy
$$
y'=f(t,y),\quad y(t_0)=y_0
$$
o\`u $f$ est de classe $C^2$ sur $[t_0,t_0+T]\times\bR$. On se donne une subdivision
$$
t_0<t_1<\dots<t_N=t_0+T
$$
de l'intervalle, et on \'etudie la m\'ethode num\'erique suivante~: si $y_n$ est la valeur
approch\'ee de la solution au temps $t_n$ et si $f_n=f(t_n,y_n)$ on pose
$$
y_{n+1}=y_{n-1}+\int^{t_{n+1}}_{t_{n-1}}p_n(t)dt\leqno{\rm (M)}
$$
o\`u $p_n$ est le polyn\^ome d'interpolation des pentes $f_n$, $f_{n-1}$ aux temps $t_n$ et
$t_{n-1}$. On note $h_n=t_{n+1}-t_n$.
\medskip

\item{(a)} Calculer explicitement $y_{n+1}$. Quelle m\'ethode obtient-on lorsque le pas $h_n=h$
est constant~?

\medskip
\item{(b)} Soit $z$ une solution exacte de l'\'equation diff\'erentielle.

\item{} D\'eterminer un \'equivalent de l'erreur de consistance $e_n$ attach\'ee \`a la m\'ethode
(M). Quel est l'ordre de cette m\'ethode~? Comment proc\'ederiez-vous pour la phase
d'initialisation~?

\medskip
\item{(c)} Soit une suite $\widetilde y_n$ v\'erifiant la relation de r\'ecurrence
$$
\widetilde y_{n+1}=\widetilde y_{n-1}+\int^{t_{n+1}}_{t_{n-1}}\widetilde
p_n(t)dt+\varepsilon_n
$$
\item{} o\`u $\widetilde p_n$ interpole les valeurs $\widetilde f_n=f(t_n,\widetilde y_n)$ et
$\widetilde y_{n-1}$. La quantit\'e $\varepsilon_n$ d\'esigne l'erreur commise \`a l'\'etape $n$.
On suppose que $f(t,y)$ est $k$-lipschizienne en $y$ et on note
$\displaystyle\theta_n=\max_{0\leq i\leq n}|\widetilde y_i-y_i|$.

\smallskip
\itemitem{$(\alpha)$} Montrer que
$$
\theta_{n+1}\leq \left(1+kh_n\left(1+{\rm max}\
\left\{\dfrac{h_n}{h_{n-1}},\dfrac{h_{n-1}}{h_n}\right\}\right)\right)
\theta_n+\varepsilon_n. 
$$

\itemitem{$(\beta)$} On suppose que le rapport de 2 pas cons\'ecutifs 
est major\'e par une constante
$\delta$ (avec disons $1\leq \delta\leq 2$). \'Etudier la stabilit\'e de 
la m\'ethode (M).

\bigskip
\noindent{\bf 5.4.} Dans cet exercice, on se propose d'\'etudier une 
m\'ethode de pr\'ediction-correction de type PEPEC pour la r\'esolution 
d'un probl\`eme de Cauchy
$$
\cases{
y'\quad\kern3pt=f(t,y),\qquad t\in[t_0,t_0+T]\cr
\noalign{\vskip3pt}
y(t_0)=y_0.\hfill\cr
}
$$
La fonction $f(t,y)$ est suppos\'ee de classe $C^4$ sur $[t_0,t_0+T]\times\bR$ et lipschitzienne de
rapport $k$ en $y$. Le pas est choisi constant~: $h=\dfrac{T}{N}$, $N\in\bN^*$.
\medskip

\item{(a)} L'objet de cette question est de d\'ecrire la m\'ethode de correction.

\smallskip
\itemitem{$(\alpha)$} Si $z$ est une solution exacte du probl\`eme de 
Cauchy, on \'ecrit
$$
z(t_{n+1})=z(t_n)+\int^{t_{n+1}}_{t_n}f(t,z(t))dt
$$
\itemitem{} et on approxime l'int\'egrale par la m\'ethode de Simpson 
\'el\'ementaire. Montrer que
l'algorithme correspondant s'\'ecrit
$$
y_{n+1}=y_n+h(\alpha f_n+\beta f_{n+\frac{1}{2}}+\gamma f_{n+1})\leqno{\rm (C)}
$$
\itemitem{} avec des coefficients $\alpha$, $\beta$, $\gamma$ que l'on 
pr\'ecisera.

\smallskip
\itemitem{$(\beta)$} On suppose que les pentes $f_n$, $f_{n+\frac{1}{2}}$, 
$f_{n+1}$ sont les d\'eriv\'ees exactes de $z$ aux points $t_n$, 
$t_n+\dfrac{h}{2}$, $t_n+h$. D\'eterminer un
\'equivalent de l'erreur de consistance 
$$
e^{*}_n=z(t_{n+1})-y_{n+1}
$$
en fonction de $h$ et de la d\'eriv\'ee $z^{(5)}(t_n)$. Quel est l'ordre de 
la m\'ethode~(C)~?

\medskip
\item{(b)} Pour pouvoir exploiter la formule (C), il est n\'ecessaire de pr\'edire des valeurs
approch\'ees $py_{n+\frac{1}{2}}$ et $py_{n+1}$ des points $y_{n+\frac{1}{2}}$ et $y_{n+1}$, ainsi
que les pentes correspondantes
$$
\pf_{n+\frac{1}{2}}=f(t_n+\dfrac{h}{2},py_{n+\frac{1}{2}}), \quad \pf_{n+1}=f(t_{n+1},py_{n+1}).
$$
\item{} Pour cela, on utilise comme pr\'edicteur (P) la m\'ethode
d'Adams-Bashforth \`a 3 pas, de pas constant $\dfrac{h}{2}$, faisant
intervenir les pentes pr\'edites ant\'erieures $\pf_n$,
$\pf_{n-\frac{1}{2}}$, $\pf_{n-1}$. La m\'ethode (P) est utilis\'ee
une premi\`ere fois pour \'evaluer $py_{n+\frac{1}{2}}$, puis une
deuxi\`eme fois pour \'evaluer $py_{n+1}$ \`a patir de
$py_{n+\frac{1}{2}}$.

\smallskip
\itemitem{$(\alpha)$} Expliciter compl\`etement l'algorithme PEPEC ainsi 
obtenu.

\smallskip
\itemitem{$(\beta)$} \'Ecrire en Pascal la boucle centrale correspondant 
\`a l'it\'eration de l'algori\-thme PEPEC (on pr\'ecisera la signification 
des variables utilis\'ees).

\smallskip 
\itemitem{$(\gamma)$} L'erreur de consistance de la m\'ethode PEPEC
est d\'efinie par $e_n=z(t_{n+1})-y_{n+1}$ o\`u l'on suppose que les
points $y_n$ et $py_{n-i}$, $i\in\big\{0,\dfrac{1}{2},1\big\}$ sont
exacts. Si $pe_n$ et $pe_{n+\frac{1}{2}}$ d\'esignent les erreurs de
consistance du pr\'edicteur sur les intervalles de temps
$\big[t_n,t_n+\dfrac{h}{2}\big]$ et
$\big[t_n+\dfrac{h}{2},t_{n+1}\big]$, montrer que
$$
|e_n|\leq|e^*_n|+\frac{4}{6}\ kh|pe_n|+\frac{1}{6}\
kh\Big(|pe_{n+\frac{1}{2}}|+|pe_n|+\frac{23}{24}\ kh|pe_n|\Big). 
$$
Le choix du pr\'edicteur est-il justifi\'e~? Quel est l'ordre de la
m\'ethode PEPEC~?  Comment proc\'ederiez-vous pour l'initialisation de
l'algorithme~?

\medskip
\item{(c)} On se propose ici de d\'eterminer une constante de stabilit\'e de l'algorithme PEPEC.
Soit $\widetilde y_n$ une suite perturb\'ee telle que la formule de correction soit entach\'ee
d'une erreur $\varepsilon_n$, toute l'erreur \'etant comptabilis\'ee au niveau de (C)~:
$$
\widetilde y_{n+1}=\widetilde y_n+h(\alpha\, p\widetilde f_n+\beta\, p\widetilde
f_{n+\frac{1}{2}}+\gamma\, p\widetilde f_{n+1})+\varepsilon_n. 
$$
\item{} Pour $i,n\in\bN$ on pose
$$
\eqalign{
\theta_n&=\max_{0\leq i\leq n}|\widetilde y_i-y_i|,\cr
p\theta_n&=\max_{0\leq i\leq n}|p\widetilde y_i-py_i|,\quad
p\theta_{n+\frac{1}{2}}=\max_{0\leq i\leq n}|p\widetilde
y_{i+\frac{1}{2}}-py_{i+\frac{1}{2}}|.\cr  
}
$$

\itemitem{$(\alpha)$} Majorer $p\theta_{n+\frac{1}{2}}$ en fonction de 
$\theta_n$, $p\theta_n$ et $p\theta_{n-\frac{1}{2}}$ puis
$p\theta_{n+1}$ en fonction de $p\theta_{n+\frac{1}{2}}$ et
$p\theta_n$. En d\'eduire successivement (pour $h$ assez petit)~:
$$
\eqalign{
p\theta_{n+1}&\leq\frac{1+\dfrac{7}{6}\ kh}{1-\dfrac{4}{6}\ kh}\ p\theta_{n+\frac{1}{2}}\leq
\frac{1}{1-\dfrac{11}{6}\ kh}\ p\theta_{n+\frac{1}{2}},\cr
p\theta_{n+\frac{1}{2}}&\leq\theta_n+kh\ \frac{\dfrac{11}{6}\
p\theta_{n-\frac{1}{2}}}{1-\dfrac{11}{6}\ kh},\cr
p\theta_{n+\frac{1}{2}}&\leq \frac{1-\dfrac{11}{6}\ kh}{1-\dfrac{11}{3}\
kh}\ \theta_n.\cr 
}
$$
En d\'eduire une majoration de $\theta_{n+1}$ en fonction de $\theta_n$,
$|\varepsilon_n|$ et une estimation de~$S$.

\medskip
\itemitem{$(\beta)$} Comparer la stabilit\'e de PEPEC \`a la stabilit\'e de
la m\'ethode PECE de m\^eme ordre ayant pour pr\'edicteur
Adams-Bashforth et pour correcteur Adams-Moulton. 
\par}

\bye
