% 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 VIII

\input equa-mac.tex

\blankline\blankline
\hfill\chtitle Chapitre VIII|%
\blankline\blankline
\hfill\title M\'ethodes num\'eriques \`a un pas|
\chaptitlerunning VIII -- M\'ethodes num\'eriques \`a un pas|
\vskip4pc

\noindent L'objectif de ce chapitre est de d\'ecrire un certain nombre de 
m\'ethodes permettant de r\'esoudre num\'eriquement le probl\`eme de Cauchy 
de condition initiale $y(t_0)=y_0$ pour une \'equation diff\'erentielle
$$
y'=f(t,y),\leqno{\rm (E)}
$$
o\`u $f:[t_0,t_0+T]\times\bR\rightarrow\bR$ est une fonction
suffisamment r\'eguli\`ere. Nous avons choisi ici d'exposer le cas des
\'equations unidimensionnelles dans le seul but de simplifier les
notations~; le cas des syst\`emes dans $\bR^m$ est tout \`a fait
identique, \`a condition de consid\'erer $y$ comme une variable
vectorielle et $f$ comme une fonction vectorielle dans les algorithmes
qui vont \^etre d\'ecrits.

\'Etant donn\'e une subdivision $t_0<t_1<\dots<t_N=t_0+T$ de
$[t_0,t_0+T]$, on cherche \`a d\'eterminer des valeurs approch\'ees
$y_0,y_1,\dots,y_N$ des valeurs $y(t_n)$ prises par la solution exacte
$y$. On notera les pas successifs
$$
h_n=t_{n+1}-t_n,\qquad 0\leq n\leq N-1,
$$
et
$$
h_{\max}=\max\,(h_n)
$$
le maximum du pas.

\medskip
On appelle {\it m\'ethode \`a un pas} une m\'ethode permettant de
calculer $y_{n+1}$ \`a partir de la seule valeur ant\'erieure $y_n$.
Une m\'ethode \`a $r$ pas est au contraire une m\'ethode qui
utilise les $r$ valeurs ant\'erieures $y_n,\,\ldots\,,y_{n-r+1}$
(valeurs qui doivent donc \^etre m\'emoris\'ees) afin de faire le
calcul de~$y_{n+1}$.

\section 1. D\'efinition des m\'ethodes \`a un pas, exemples|

\subsection 1.1. D\'efinitions|

\noindent Les m\'ethodes \`a un pas sont les m\'ethodes de r\'esolution 
num\'erique qui peuvent s'\'ecrire sous la forme
$$
y_{n+1}=y_n+h_n\Phi(t_n,y_n,h_n),\qquad 0\leq n<N,
$$
o\`u $\Phi:[t_0,t_0+T]\times\bR\times\bR\rightarrow\bR$ est une
fonction que l'on supposera continue.  Dans la pratique, la fonction
$\Phi(t,y,h)$ peut n'\^etre d\'efinie que sur une partie de la forme
$[t_0,t_0+T]\times J\times[0,\delta]$ o\`u $J$ est un intervalle de
$\bR$ (de sorte en particulier que $[t_0,t_0+T]\times J$ soit contenu
dans le domaine de d\'efinition de l'\'equation diff\'erentielle).
\bigskip

\state Exemple -- | La m\'ethode d'Euler est la m\'ethode \`a un pas 
associ\'ee \`a la fonction $\Phi(t,y,h)=f(t,y)$, et
d\'efinie par la formule de r\'ecurrence $y_{n+1}=y_n+h_nf(t_n,y_n)$
(voir chapitre V, \S$\,$2.3).

\bigskip
\majdeb
\state D\'efinition -- |  {\it L'erreur de consistance $e_n$ relative \`a 
une solution exacte $z$ est l'erreur
$$
e_n=z(t_{n+1})-y_{n+1},\quad 0\leq n<N
$$
produite par application de l'algorithme $y_{n+1}=y_n+h_n\Phi(t_n,y_n,h_n)$
\`a partir de la valeur $y_n=z(t_n)$. Autrement dit, 
cette erreur mesure l'\'ecart entre la valeur exacte $z(t_{n+1})$ au
temps $t_{n+1}$, et la valeur approch\'ee $y_{n+1}$ issue de la valeur
$y_n=z(t_n)$ prise comme valeur initiale au temps~$t_n$ $($une seule
\'etape de l'algorithme est donc mise en jeu$)$. En termes de la 
fonction~$\Phi$, on~a
$$
e_n=z(t_{n+1})-z(t_n)-h_n\Phi(t_n,z(t_n),h_n).
$$
}
\vskip-\baselineskip 
\majfin
\vskip\baselineskip 

\InsertFig 25.000 58.000 {
1 mm unit
  0.000  10.000 moveto 
 80.000   0.000 vector 
 10.000   5.000 moveto 
 45.000  90.000 vector 
 25.000  18.000 moveto
  0.340 setlinewidth
[  25.000  18.000   35.000  20.000   45.000  24.000   55.000  29.000  
 65.000  36.000   70.000  41.000  ] curve 
stroke 
 10.000  32.000 moveto 
  0.700   0.000 indent 
 10.000  29.000 moveto 
  0.700   0.000 indent 
 10.000  20.000 moveto 
  0.700   0.000 indent 
 35.000  10.000 moveto 
  0.700  90.000 indent 
 60.000  10.000 moveto 
  0.700  90.000 indent 
  0.113 setlinewidth
 35.000  20.000 moveto 
  0.700 disk 
 60.000  29.000 lineto 
stroke 
 60.000  29.000 moveto 
  0.700 disk 
 60.000  32.000 moveto 
  0.700 disk 
 69.000  29.000 moveto 
  0.000  90.000 vector 
 69.000  32.000 moveto 
  0.000 -90.000 vector 
 69.000  29.000 lineto 
stroke 
[ 0.300   0.600 ] 0 setdash 
 11.500  29.000 moveto 
 69.000  29.000 lineto 
stroke 
 11.500  32.000 moveto 
 69.000  32.000 lineto 
stroke 
 11.500  20.000 moveto 
 35.000  20.000 lineto 
stroke 
 35.000  11.500 moveto 
 35.000  20.000 lineto 
stroke 
 60.000  11.500 moveto 
 60.000  32.000 lineto 
stroke}
\LabelTeX   67.000  41.000 $z$\ELTX
\LabelTeX   71.000  29.700 $e_n$\ELTX
\LabelTeX   34.000   6.200 $t_n$\ELTX
\LabelTeX   51.500   6.200 $t_{n+1}=t_n+h_n$\ELTX
\LabelTeX   78.500   6.200 $t$\ELTX
\LabelTeX    5.000  19.000 $y_n$\ELTX
\LabelTeX    0.800  28.000 $y_{n+1}$\ELTX
\LabelTeX   -2.700  32.000 $z(t_{n+1})$\ELTX
\LabelTeX    6.000  48.000 $y$\ELTX
\EndFig
\bigskip

\noindent Comme le montre le sch\'ema ci-dessous, l'erreur de consistance
n'a a priori que peu de rapport avec l'erreur globale
$\theta_n=\displaystyle\max_{0\le j\le n}|z(t_j)-y_j|$ r\'esultant
d'un calcul de $n$ valeurs successives $y_1,\,\ldots,\,y_n$ \`a partir
de la donn\'ee initiale $y_0=z(t_0)$ (qui est \`a vrai dire la seule
erreur int\'eressant r\'eellement le num\'ericien). On imagine
cependant, et nous reviendrons l\`a-dessus plus en d\'etail au \S~2,
que $|\theta_n|$ sera de l'ordre de grandeur de
$|e_0|+|e_1|+\cdots+|e_{n-1}|$, sous des hypoth\`eses convenables de
r\'egularit\'e pour la fonction~$f$. C'est pourquoi l'\'evaluation de
$e_n$ va gouverner l'\'evaluation de l'erreur globale.

\InsertFig  13.000  73.000  {
1 mm unit /nbpt 12 def 
  0.000   0.000 translate 
gsave 
  0.000  10.000 moveto 
 95.000   0.000 vector 
 10.000   0.000 moveto 
 65.000  90.000 vector 
 20.000  17.000 moveto 
  0.700 disk 
  0.226 setlinewidth
 20.000  17.000 moveto 
 32.000  20.000 lineto 
stroke 
 32.000  20.000 moveto 
  0.500 disk 
 32.000  20.000 moveto 
 44.200  24.000 lineto 
stroke 
 44.200  24.000 moveto 
  0.500 disk 
 44.200  24.000 moveto 
 56.200  29.500 lineto 
stroke 
 56.200  29.500 moveto 
  0.500 disk 
 56.200  29.500 moveto 
 68.200  37.500 lineto 
stroke 
 68.200  37.500 moveto 
  0.500 disk 
  0.339 setlinewidth 
 20.000  17.050 moveto 
begincurve 
[  20.000  17.050    20.079  17.076    20.179  17.109    20.303  17.150 
   20.458  17.202    20.652  17.267    20.895  17.348    21.198  17.452 
   21.577  17.582    21.956  17.714    22.335  17.848    22.714  17.984 
   23.093  18.121    23.472  18.260    23.851  18.401    24.229  18.544 
   24.608  18.688    24.987  18.835    25.366  18.983    25.745  19.133 
   26.124  19.285    26.503  19.439    26.882  19.595    27.261  19.754 
   27.640  19.914    28.019  20.076    28.398  20.240    28.777  20.407 
   29.156  20.575    29.535  20.746    29.914  20.919    30.293  21.094 
   30.672  21.271    31.050  21.451    31.429  21.632    31.808  21.817 
   32.187  22.003    32.566  22.192    32.945  22.383    33.324  22.577 
   33.703  22.773    34.082  22.972    34.461  23.173    34.840  23.377 
   35.219  23.583    35.598  23.792    35.977  24.004    36.356  24.218 
   36.735  24.435    37.114  24.654    37.493  24.876    37.871  25.101 
   38.250  25.329    38.629  25.560    39.008  25.794    39.387  26.030 
   39.766  26.269    40.145  26.512    40.524  26.757    40.903  27.005 
   41.282  27.257    41.661  27.511    42.040  27.768    42.419  28.029 
   42.798  28.293    43.177  28.560    43.556  28.830    43.935  29.104 
   44.314  29.380    44.692  29.661    45.071  29.944    45.450  30.231 
   45.829  30.521    46.208  30.815    46.587  31.113    46.966  31.413 
   47.345  31.718    47.724  32.026    48.103  32.338    48.482  32.653 
   48.861  32.972    49.240  33.295    49.619  33.621    49.998  33.952 
  ] piececurve 
[  50.377  34.286    50.756  34.624    51.135  34.966    51.513  35.312 
   51.892  35.662    52.271  36.016    52.650  36.374    53.029  36.737 
   53.408  37.103    53.787  37.474    54.166  37.848    54.545  38.228 
   54.924  38.611    55.303  38.999    55.682  39.391    56.061  39.787 
   56.440  40.189    56.819  40.594    57.198  41.004    57.577  41.419 
   57.956  41.838    58.334  42.262    58.713  42.691    59.092  43.124 
   59.471  43.563    59.850  44.006    60.229  44.454    60.608  44.907 
   60.987  45.364    61.366  45.827    61.745  46.295    62.124  46.768 
   62.503  47.246    62.882  47.730    63.261  48.218    63.640  48.712 
   64.019  49.212    64.398  49.716    64.777  50.226    65.155  50.742 
   65.534  51.263    65.913  51.789    66.292  52.321    66.671  52.859 
   67.050  53.403    67.429  53.952    67.808  54.507    68.187  55.068 
   68.566  55.635    68.945  56.208    69.324  56.786    69.703  57.371 
   70.082  57.962    70.461  58.559    70.840  59.162    71.219  59.772 
   71.418  60.095   ] piececurve 
endcurve 
  0.113 setlinewidth
 32.080  20.050 moveto 
begincurve 
[  32.080  20.050    32.157  20.083    32.254  20.123    32.374  20.175 
   32.525  20.239    32.713  20.319    32.948  20.421    33.243  20.548 
   33.610  20.709    33.978  20.872    34.346  21.037    34.713  21.203 
   35.081  21.371    35.449  21.541    35.817  21.713    36.184  21.886 
   36.552  22.062    36.920  22.239    37.288  22.418    37.655  22.599 
   38.023  22.782    38.391  22.967    38.759  23.154    39.126  23.343 
   39.494  23.534    39.862  23.727    40.230  23.922    40.597  24.119 
   40.965  24.319    41.333  24.520    41.701  24.724    42.068  24.930 
   42.436  25.138    42.804  25.349    43.172  25.562    43.539  25.777 
   43.907  25.994    44.275  26.214    44.643  26.436    45.010  26.660 
   45.378  26.887    45.746  27.117    46.114  27.349    46.481  27.583 
   46.849  27.820    47.217  28.060    47.585  28.302    47.952  28.547 
   48.320  28.794    48.688  29.044    49.056  29.297    49.423  29.552 
   49.791  29.811    50.159  30.072    50.527  30.336    50.894  30.602 
   51.262  30.872    51.630  31.144    51.998  31.420    52.365  31.698 
   52.733  31.980    53.101  32.264    53.469  32.551    53.836  32.842 
   54.204  33.136    54.572  33.432    54.940  33.732    55.307  34.035 
   55.675  34.342    56.043  34.651    56.411  34.964    56.778  35.281 
   57.146  35.600    57.514  35.923    57.882  36.249    58.249 36.579 
   58.617  36.913    58.985  37.250    59.353  37.590    59.720  37.934 
   60.088  38.281    60.456  38.633    60.824  38.988    61.191  39.346 
  ] piececurve 
[  61.559  39.709    61.927  40.075    62.295  40.445    62.662  40.819 
   63.030  41.196    63.398  41.578    63.766  41.963    64.133  42.353 
   64.501  42.747    64.869  43.144    65.237  43.546    65.604  43.952 
   65.972  44.362    66.340  44.776    66.708  45.195    67.075  45.618 
   67.443  46.045    67.811  46.476    68.179  46.912    68.546  47.352 
   68.914  47.797    69.282  48.247    69.650  48.700    70.017  49.159 
   70.385  49.622    70.753  50.090    71.121  50.562    71.488  51.039 
   71.856  51.521    72.224  52.008    72.591  52.500    72.959  52.996 
   73.327  53.498    73.695  54.004    74.062  54.516    74.430  55.033 
   74.798  55.554    75.166  56.081    75.193  56.120   ] piececurve 
endcurve 
 44.200  23.950 moveto 
begincurve 
[  44.200  23.950    44.273  23.991    44.364  24.043    44.478  24.107 
   44.621  24.188    44.799  24.290    45.022  24.418    45.301  24.579 
   45.579  24.742    45.858  24.906    46.136  25.072    46.415  25.238 
   46.693  25.407    46.972  25.576    47.251  25.747    47.529  25.920 
   47.808  26.094    48.086  26.269    48.365  26.446    48.643  26.624 
   48.922  26.804    49.200  26.986    49.479  27.168    49.758  27.353 
   50.036  27.539    50.315  27.726    50.593  27.915    50.872  28.106 
   51.150  28.298    51.429  28.492    51.708  28.688    51.986  28.885 
   52.265  29.084    52.543  29.284    52.822  29.486    53.100  29.690 
   53.379  29.895    53.657  30.103    53.936  30.312    54.215  30.522 
   54.493  30.735    54.772  30.949    55.050  31.165    55.329  31.382 
   55.607  31.602    55.886  31.823    56.164  32.046    56.443  32.271 
   56.722  32.498    57.000  32.726    57.279  32.957    57.557  33.189 
   57.836  33.424    58.114  33.660    58.393  33.898    58.671  34.138 
   58.950  34.380    59.229  34.624    59.507  34.870    59.786  35.119 
   60.064  35.369    60.343  35.621    60.621  35.875    60.900  36.131 
   61.178  36.389    61.457  36.650    61.736  36.912    62.014  37.177 
   62.293  37.443    62.571  37.712    62.850  37.983    63.128  38.256 
   63.407  38.532    63.685  38.809    63.964  39.089    64.243  39.371 
   64.521  39.655    64.800  39.942    65.078  40.231    65.357  40.522 
   65.635  40.815    65.914  41.111    66.192  41.409    66.471  41.709 
  ] piececurve 
[  66.750  42.012    67.028  42.317    67.307  42.625    67.585  42.935 
   67.864  43.247    68.142  43.562    68.421  43.880    68.700  44.200 
   68.978  44.522    69.257  44.847    69.535  45.174    69.814  45.504 
   70.092  45.837    70.371  46.172    70.649  46.510    70.928  46.850 
   71.207  47.193    71.485  47.539    71.764  47.887    72.042  48.238 
   72.321  48.592    72.599  48.948    72.878  49.308    73.156  49.670 
   73.435  50.034    73.714  50.402    73.992  50.772    74.271  51.145 
   74.549  51.521    74.828  51.900    75.106  52.282    75.385  52.667 
   75.663  53.054    75.942  53.445    76.221  53.838    76.499  54.235 
   76.778  54.634    77.056  55.037    77.335  55.442    77.613  55.851 
   77.796  56.120   ] piececurve 
endcurve 
 56.200  29.500 moveto 
begincurve 
[  56.200  29.500    56.267  29.551    56.350  29.615    56.454  29.695 
   56.585  29.795    56.748  29.921    56.951  30.079    57.206  30.278 
   57.460  30.479    57.715  30.681    57.970  30.885    58.224  31.091 
   58.479  31.298    58.733  31.506    58.988  31.716    59.242  31.928 
   59.497  32.141    59.752  32.356    60.006  32.573    60.261  32.791 
   60.515  33.011    60.770  33.233    61.024  33.456    61.279  33.681 
   61.534  33.907    61.788  34.136    62.043  34.366    62.297  34.597 
   62.552  34.831    62.806  35.066    63.061  35.303    63.316  35.542 
   63.570  35.783    63.825  36.025    64.079  36.269    64.334  36.515 
   64.588  36.763    64.843  37.013    65.097  37.264    65.352  37.518 
   65.607  37.773    65.861  38.030    66.116  38.289    66.370  38.550 
   66.625  38.813    66.879  39.078    67.134  39.345    67.389  39.613 
   67.643  39.884    67.898  40.157    68.152  40.431    68.407  40.708 
   68.661  40.987    68.916  41.267    69.171  41.550    69.425  41.835 
   69.680  42.122    69.934  42.411    70.189  42.702    70.443  42.995 
   70.698  43.291    70.953  43.588    71.207  43.888    71.462  44.189 
   71.716  44.493    71.971  44.799    72.225  45.108    72.480  45.418 
   72.735  45.731    72.989  46.046    73.244  46.363    73.498  46.682 
   73.753  47.004    74.007  47.328    74.262  47.654    74.516  47.983 
   74.771  48.314    75.026  48.647    75.280  48.983    75.535  49.321 
   75.789  49.662    76.044  50.004    76.298  50.350    76.553  50.697 
  ] piececurve 
[  76.808  51.047    77.062  51.400    77.317  51.755    77.571  52.112 
   77.826  52.472    78.080  52.835    78.335  53.200    78.590  53.567 
   78.844  53.938    79.099  54.310    79.353  54.685    79.608  55.063 
   79.862  55.444    80.117  55.827    80.311  56.120   ] piececurve 
endcurve 
  0.226 setlinewidth
 32.000  22.000 moveto 
  0.500 disk 
 32.000  22.000 moveto 
 44.200  27.500 lineto 
stroke 
 44.200  27.500 moveto 
  0.500 disk 
 44.200  29.100 moveto 
  0.500 disk 
 44.200  29.100 moveto 
 56.200  37.500 lineto 
stroke 
 56.200  37.500 moveto 
  0.500 disk 
 56.200  39.800 moveto 
  0.500 disk 
 56.200  39.800 moveto 
 68.200  52.000 lineto 
stroke 
 68.200  52.000 moveto 
  0.500 disk 
 68.200  55.000 moveto 
  0.500 disk
  0.113 setlinewidth
 32.000  20.000 moveto 
 32.000  22.000 lineto 
stroke 
 44.200  27.500 moveto 
 44.200  29.100 lineto 
stroke 
 56.200  37.500 moveto
 56.200  39.500 lineto 
stroke  
 68.200  52.000 moveto 
 68.200  55.000 lineto 
stroke 
[ 0.400   0.600 ] 0 setdash
 10.000  20.000 moveto 
 32.000  20.000 lineto 
stroke
 10.000  22.000 moveto 
 32.000  22.000 lineto
stroke 
 10.000  27.500 moveto 
 44.200  27.500 lineto
stroke
 10.000  29.100 moveto 
 44.200  29.100 lineto
stroke
 10.000  37.500 moveto 
 56.200  37.500 lineto
stroke
 10.000  39.800 moveto 
 56.200  39.800 lineto
stroke
 10.000  52.000 moveto  
 68.200  52.000 lineto 
stroke
 10.000  55.000 moveto  
 90.000  55.000 lineto 
stroke
 10.000  37.500 moveto  
 90.000  37.500 lineto 
stroke
 3 setdashtype 
 20.000  10.000 moveto 
 20.000  17.000 lineto 
stroke 
 32.000  10.000 moveto 
 32.000  20.000 lineto 
stroke 
 44.200  10.000 moveto 
 44.200  27.500 lineto 
stroke 
 56.200  10.000 moveto 
 56.200  37.500 lineto 
stroke 
 68.200  10.000 moveto 
 68.200  52.000 lineto 
stroke 
 20.000  10.000 moveto 
  0.500 disk 
 32.000  10.000 moveto 
  0.500 disk 
 44.200  10.000 moveto 
  0.500 disk 
 56.200  10.000 moveto 
  0.500 disk 
 68.200  10.000 moveto 
  0.500 disk 
0 setdashtype
  1.400 /arrowlgth exch def 
 90.000  47.000 moveto 
  8.000  90.000 vector 
 90.000  47.000 moveto 
  9.500 -90.000 vector
grestore 
}
\LabelTeX  19.200   6.400 $t_0$\ELTX
\LabelTeX  31.200   6.400 $t_1$\ELTX
\LabelTeX  43.200   6.400 $t_2$\ELTX
\LabelTeX  55.200   6.400 $t_3$\ELTX
\LabelTeX  67.200   6.400 $t_4$\ELTX
\LabelTeX  94.000   6.400 $t$\ELTX
\LabelTeX   6.200  20.600 $e_0$\ELTX
\LabelTeX   6.200  27.900 $e_1$\ELTX
\LabelTeX   6.200  38.400 $e_2$\ELTX
\LabelTeX   6.200  53.000 $e_3$\ELTX
\LabelTeX   6.200  63.000 $y$\ELTX
\LabelTeX  16.000  15.500 $y_0$\ELTX
\LabelTeX  33.200  17.900 $y_1$\ELTX
\LabelTeX  45.200  21.900 $y_2$\ELTX
\LabelTeX  57.200  27.900 $y_3$\ELTX
\LabelTeX  69.200  34.900 $y_4$\ELTX
\LabelTeX  70.700  62.000 $z$\ELTX
\LabelTeX  73.000  57.900 $z_1$\ELTX
\LabelTeX  76.800  57.900 $z_2$\ELTX
\LabelTeX  81.000  57.900 $z_3$\ELTX
\LabelTeX  91.200  45.900 $\theta_4$\ELTX
\EndFig

\smallskip\noindent
[Les fonctions $z$, $z_1$, $z_2$, $z_3$ repr\'esentent ici les solutions
exactes passant par les points $(t_0,y_0)$ et $(t_j,y_j)$, $j=1,\,2,\,3$].

\subsection 1.2. Retour sur la m\'ethode d'Euler|

\noindent Soit $z$ une solution exacte de l'\'equation (E). On a au premier ordre l'approximation
$$
z(t_{n+1})=z(t_n+h_n)\simeq z(t_n)+h_nz'(t_n)=z(t_n)+h_nf(t_n,z(t_n)).
$$
Comme on l'a d\'ej\`a vu au chapitre V, ceci conduit \`a l'algorithme
$$
\cases{
y_{n+1}=y_n+h_nf(t_n,y_n)\hfill\cr
\noalign{\vskip3pt}
t_{n+1}=t_n+h_n.\hfill\cr
}
$$
Par d\'efinition de l'erreur de consistance, on a $e_n=z(t_n+h_n)-y_{n+1}$
o\`u 
$$
y_{n+1}=z(t_n)+h_nf(t_n,z(t_n))=z(t_n)+h_nz'(t_n).
$$
La formule de Taylor-Lagrange donne
$$
e_n=z(t_n+h_n)-(z(t_n)+h_nz'(t_n))=\frac{1}{2}\ h^2_nz''(t_n)+o(h_n^2),
$$
pourvu que $z$ soit de classe $C^2$. C'est bien le cas si $f$ 
est de classe $C^1$, et on sait alors que
$$
z''(t)=f^{[1]}(t,z(t))\qquad\hbox{o\`u}\qquad f^{[1]}=f'_t+f'_y f.
$$
On en d\'eduit par cons\'equent
$$
e_n=\frac{1}{2}\ h^2_nf^{[1]}(t_n,y_n)+o(h_n^2).
$$
Cette erreur en $h^2_n$ est relativement importante, \`a moins que le
pas $h_n$ ne soit choisi tr\`es petit, ce qui augmente
consid\'erablement le volume des calculs \`a effectuer. On va donc
essayer de construire des m\'ethodes permettant de r\'eduire l'erreur
de consistance $e_n$.

\subsection 1.3. M\'ethode de Taylor d'ordre $p$|

\noindent Supposons que $f$ soit de classe $C^p$. Alors toute solution
exacte $z$ est de classe $C^{p+1}$, et sa d\'eriv\'ee $k$-i\`eme est
$z^{(k)}(t)=f^{[k-1]}(t,z(t))$. La formule de Taylor d'ordre $p$
implique
$$
z(t_n+h_n)=z(t_n)+\sum_{k=1}^p\frac{1}{k!}\ h_n^kf^{[k-1]}(t_n,z(t_n))
+o(h_n^p).
$$
Lorsque $h_n$ est assez petit, l'approximation est d'autant meilleure que
$p$ est plus grand. On est donc amen\'e \`a consid\'erer l'algorithme 
suivant, appel\'e m\'ethode de Taylor d'ordre $p$~:
$$
\cases{
\displaystyle y_{n+1}=y_n+\sum_{k=1}^p\dfrac{1}{k!}\ h^k_nf^{[k-1]}(t_n,y_n)
\hfill\cr
\noalign{\vskip8pt}
t_{n+1}=t_n+h_n.\hfill\cr}
$$
D'apr\`es la d\'efinition g\'en\'erale des m\'ethodes \`a un pas,
cet algorithme correspond au choix $\Phi(t,y,h)=\sum_{k=1}^p
{1\over k!}\ h^{k-1}f^{[k-1]}(t,y)$.  Calculons l'erreur de consistance
$e_n$. En supposant $y_n=z(t_n)$, la formule de Taylor d'ordre $p+1$
donne
$$
\eqalign{
e_n&=z(t_{n+1})-y_{n+1}=z(t_n+h_n)-\sum_{k=0}^p\frac{1}{k!}\ h^k_n
z^{(k)}(t_n)\cr
&=\frac{1}{(p+1)!}\ h^{p+1}_n f^{[p]}(t_n,y_n)+o(h_n^{p+1}).}
$$
L'erreur est donc maintenant de l'ordre de $h^{p+1}_n$. On dira d'une
mani\`ere g\'en\'erale qu'une m\'ethode est {\it d'ordre $p$} si
l'erreur de consistance est en $h^{p+1}_n$ chaque fois que $f$ est de
classe $C^p$ au moins. La m\'ethode d'Euler est le cas particulier
$p=1$ de la m\'ethode de Taylor.

\bigskip
\state Remarque -- |  Dans la pratique, la m\'ethode de Taylor souffre
de deux inconv\'enients graves qui en font g\'en\'eralement
d\'econseiller l'utilisation pour $p\geq 2$~:

\medskip
{\parindent1pc
\item{$\bullet$} Le calcul des quantit\'es $f^{[k]}$ est souvent complexe et
co\^uteux en temps machine. Il faut aussi pouvoir \'evaluer explicitement
$f^{[k]}$, ce qui n'est pas toujours le cas (par exemple, si $f$ est 
une donn\'ee exp\'erimentale discr\'etis\'ee).

\item{$\bullet$} La m\'ethode suppose {\it a priori} que $f$ soit
tr\`es r\'eguli\`ere~; les erreurs risquent donc de ne pas pouvoir
\^etre contr\^ol\'ees si certaines d\'eriv\'ees de $f$ pr\'esentent des
discontinuit\'es ou une mauvaise continuit\'e (pentes \'elev\'ees).
\par}

\subsection 1.4. M\'ethode du point milieu|

\noindent Cette m\'ethode est d\'ecrite par le sch\'ema suivant.

\InsertFig 19.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 
  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 
 75.000  34.000 moveto 
  0.000  90.000 vector 
 75.000  36.000 lineto 
stroke 
 75.000  36.000 moveto 
  0.000 -90.000 vector 
[ 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 
 11.500  36.000 moveto 
 75.000  36.000 lineto 
stroke 
 60.000  34.000 moveto 
 75.000  34.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}
\LabelTeX    6.000  56.000 $y$\ELTX
\LabelTeX   -3.600  35.000 $z(t+h)$\ELTX
\LabelTeX    2.700  24.000 $z(t)$\ELTX
\LabelTeX   29.200   6.200 $t$\ELTX
\LabelTeX   41.200   6.200 $t+h/2$\ELTX
\LabelTeX   58.000   6.200 $t+h$\ELTX
\LabelTeX   88.000   6.200 $t$\ELTX
\LabelTeX   76.200  34.000 $\varepsilon$\ELTX
\LabelTeX   71.000  50.500 $z$\ELTX
\EndFig
\vskip3mm

\noindent L'id\'ee est que la corde de la fonction $z$ sur $[t,t+h]$ a
une pente voisine de $z'\big(t+\dfrac{h}{2}\big)$, alors que dans la
m\'ethode d'Euler on approxime brutalement cette pente par $z'(t)$. On
\'ecrit donc~:
$$
z(t+h)\simeq z(t)+hz'\Big(t+\frac{h}{2}\Big).\eqno(*)
$$
Si $z$ est de classe $C^3$, il vient
$$
\eqalign{
z(t+h)&=z(t)+hz'(t)+\frac{1}{2}\ h^2z''(t)+\frac{1}{6}\ h^3z'''(t)+o(h^3),\cr
z'\Big(t+\frac{h}{2}\Big)&=\qquad\quad\kern6pt z'(t)+\frac{1}{2}\ hz''(t)+\frac{1}{8}\
h^2z'''(t)+o(h^2).\cr }
$$
L'erreur commise est donc
$$
\varepsilon=z(t+h)-z(t)-hz'\Big(t+\frac{h}{2}\Big)=\frac{1}{24}\ h^3z'''(t)+o(h^3),
$$
soit une erreur en $h^3$ au lieu de $h^2$ dans la m\'ethode d'Euler. On par ailleurs
$$
z'\Big(t+\frac{h}{2}\Big)=f\Big(t+\frac{h}{2},z\Big(t+\frac{h}{2}\Big)\Big).
$$
Comme la valeur de $z\Big(t+\dfrac{h}{2}\Big)$ n'est pas connue, on l'approxime par
$$
z\Big(t+\frac{h}{2}\Big)\simeq z(t)+\frac{h}{2}\ f(t,z(t)),\eqno(**)
$$
d'o\`u en d\'efinitive
$$
z(t+h)\simeq z(t)+hf\Big(t+\frac{h}{2},z(t)+\frac{h}{2}\,f(t,z(t))\Big).
$$
L'algorithme du point milieu est associ\'e au choix
$$
\Phi(t,y,h)=f\Big(t+\frac{h}{2},y+\frac{h}{2}\,f(t,y)\Big)
$$
et donne lieu au sch\'ema num\'erique
$$
\cases{
y_{n+\frac{1}{2}}=y_n+\displaystyle\dfrac{h_n}{2}\ f(t_n,y_n)\hfill\cr
\noalign{\vskip3pt}
p_n=f\Big(t_n+\displaystyle\dfrac{h_n}{2},y_{n+\fracss{1}{2}}\Big)\cr
\noalign{\vskip8pt}
y_{n+1}=y_n+h_np_n\hfill\cr
\noalign{\vskip11pt}
t_{n+1}=t_n+h_n.\hfill\cr
\noalign{\vskip4pt}
}
$$
\vskip0pt\noindent
Calculons l'erreur de consistance~: $e_n=z(t_{n+1})-y_{n+1}$, avec $y_n=z(t_n)$. On a
$e_n=\varepsilon_n+\varepsilon'_n$ o\`u les erreurs 
$$
\eqalign{
\varepsilon_n&=z(t_{n+1})-z(t_n)-h_nz'\Big(t_n+\frac{h_n}{2}\Big),\cr
\varepsilon'_n&=h_nz'\Big(t_n+\frac{h_n}{2}\Big)-(y_{n+1}-z(t_n))\cr
&=h_n\left(f\Big(t_n+\frac{h_n}{2},z\Big(t_n+\frac{h_n}{2}\Big)\Big)-f\Big(t_n+\frac{h_n}{2},
y_{n+\frac{1}{2}}\Big)\right)\cr
}
$$
proviennent respectivement des approximations $(*)$ et $(**)$. D'apr\`es le calcul fait plus haut
$$
\varepsilon_n=\frac{1}{24}\ h^3_nz'''(t_n)+o(h_n^3)=\frac{1}{24}\ h^3_nf^{[2]}(t_n,y_n)+o(h^3_n).
$$
D'autre part
$$
\eqalign{
z\Big(t_n+\frac{h_n}{2}\Big)-y_{n+\frac{1}{2}}&=z\Big(t_n+\frac{h_n}{2}\Big)-\Big(z(t_n)+
\frac{h_n}{2}\,z'(t_n)\Big)\cr
&=\frac{1}{8}\ h^2_nz''(t_n)+o(h^2_n)=\frac{1}{8}\ h^2_nf^{[1]}(t_n,y_n)+o(h^2_n).\cr
}
$$
D'apr\`es le th\'eor\`eme des accroissements finis appliqu\'e en $y$, on a
$$
\eqalign{
f\Big(t_n+\frac{h_n}{2},z\Big(t_n+\frac{h_n}{2}\Big)\Big)&-\
f\Big(t_n+\frac{h_n}{2},y_{n+\frac{1}{2}} \Big)\cr
&=f'_y\Big(t_n+\frac{h_n}{2},c_n\Big)\Big(z\Big(t_n+\frac{h_n}{2}\Big)-y_{n+\frac{1}{2}}\Big)\cr
&=\Big(f'_y(t_n,y_n)+o(h_n)\Big)\Big(\frac{1}{8}\ h^2_nf^{[1]}(t_n,y_n)+o(h^2_n)\Big)\cr
&=\frac{1}{8}\ h^2_n f'_yf^{[1]}(t_n,y_n)+o(h^2_n),\cr
}
$$
d'o\`u
$$
\varepsilon'_n=\frac{1}{8}\ h^3_nf'_yf^{[1]}(t_n,y_n)+o(h^3_n).
$$
On en d\'eduit
$$
e_n=\varepsilon_n+\varepsilon'_n=\frac{1}{24}\ h^3_n\Big(f^{[2]}+3f'_yf^{[1]}\Big)(t_n,y_n)+o(h^3_n).
$$
La m\'ethode du point milieu est donc d'ordre 2.

\subsection 1.5.* M\'ethode du point milieu modifi\'e|

\noindent Si on observe les algorithmes pr\'ec\'edents, on voit que la seule
op\'eration \'eventuel\-lement co\^uteuse en temps de calcul est
l'\'evaluation de la fonction $f(t,y)$, le reste consistant en un
petit nombre d'additions ou de multiplications. On mesure donc le
co\^ut d'une m\'ethode d'ordre donn\'e {\it par le nombre
d'\'evaluations de la fonction} $f$ qu'elle r\'eclame \`a chaque
pas. Pour des m\'ethodes d'ordres diff\'erents la comparaison ne tient
pas, puisqu'une m\'ethode d'ordre plus \'elev\'e exige \`a pr\'ecision
\'egale un nombre de pas nettement inf\'erieur.
\smallskip

Dans la m\'ethode du point milieu, on va modifier le calcul successif
de $f(t_n,y_n)$ et de la pente interm\'ediaire $p_n=f\big(t_n+
\dfrac{h_n}{2},y_{n+\frac{1}{2}}\big)$ en introduisant l'algorithme 
suivant, qui fait l'\'economie de l'\'evaluation de $f(t_n,y_n)$~:
$$
\cases{
\widetilde y_{n+\frac{1}{2}}=\widetilde y_n+
\displaystyle\dfrac{h_n}{2}\ \widetilde p_{n-1}\hfill\cr
\noalign{\vskip3pt}
\widetilde p_n=f\Big(t_n+
\displaystyle\dfrac{h_n}{2},\widetilde y_{n+\frac{1}{2}}\Big)
\hfill\cr
\noalign{\vskip8pt}
\widetilde y_{n+1}=\widetilde y_n+h_n\widetilde p_n\hfill\cr
\noalign{\vskip11pt}
t_{n+1}=t_n+h_n.\hfill\cr
\noalign{\vskip4pt}
}
$$
On a donc modifi\'e l\'eg\`erement le calcul de $y_{n+\frac{1}{2}}$
en rempla\c cant la pente $f(t_n,y_n)$ par la pente $\widetilde
p_{n-1}$ calcul\'ee \`a l'\'etape ant\'erieure. Il s'ensuit
naturellement que les valeurs $y_n$ sont elles aussi modifi\'ees en
des valeurs $\widetilde y_n$.

\medskip
\state Remarque -- |  Le d\'emarrage (\'etape $n=0$) pr\'esente une 
difficult\'e car la pente $\widetilde p_{-1}$ n'a pas \'et\'e 
\'evalu\'ee. On r\'esout cette difficult\'e en initialisant
$\widetilde p_{-1}=f(t_0,y_0)$. On observera que la m\'ethode du point
milieu modifi\'e est en fait une m\'ethode \`a 2 pas (les \'etapes $n$
et $n-1$ sont utilis\'ees pour calculer $\widetilde y_{n+1}$).\hfill
\blacksquare
\bigskip

\'Evaluons maintenant l'erreur de consistance $\widetilde
e_n=z(t_{n+1})-\widetilde y_{n+1}$, en supposant $\widetilde
y_n=z(t_n)$. On peut \'ecrire
$$
\widetilde e_n=(z(t_{n+1})-y_{n+1})+(y_{n+1}-\widetilde y_{n+1})=e_n+\varepsilon''_n,
$$
o\`u $e_n$ est l'erreur de consistance de la m\'ethode du point
milieu standard (on suppose donc aussi $y_n=z(t_n)$ pour la calculer).
Il vient
$$
\eqalign{
\varepsilon''_n=y_{n+1}-\widetilde y_{n+1}&=h_n\Big(f\Big(t_n+\frac{h_n}{2},y_{n+\frac{1}{2}}\Big)
-f\Big(t_n+\frac{h_n}{2},\widetilde y_{n+\frac{1}{2}}\Big)\Big).\cr
y_{n+\frac{1}{2}}-\widetilde y_{n+\frac{1}{2}}&=\frac{h_n}{2}\ \Big(f(t_n,y_n)-\widetilde
p_{n-1}\Big)\cr
&=\frac{h_n}{2}\ \Big(f(t_n,y_n)-f\Big(t_{n-1}+\frac{h_{n-1}}{2},\widetilde
y_{n-\frac{1}{2}}\Big)\Big).\cr
}
$$
Or\quad $t_n-\big(t_{n-1}+\dfrac{h_{n-1}}{2}\big)=\dfrac{h_{n-1}}{2}$\quad et
$$
\eqalign{
y_n-\widetilde y_{n-\frac{1}{2}}&=y_n-\Big(\widetilde y_{n-1}+\frac{h_{n-1}}{2}\ \widetilde
p_{n-2}\Big)\cr 
&=y_n-\Big(\widetilde y_n-h_{n-1}\widetilde p_{n-1}+\frac{h_{n-1}}{2}\ \widetilde p_{n-2}\Big)\cr
&=h_{n-1}\Big(\widetilde p_{n-1}-\frac{1}{2}\ \widetilde p_{n-2}\Big)\cr
&=\frac{1}{2}\ h_{n-1}f(t_n,y_n)+o(h_{n-1})~;\cr
}
$$
\`a la troisi\`eme ligne on utilise le fait que $y_n=\widetilde y_n=z(t_n)$, et \`a la quatri\`eme
le fait que $\widetilde p_{n-i}=f(t_{n-i}+h_{n-i}/2,\widetilde y_{n-i}+1/2)$ converge vers
$f(t_n,y_n)$ pour $i=1,2$ lorsque $h_{\max}$ tend vers 0.
Gr\^ace \`a la formule de Taylor pour les fonctions de 2 variables, il vient
$$
\eqalign{
f(t_n,y_n)-f\Big(t_{n-1}&{}+\frac{h_{n-2}}{2},\widetilde y_{n-\frac{1}{2}}\Big)
\cr
&=\frac{h_{n-1}}{2}\ f'_t(t_n,y_n)+\frac{1}{2}\ h_{n-1}f(t_n,y_n)
f'_y(t_n,y_n)+o(h_{n-1})\cr
&=\frac{1}{2}\ h_{n-1}(f'_t+ff'_y)(t_n,y_n)+o(h_{n-1})\cr
&=\frac{1}{2}\ h_{n-1}f^{[1]}(t_n,y_n)+o(h_{n-1}),\cr}
$$
d'o\`u
$$
y_{n+\frac{1}{2}}-\widetilde y_{n+\frac{1}{2}}=\frac{1}{4}
h_nh_{n-1}f^{[1]}(t_n,y_n)+o(h_nh_{n-1}). 
$$
On en d\'eduit finalement
$$
\varepsilon''_n=h_nf'_y\Big(t_n+\frac{h_n}{2},c_n\Big)
\big(y_{n+\frac{1}{2}}-\widetilde
y_{n+\frac{1}{2}}\big)
=\frac{1}{4}\ h^2_nh_{n-1}\big(f'_yf^{[1]}\big)(t_n,y_n)+o(h^2_nh_{n-1}),
$$
d'o\`u l'erreur de consistance
$$
\widetilde e_n=\frac{1}{24}\ h^3_n\Big(f^{[2]}+3f'_yf^{[1]}\Big)
(t_n,y_n)+\frac{1}{4}\ h^2_n h_{n-1}
(f'_yf^{[1]})(t_n,y_n)+o(h^3_n+h^2_nh_{n-1}).
$$
La m\'ethode du point milieu modifi\'e est donc encore une m\'ethode d'ordre 2 (mais ce n'est pas
une m\'ethode \`a un pas~!).

\section 2. \'Etude g\'en\'erale des m\'ethodes \`a un pas|

\subsection 2.1. M\'ethodes consistantes, stables et convergentes|

La premi\`ere notion que nous introduisons a trait au probl\`eme
de l'accumulation des erreurs de consistance (accumulation purement
th\'eorique dans le sens o\`u on tient pas compte du fait que la
solution calcul\'ee s'\'ecarte de la solution exacte, cf.\ sch\'emas du
\S~1.1).

\bigskip
\majdeb
\state D\'efinition 1 -- |  
{\it On dit que la m\'ethode est {\rm consistante} si pour toute
solution exacte $z$ la somme des erreurs de consistance relatives
\`a $z$, soit $\displaystyle\sum_{0\leq n\leq N}|e_n|$, tend vers 0
quand $h_{\max}$ tend vers 0.}
\par\majfin\par

\bigskip Une autre notion fondamentale est la notion de stabilit\'e.
Dans la pratique, le calcul r\'ecurrent des points $y_n$ est en effet
ent\^ach\'e d'erreurs d'arrondi $\varepsilon_n$. Pour que les calculs
soient significatifs, il est indispensable que la propagation de ces
erreurs reste contr\^olable. On est amen\'e \`a la d\'efinition
suivante.

\bigskip
\majdeb
\state D\'efinition 2 -- |  {\it On dit que la m\'ethode est {\rm
stable} s'il existe une constante $S\geq 0$, appel\'ee {\rm constante de 
stabilit\'e}, telle que pour toutes suites $(y_n)$,
$(\widetilde y_n)$ d\'efinies par
$$
\eqalign{
y_{n+1}&=y_n+h_n\Phi(t_n,y_n,h_n),\qquad\quad\kern2pt 0\leq n<N\cr
\widetilde y_{n+1}&=\widetilde y_n+h_n\Phi(t_n,\widetilde y_n,h_n)
+\varepsilon_n,\quad 0\leq n<N\cr}
$$
on ait
$$
\max_{0\leq n\leq N}|\widetilde y_n-y_n|\leq 
S\Big(|\widetilde y_0-y_0|+\sum_{0\leq n<N}|\varepsilon_n|\Big). 
$$} 
\vskip-\baselineskip
\majfin
\vskip\baselineskip

Autrement dit, une petite erreur initiale $|\widetilde
y_0-y_0|$ et de petites erreurs d'arrondi $\varepsilon_n$ dans le
calcul r\'ecurrent des $\widetilde y_n$ provoquent une erreur finale
$\max|\widetilde y_n-y_n|$ contr\^olable. Une derni\`ere notion
importante en pratique est la suivante.

\bigskip
\majdeb
\state D\'efinition 3 -- |  {\it On dit que la m\'ethode est {\rm convergente}
si pour toute solution exacte $z$, la suite $y_n$ telle que 
$y_{n+1}=y_n+h_n\Phi(t_n,y_n,h_n)$ v\'erifie
$$
\max_{0\leq n\leq N}|y_n-z(t_n)|\rightarrow 0
$$
quand $y_0\rightarrow z(t_0)$ et quand $h_{\max}\rightarrow 0$.
}
\par\majfin\par

\bigskip La quantit\'e $\displaystyle\max_{0\leq n\leq N}|y_n-z(t_n)|$
s'appelle {\it l'erreur globale} (de la suite $y_n$ calcul\'ee par
rapport \`a la solution exacte $z$). C'est \'evidemment cette erreur
qui importe dans la pratique.

\bigskip
\state Calcul de l'erreur globale -- |  Posons $\widetilde y_n=z(t_n)$. Par
d\'efinition de l'erreur de consistance (cf. \S\ 2.1) on a 
$$
\widetilde y_{n+1}=\widetilde y_n+h_n\Phi(t_n,\widetilde y_n,h_n)+e_n.
$$
Si la m\'ethode est stable, de constante de stabilit\'e $S$,
l'erreur globale est donc

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

\bigskip
\majdeb
\state Corollaire -- |  {\it Si la m\'ethode est stable et consistante, elle
est convergente.}
\par\majfin\par

\bigskip
En effet $\sum_{0\leq n<N}|e_n|$ tend vers 0 quand $h_{\max}$ tend vers 0, 
puisque la m\'ethode est consistante par hypoth\`ese.

\subsection 2.2. Condition n\'ecessaire et suffisante de consistance|

\noindent Soit $z$ une solution exacte de l'\'equation (E) et soient 
$$
e_n=z(t_{n+1})-z(t_n)-h_n\Phi(t_n,z(t_n),h_n)
$$
les erreurs de consistance correspondantes. D'apr\`es le th\'eor\`eme des
accroissements finis, il existe $c_n\in\ ]t_n,t_{n+1}[$ tel que
$$
z(t_{n+1})-z(t_n)=h_nz'(c_n)=h_nf(c_n,z(c_n))
$$
d'o\`u
$$
e_n=h_n(f(c_n,z(c_n))-\Phi(t_n,z(t_n),h_n))=h_n(\alpha_n+\beta_n)
$$
avec
$$
\eqalign{
\alpha_n&=f(c_n,z(c_n))-\Phi(c_n,z(c_n),0),\cr
\beta_n&=\Phi(c_n,z(c_n),0)-\Phi(t_n,z(t_n),h_n).\cr
}
$$
Comme la fonction $(t,h)\mapsto\Phi(t,z(t),h)$ est continue sur
$[t_0,t_0+T]\times[0,\delta]$ qui est compact, elle y est
uniform\'ement continue. Par cons\'equent, pour tout $\varepsilon>0$,
il existe $\eta>0$ tel que $h_{\max}\leq\eta\rightarrow|\beta_n|\leq
\varepsilon$. Pour $h_{\max}\leq\eta$ on a donc
$$
\left|\sum_{0\leq n<N}|e_n|-\sum_{0\leq n<N}h_n|\alpha_n|\right|\leq 
\sum_{0\leq n<N}h_n|\beta_n|\leq\varepsilon\sum h_n=T{\varepsilon}.
$$
On en d\'eduit
$$
\eqalign{
\lim_{h_{\max}\rightarrow 0}\sum_{0\leq n<N}|e_n|&=\lim_{h_{\max}\rightarrow 0}
\sum_{0\leq n<N}h_n|\alpha_n|\cr
&=\int^{t_0+T}_{t_0}|f(t,z(t))-\Phi(t,z(t),0)|dt\cr
}
$$
car $\displaystyle \sum h_n|\alpha_n|$ est une somme de Riemann de
l'int\'egrale pr\'ec\'edente. Par d\'efinition, la m\'ethode est
consistante si et seulement si $\displaystyle\lim\sum|e_n|=0$ pour
toute solution exacte $z$. On en d\'eduit~:

\bigskip
\majdeb
\state Th\'eor\`eme -- |  {\it La m\'ethode \`a $1$ pas d\'efinie par la 
fonction $\Phi$ est consistante si et seulement si
$$
\forall(t,y)\in [t_0,t_0+T]\times\bR,\quad\Phi(t,y,0)=f(t,y).
$$
} 
\vskip-\baselineskip
\majfin
\vskip\baselineskip

\noindent Il r\'esulte de ce th\'eor\`eme que les m\'ethodes \`a un pas
d\'ej\`a mentionn\'ees sont bien consistantes. 

\subsection 2.3. Condition suffisante de stabilit\'e|

\noindent Pour pouvoir majorer l'erreur globale d\'ecrite au \S\ 2.1, il faut
savoir estimer d'une part la constante de stabilit\'e $S$, et d'autre
part la somme $\sum_{0\leq n<N}|e_n|$.  Le r\'esultat
suivant permet d'\'evaluer $S$.

\bigskip
\majdeb
\state Th\'eor\`eme -- |  {\it Pour que la m\'ethode soit stable, il suffit
que la fonction $\Phi$ soit lipschitzienne en $y$, c'est-\`a-dire qu'il 
existe une constante $\Lambda\geq 0$ telle que
\vskip6pt
$\forall t\in[t_0,t_0+T]$, $\forall(y_1,y_2)\in\bR^2$, $\forall h\in\bR$ on ait
$$
|\Phi(t,y_1,h)-\Phi(t,y_2,h)|\leq \Lambda|y_1-y_2|.
$$
Dans ce cas, on peut prendre pour constante de stabilit\'e~
$S=e^{\Lambda T}$.} 
\par\majfin\bigskip

\noindent{\it D\'emonstration.} Consid\'erons deux suites $(y_n)$, 
$(\widetilde y_n)$ telles que
$$
\eqalign{
y_{n+1}&=y_n+h_n\Phi(t_n,y_n,h_n),\cr
\widetilde y_{n+1}&=\widetilde y_n+h_n\Phi(t_n,\widetilde y_n,h_n)
+\varepsilon_n.\cr}
$$
Par diff\'erence, on obtient
$$
|\widetilde y_{n+1}-y_{n+1}|\leq |\widetilde y_n-y_n|+h_n\Lambda|\widetilde
y_n-y_n|+|\varepsilon_n|. 
$$
En posant $\theta_n=|\widetilde y_n-y_n|$, il vient
$$
\theta_{n+1}\leq (1+\Lambda h_n)\theta_n+|\varepsilon_n|.
$$

\bigskip
\majdeb
\state Lemme de Gronwall {\rm(cas discret)} -- | {\it Soient des suites 
$h_n,\theta_n\geq 0$ et
$\varepsilon_n\in\bR$ telles que $\theta_{n+1}\leq (1+\Lambda h_n)\theta_n
+|\varepsilon_n|$. Alors
$$
\theta_n\leq e^{\Lambda(t_n-t_0)}\theta_0+\sum_{0\leq i\leq n-1}
e^{\Lambda(t_n-t_{i+1})}|\varepsilon_i|. 
$$}
\vskip-\baselineskip
\majfin
\vskip\baselineskip

Le lemme se v\'erifie par r\'ecurrence sur $n$. Pour $n=0$,
l'in\'egalit\'e se r\'eduit \`a $\theta_0\leq \theta_0$. Supposons
maintenant l'in\'egalit\'e vraie \`a l'ordre $n$. On observe que
$$
1+\Lambda h_n\leq e^{\Lambda h_n}=e^{\Lambda(t_{n+1}-t_n)}.
$$
Par hypoth\`ese on a
$$
\eqalign{
\theta_{n+1}&\leq e^{\Lambda(t_{n+1}-t_n)}\theta_n+|\varepsilon_n|\cr
&\leq e^{\Lambda(t_{n+1}-t_0)}\theta_0+\sum_{0\leq i\leq
n-1}e^{\Lambda(t_{n+1}-t_{i+1})}|\varepsilon_i|+|\varepsilon_n|.\cr 
}
$$
L'in\'egalit\'e cherch\'ee s'ensuit \`a l'ordre $n+1$.\hfill\blacksquare

\bigskip
Comme $t_n-t_0\leq T$ et $t_{n}-t_{i+1}\leq T$, le lemme de Gronwall implique
$$
\max_{0\leq n\leq N}\theta_n\leq e^{\Lambda T}\Big(\theta_0+\sum_{0\leq i\leq N-1}|\varepsilon_i|\Big)
$$
Par d\'efinition de $\theta_n$, on a donc
$$
\max_{0\leq n\leq N}|\widetilde y_n-y_n|\leq e^{\Lambda T}\Big(|\widetilde y_0-y_0|+\sum_{0\leq n\leq
N} |\varepsilon_n|\Big),
$$
et le th\'eor\`eme est d\'emontr\'e.\hfill\blacksquare

\bigskip
\state Remarque -- |  Dans la pratique, l'hypoth\`ese lipschitzienne faite 
sur $\Phi$ est tr\`es rarement satisfaite globalement pour
$y_1,y_2\in\bR$ et $h\in\bR$ (ne serait-ce que parce que le domaine de
d\'efinition de $\Phi$ est peut-\^etre plus petit). Par contre cette
hypoth\`ese est souvent satisfaite si on se restreint \`a des
$y_1,y_2\in J$ et $|h|\leq\delta$ o\`u $J$ est un intervalle ferm\'e
born\'e assez petit. Dans ce cas, la constante de stabilit\'e
$S=e^{\Lambda T}$ est valable pour des suites $y_n,\widetilde y_n\in
J$ et pour $h_{\max}\leq\delta$.

\bigskip
\state Exemples -- |  Supposons que la fonction $f$ soit lipschitzienne 
de rapport $k$ en $y$ et calculons $\Lambda,S$ pour les diff\'erentes
m\'ethodes d\'ej\`a pr\'esent\'ees.

\medskip
{\parindent1pc
\item{$\bullet$} Dans la m\'ethode d'Euler, $\Phi(t,y,h)=f(t,y)$.
On peut prendre $\Lambda=k$, $S=e^{kT}$.
\smallskip

\item{$\bullet$} Dans la m\'ethode du point milieu, on a
$$
\Phi(t,y,h)=f\Big(t+\frac{h}{2},y+\frac{h}{2}\ f(t,y)\Big).
$$
\par}

\noindent On en d\'eduit
$$
\eqalign{
&|\Phi(t,y_1,h)-\Phi(t,y_2,h)|\leq k\Big|y_1+\frac{h}{2}\ f(t,y_1)-\Big(y_2+\frac{h}{2}\
f(t,y_2)\Big)\Big|\cr
&\qquad\leq k\Big(|y_1-y_2|+\frac{h}{2}\ |f(t,y_1)-f(t,y_2)|\Big)\cr
&\qquad\leq k\Big(|y_1-y_2|+\frac{h}{2}\ k|y_1-y_2|\Big) = k\Big(1+\frac{1}{2}\ hk\Big)|y_1-y_2|.\cr 
}
$$
On peut prendre ici $\Lambda=k\Big(1+\frac{1}{2}\ h_{\max}\,k\Big)$, d'o\`u
$$
S=\exp\Big(kT\Big(1+\frac{1}{2}\ h_{\max}k\Big)\Big).
$$
Si $h_{\max}$ est petit (par rapport \`a $1/k^2 T$), cette constante est
du m\^eme ordre de grandeur que dans la m\'ethode d'Euler.

\bigskip
\majdeb
\state Corollaire -- |  {\it Lorsque $f$ est lipschitzienne en $y$, les 
m\'ethodes d'Euler et du point milieu sont convergentes.} 
\par\majfin\par

\bigskip

\noindent Le lecteur observera l'analogie des techniques utilis\'ees pour
obtenir ce corollaire avec celles utilis\'ees au chapitre V, \S$\,$3.1.

\subsection 2.4. Influence de l'ordre de la m\'ethode sur l'erreur globale|

Nous allons voir ici que l'ordre d'une m\'ethode num\'erique a une influence
d\'etermi\-nante sur la pr\'ecision que cette m\'ethode permet d'atteindre.
La d\'efinition de l'ordre que nous allons donne prend l'erreur de
consistance comme point de d\'epart (le lecteur prendra garde au d\'ecalage 
d'une unit\'e entre l'ordre et l'exposant de $h_n$ dans l'erreur de 
consistance~!$\,$)

\bigskip
\majdeb
\state D\'efinition -- |  {\it On dit qu'une m\'ethode \`a $1$ pas est 
d'ordre${}\geq p$ si pour toute solution exacte $z$ d'une \'equation
diff\'erentielle
$$
y'=f(t,y)\quad\hbox{o\`u $f$ est de classe $C^p$,}\leqno{\rm (E)}
$$
il existe une constante $C\geq 0$ telle que l'erreur de consistance
relative \`a $z$ v\'erifie
$$
|e_n|\leq Ch_n^{p+1},\quad\forall n,\quad 0\leq n<N.
$$
Elle est dite d'ordre $p$ $($exactement$)$ si elle est d'ordre${}\geq p$
mais pas d'ordre${}\geq p+1$.
}
\par\majfin\bigskip

Rappelons que l'erreur de consistance est donn\'ee par
$$
e_n=z(t_{n+1})-y_n-h_n\Phi(t_n,y_n,h_n)\quad\hbox{o\`u}\quad y_n=z(t_n).
$$
Supposons que $\Phi$ soit de classe $C^p$. La formule de Taylor donne alors
$$
\Phi(t_n,y_n,h_n)=\sum_{l=0}^p\frac{1}{l!}\ h^l_n\frac{\partial^l\Phi}{\partial h^l}\
(t_n,y_n,0)+o(h^p_n). 
$$
Si $f$ est de classe $C^p$, la solution $z$ est de classe $C^{p+1}$ donc 
$$
\eqalign{
z(t_{n+1})-y_n&=z(t_n+h)-z(t_n)\cr
&=\sum_{k=1}^{p+1}\frac{1}{k!}\ h^k_nz^{(k)}(t_n)+o(h^{n+1}_n)\cr
&=\sum_{l=0}^p \frac{1}{(l+1)!}\ h^{l+1}_nf^{[l]}(t_n,y_n)+o(h^{p+1}_n).\cr
}
$$
On en d\'eduit aussit\^ot
$$
e_n=\sum_{l=0}^p\frac{1}{l!}\ h^{l+1}_n\Big(\frac{1}{l+1}\
f^{[l]}(t_n,y_n)-\frac{\partial^l\Phi}{\partial h^l}\ \Big(t_n,y_n,0)\Big) +o(h_n^{p+1})
$$

\bigskip
\majdeb
\state Cons\'equence -- |  {\it La m\'ethode est d'ordre $\geq p$ si et 
seulement si $\Phi$ est telle que}
$$
\frac{\partial^l\Phi}{\partial h^l}\ (t,y,0)=\frac{1}{l+1}\ f^{[l]}(t,y),\quad 0\leq l\leq p-1.
$$ 
\vskip-\baselineskip
\majfin
\vskip\baselineskip

\noindent
Sous cette hypoth\`ese, l'erreur $e_n$ se r\'eduit \`a
$$
e^n=\frac{1}{p!}\ h^{p+1}_n\Big(\frac{1}{p+1}\ f^{[p]}(t_n,y_n)-\frac{\partial^p\Phi}{\partial h^p}\
(t_n,y_n,0)\Big)+o(h^{p+1}_n).
$$

\bigskip
\majdeb
\state Remarque -- |  {\it De ce qui pr\'ec\`ede on d\'eduit les 
\'equivalences 
$$
\hbox{M\'ethode consistante}\quad\Leftrightarrow\quad
\Phi(t,y,0)=f(t,y)\quad\Leftrightarrow\quad\hbox{M\'ethode d'ordre}\quad\geq 1.
$$}
\vskip-\baselineskip
\majfin
\vskip\baselineskip

L'utilisation de la formule de Taylor avec reste de Lagrange implique 
l'existence de points
$\tau_n\in\ ]t_n,t_{n+1}[$ et $\eta_n\in\ ]0,h_{\max}[$ tels que
$$
e_n=h_n^{p+1}\left(\frac{1}{(p+1)!}\ f^{[p]}(\tau_n,z(\tau_n))-\frac{1}{p!}\
\frac{\partial^p\Phi}{\partial h^p}\ (t_n,z(t_n),\eta_n)\right).
$$
Ceci permet (au moins th\'eoriquement) de trouver une constante $C$ dans la
majoration de l'erreur de consistance~: on peut prendre
$$
C=\frac{1}{(p+1)!}\ \Vert f^{[p]}(t,z(t))\Vert_{\infty}+\frac{1}{p!}\
\Vert\frac{\partial^p\Phi}{\partial h^p}\ (t,z(t),h)\Vert_{\infty} 
$$
o\`u les normes $\Vert\quad\Vert_{\infty}$ sont \'etendues aux\quad 
$(t,h)\in[t_0,t_0+T]\times[0,h_{\max}]$.

\bigskip
\state Majoration de l'erreur globale -- |  Compte tenu de la majoration
suppos\'ee satisfaite pour $e_n$, on a
$$
\sum_{0\leq n<N}|e_n|\leq\sum Ch^{p+1}_n\leq C\sum h_nh_{\max}^p\leq 
CT h^p_{\max}.
$$
Si la m\'ethode est stable avec constante de stabilit\'e $S$, on obtient 
donc la majoration

\bigskip
\majdeb
\hfil$\displaystyle\max_{0\leq n\leq N}|y_n-z(t_n)|\leq S(|y_0-z(t_0)|+CTh_{\max}^p)$.\hfill
\par\majfin\par

\bigskip
L'erreur initiale $|y_0-z(t_0)|$ est g\'en\'eralement n\'egligeable.
L'erreur globale donn\'ee par une m\'ethode stable d'ordre $p$ est
donc de l'ordre de grandeur de $h^p_{\max}$ avec une constante de
proportionnalit\'e SCT (on retiendra que l'ordre est \'egal \`a 
l'exposant de $h_{\max}$ dans la majoration de l'erreur globale,
alors que l'erreur de consistance, elle, est en $h_n^{p+1}$).

Si la constante SCT n'est pas trop grande (disons $\leq 10^2$), une
m\'ethode d'ordre 3 avec pas maximum $h_{\max}=10^{-2}$ permet
d'atteindre une pr\'ecision globale de l'ordre de $10^{-4}$.

\subsection 2.5. Influence des erreurs d'arrondi|

\noindent L'erreur globale calcul\'ee au \S\ 2.4 est une erreur th\'eorique,
c'est-\`a-dire qu'elle ne tient pas compte des erreurs d'arrondi qui
se produisent in\'evitablement en pratique. Dans la r\'ealit\'e
l'ordinateur va calculer non pas la suite r\'ecurrente $y_n$, mais une
valeur approch\'ee $\widetilde y_n$ de $y_n$ dans laquelle
interviendront

\medskip
\noindent$\bullet$ une erreur d'arrondi $\rho_n$ sur 
$\Phi(t_n,\widetilde y_n,h_n)$,

\noindent$\bullet$ une erreur d'arrondi $\sigma_n$ sur le calcul de
$\widetilde y_{n+1}$.

\medskip
En d\'efinitive, on aura
$$
\eqalign{
\widetilde y_{n+1}&=\widetilde y_n+h_n(\Phi(t_n,\widetilde y_n,h_n)+\rho_n)
+\sigma_n\cr
&=\widetilde y_n+h_n\Phi(t_n,\widetilde y_n,h_n)+h_{n}\rho_n+\sigma_n.\cr}
$$
Il se peut \'egalement que $\widetilde y_0$ diff\`ere
l\'eg\`erement de la valeur th\'eorique $y_0$~: $\widetilde
y_0=y_0+\varepsilon_0$.

\medskip
\majdeb
\state Hypoth\`ese -- | \quad $\forall n$,\quad $|\rho_n|\leq\rho$,\quad$|\sigma_n|\leq\sigma$.
\par\majfin\par

\medskip
Les constantes $\rho,\sigma$ d\'ependent des caract\'eristiques de l'ordinateur et de la pr\'ecision
des op\'erations arithm\'etiques (si les r\'eels sont cod\'es sur 6 octets, on a typiquement
$\rho=10^{-9}$, $\sigma=10^{-10}$, $|\varepsilon_0|\leq 10^{-10}$).

Si la m\'ethode est stable avec constante de stabilit\'e $S$, on en d\'eduit
$$
\eqalign{
\max_{0\leq n\leq N}|\widetilde y_n-y_n|&\leq S\Big(|\varepsilon_0|+\sum_{0\leq
n<N}(h_n|\rho_n|+|\sigma_n| )\Big)\cr
&\leq S(|\varepsilon_0|+T\rho+N\sigma).\cr
}$$
A cette erreur due aux arrondis s'ajoute l'erreur globale th\'eorique
$$
\max_{0\leq n\leq N}|y_n-z(t_n)|\leq SCT h^p_{\max}\qquad{\rm si}\qquad y_0=z(t_0).
$$
L'erreur totale commise est donc
$$
\max_{0\leq n\leq N}|\widetilde y_n-z(t_n)|\leq S(|\varepsilon_0|+T\rho+N\sigma+CT h^p_{\max})
$$
Supposons le pas $h_n=h$ constant pour simplifier. On a alors $N=\dfrac{T}{h}$ et l'erreur est
major\'ee par
$$
E(h)=S(|\varepsilon_0|+T\rho+\frac{T}{h}\,\sigma+CT h^p)
=S(|\varepsilon_0|+T\rho)+ST\Big(\frac{\sigma}{h}+Ch^p\Big)
$$
L'\'etude de $E(h)$ donne la courbe suivante, avec un minimum
de $E(h)$ r\'ealis\'e en 
$h_{\rm opt}=\big(\dfrac{\sigma}{pC}\big)^{\fracss{1}{p+1}}$. 

\InsertFig -37.000 54.000 {
1 mm unit
  0.000 -47.000 translate
 60.000  55.000 moveto 
 80.000   0.000 vector 
 70.000  50.000 moveto 
 45.000  90.000 vector 
 71.990  95.246 moveto 
[  71.990  95.246    72.028  94.497    72.067  93.751    72.107  93.006   
 72.150  92.263    72.193  91.521    72.239  90.782    72.286  90.045   
 72.335  89.309    72.387  88.577    72.440  87.846    72.496  87.118   
 72.554  86.393    72.614  85.670    72.677  84.951    72.744  84.235   
 72.813  83.522    72.885  82.812    72.961  82.107    73.040  81.405   
 73.124  80.708    73.211  80.015    73.303  79.326    73.400  78.643   
 73.502  77.966    73.610  77.294    73.723  76.628    73.842  75.969   
 73.968  75.317    74.102  74.672    74.243  74.036    74.392  73.408   
 74.550  72.789    74.718  72.180    74.896  71.581    75.084  70.993   
 75.285  70.417    75.498  69.854    75.724  69.305    75.964  68.769   
 76.220  68.249    76.492  67.745    76.781  67.258    77.088  66.789   
 77.414  66.340    77.761  65.910    78.129  65.502    78.520  65.116   
 78.934  64.753    79.373  64.414    79.837  64.100    80.328  63.813   
 80.846  63.552    81.393  63.320    81.969  63.116    82.576  62.943   
 83.213  62.800    83.884  62.690    84.587  62.612    85.325  62.569   
 86.098  62.561    86.888  62.589    87.645  62.647    88.375  62.730   
 89.082  62.834    89.771  62.955    90.442  63.092    91.100  63.243   
 91.744  63.407    92.376  63.582    92.997  63.767    93.608  63.962   
 94.210  64.166    94.804  64.378    95.390  64.597    95.968  64.824   
 96.539  65.058    97.103  65.298    97.662  65.544    98.214  65.796   
 98.760  66.053    99.301  66.316    99.837  66.584   100.368  66.857   
100.894  67.134   101.416  67.416   101.933  67.702   102.446  67.993   
102.954  68.287   103.459  68.586   103.959  68.888   104.456  69.194   
104.950  69.504   105.439  69.817   105.926  70.133   106.409  70.453   
106.888  70.776   107.365  71.102   107.838  71.432   108.308  71.764   
108.776  72.099   109.240  72.437   109.702  72.778   110.161  73.121   
110.617  73.467   111.071  73.816   111.522  74.167   111.970  74.521   
112.416  74.877   112.859  75.236   113.301  75.597   113.739  75.960   
114.176  76.326   114.610  76.694   115.042  77.064   115.472  77.436   
115.899  77.810   116.325  78.187   116.748  78.565   117.170  78.946   
117.589  79.328   118.007  79.713   118.422  80.099   118.836  80.487   
119.247  80.877   119.657  81.269   120.065  81.663   120.471  82.059   
120.876  82.456   121.279  82.855   121.680  83.256   122.079  83.658   
122.476  84.062   122.872  84.468   123.267  84.875   123.659  85.284   
124.050  85.695   124.440  86.107   124.828  86.520   125.215  86.935   
125.599  87.352   125.983  87.770   126.365  88.189   126.746  88.610   
127.125  89.033   127.502  89.457   127.879  89.882   128.254  90.308   
128.627  90.736   128.999  91.165   129.370  91.596   129.740  92.028   
130.108  92.461   130.475  92.895   ] scurve 
 85.874  55.000 moveto 
[   2.000   1.000 ] 0 setdash 
 85.874  56.500 moveto 
 85.874  62.500 lineto 
stroke 
 71.500  62.560 moveto 
 85.800  62.560 lineto 
stroke 
 85.874  55.000 moveto 
  0.340 setlinewidth 
  0.700  90.000 indent 
 85.874  62.560 moveto 
  0.700 disk 
 70.000  62.560 moveto 
  0.700   0.000 indent}
\LabelTeX   61.000  46.000 $E(h)$\ELTX
\LabelTeX   83.500   4.000 $h_{\rm opt}$\ELTX
\LabelTeX  139.000   4.000 $h$\ELTX
\EndFig

Typiquement, pour une m\'ethode d'ordre $p=2$ o\`u $pC\simeq 10$, on
obtient $h_{\rm opt}\simeq 10^{-3}$. Si on prend un pas plus petit,
l'erreur augmente~!  Ceci est d\^u au fait que le nombre de pas
$N=\dfrac{T}{h}$ augmente, et avec lui les erreurs d'arrondi, lorsque
le pas $h$ diminue. Les erreurs d'arrondi l'emportent alors sur
l'erreur globale th\'eorique $SCT h^p$. L'exp\'erience num\'erique
suivante confirme ces pr\'evisions th\'eoriques.

\bigskip
\state Exemple -- |  Consid\'erons le probl\`eme  de Cauchy $y'=y$ avec
donn\'ee $y_0=1$ en $t_0=0$. La solution exacte est $y(t)=e^t$, d'o\`u
$y(1)=e\simeq 2,7182818285$. Si l'on utilise la m\'ethode du point
milieu avec pas constant $h$, on obtient l'algorithme
$$
y_{n+1}=(1+h+h^2/2)y_n,\quad y_0=1.
$$
L'erreur de consistance est donn\'ee par $e_n\sim h^3 y_n/6$, d'o\`u
$e_n\leq C h^3$ avec $C=e/6$ sur $[0,T]=[0,1]$.
Par ailleurs, on peut au mieux esp\'erer $\sigma=10^{-11}$, ce qui donne
$$
h_{\rm opt}\geq\left(\frac{10^{-11}}{2\cdot e/6}\right)^{1/3}\simeq 
2,224\cdot 10^{-4},\qquad
N_{\rm opt}=\frac{T}{h_{\rm opt}}<4500.
$$
Un calcul en Turbo-Pascal donne en fait les r\'esultats suivants~:
$$
\vbox{\offinterlineskip\halign{
\quad#\hfill\quad&
\tv\quad\hfill#\hfill\quad&
\tv\quad\hfill#\hfill\quad\cr
\hfill Nombre de pas&Valeur $y_N$ associ\'ee&Erreur $y(1)-y_N$\cr
\noalign{\hrule}
$N=\quad\enspace 10$&2,7140808465&$4,2\cdot 10^{-3}$\cr
$N=\quad 100$&2,7182368616&$4,5\cdot 10^{-5}$\cr
$N=\quad 500$&2,7182800146&$4,8\cdot 10^{-6}$\cr
$N=\enspace 1000$&2,7182813650&$4,6\cdot 10^{-7}$\cr
$N=\enspace 2000$&2,7182816975&$1,3\cdot 10^{-7}$\cr
$N=\enspace 3000$&2,7182817436&$8,5\cdot 10^{-8}$\cr
$N=\enspace 4000$&2,7182817661&$6,2\cdot 10^{-8}$\cr
$N=\enspace 4400$&2,7182817882&$4,0\cdot 10^{-8}$\cr
$N=\enspace 5000$&2,7182817787&$5,0\cdot 10^{-8}$\cr
$N=\enspace 6000$&2,7182817607&$6,8\cdot 10^{-8}$\cr
$N=\enspace 7000$&2,7182817507&$7,8\cdot 10^{-8}$\cr
$N=10000$&2,7182817473&$8,1\cdot 10^{-8}$\cr
$N=20000$&2,7182817014&$1,3\cdot 10^{-7}$\cr
}}
$$
Le nombre de pas optimal observ\'e est $N_{\rm opt}\simeq 4400$.

\subsection 2.6. Probl\`emes bien pos\'es, bien conditionn\'es, 
probl\`emes raides|

\noindent 
L'objet de ce paragraphe est de mettre en \'evidence les difficult\'es
qui peuvent appara\^\i tre dans la mise en {\oe}uvre des algorithmes de
r\'esolution num\'erique.

\bigskip
\majdeb
\state D\'efinition 1 -- |  {\it On dit qu'un probl\`eme de Cauchy est
math\'ematiquement bien pos\'e si la solution est unique et d\'epend
contin\^ument de la donn\'ee initiale.}
\par\majfin\par

\bigskip
\state Exemple -- |  Consid\'erons le probl\`eme de Cauchy
$$
\cases{
y'=2\sqrt{|y|},\qquad t\in[0,+\infty[\cr
\noalign{\vskip3pt}
y(0)=0.\hfill\cr
}
$$
Ce probl\`eme admet les solutions $y(t)=0$, $y(t)=t^2$ et plus g\'en\'eralement
$$
\cases{
y(t)=0,\hfill&$t\in[0,a]$\hfill\cr
\noalign{\vskip3pt}
y(t)=(t-a)^2,&$t\in[a,+\infty[$.\hfill\cr
}
$$
L'utilisation de la m\'ethode d'Euler 
$y_{n+1}=y_n+2h_n\sqrt{y_n}$ va conduire aux solutions
approch\'ees suivantes~:

\medskip
\noindent$\bullet$ si $y_0=0$\qquad $y(t)=0$

\noindent$\bullet$ si $y_0=\varepsilon$\qquad$y(t)\simeq (t+\sqrt{\varepsilon})^2$ quand $h_{\rm
max}\rightarrow 0$.

\medskip
Il n'y a ici ni unicit\'e, ni continuit\'e de la solution. Le probl\`eme de Cauchy est donc
math\'ematiquement mal pos\'e.

Les r\'esultats du chapitre V \S\ 3.2 (et ceux \`a venir du chapitre
XI \S\ 1.2) montrent que le probl\`eme de Cauchy est {\it
math\'ematiquement bien pos\'e d\`es que $f(t,y)$ est localement
lipschitzienne} en $y$.

\bigskip 
\majdeb 
\state D\'efinition 2 -- |  {\it On dit qu'un
probl\`eme de Cauchy est num\'eriquement bien pos\'e si la
continuit\'e de la solution par rapport \`a la donn\'ee initiale est
suffisamment bonne pour que la solution ne soit pas perturb\'ee par
une erreur initiale ou des erreurs d'arrondi faibles.}
\par\majfin\par

\bigskip
\noindent 
Par l'expression \og continuit\'e suffisamment bonne\fg, on entend en
g\'en\'eral l'existence d'une constante de Lipschitz petite en regard
de la pr\'ecision des calculs. On notera que la d\'efinition 2 ne fait
pas r\'ef\'erence \`a la m\'ethode de calcul utilis\'ee.

\bigskip
\state Exemple 2 -- |  Soit le probl\`eme de Cauchy
$$
\cases{
y'=3y-1,\qquad t\in [0,10]\cr
\noalign{\vskip3pt}
y(0)=\dfrac{1}{3}.\hfill\cr
}
$$
La solution exacte est $y(t)=t+\dfrac{1}{3}$. La donn\'ee initiale $\widetilde
y(0)=\dfrac{1}{3}+\varepsilon$ fournit $\widetilde y(t)=t+\dfrac{1}{3}+\varepsilon e^{3t}$. On a
alors 
$$
\widetilde y{(10)}-y{(10)}=\varepsilon\cdot e^{30}\simeq 10^{13}\varepsilon.
$$
Le probl\`eme est ici math\'ematiquement bien pos\'e, mais num\'eriquement mal pos\'e si la
pr\'ecision des calculs est seulement de $10^{-10}$. Le probl\`eme redevient num\'eriquement bien
pos\'e si la pr\'ecision des calculs est de $10^{-20}$.

\bigskip
\state Exemple 3 -- |  L'exemple suivant montre que m\^eme un probl\`eme num\'eriquement bien
pos\'e peut soulever des difficult\'es inattendues~:


$$
\cases{
y'=-150y+30,\quad  t\in[0,1],\cr
\noalign{\vskip3pt}
y(0)=\dfrac{1}{5}.\hfill\cr
}
$$
La solution exacte est $y(t)=\dfrac{1}{5}$ et la donn\'ee initiale $\widetilde
y(0)=\dfrac{1}{5}+\varepsilon$ fournit $\widetilde y(t)=\dfrac{1}{5}+\varepsilon e^{-150t}$. Comme
$0\leq e^{-150 t}\leq 1$ sur $[0,1]$, le probl\`eme est num\'eriquement bien pos\'e. La m\'ethode
d'Euler avec pas constant $h$ donne 
$$
\leqalignno{
&y_{n+1}=y_n+h(-150 y_n+30)=(1-150 h)y_n+30 h,\cr
&y_{n+1}-\frac{1}{5}=(1-150h)(y_n-\frac{1}{5}),\cr
&y_n-\frac{1}{5}=(1-150h)^n\Big(y_0-\frac{1}{5}\Big).&\hbox{d'o\`u}\cr
}
$$
Supposons $h=\dfrac{1}{50}$. Une erreur initiale $y_0=\dfrac{1}{5}+\varepsilon$ conduit \`a
$y_n=\dfrac{1}{5}+(-2)^m\varepsilon$, d'o\`u pour $t=1$
$$
y_{50}=\frac{1}{5}+2^{50}\varepsilon\simeq\frac{1}{50}+10^{15}\varepsilon!
$$
Pour que $|y_n|$ ne diverge pas vers $+\infty$, il est n\'ecessaire de prendre $|1-150h|\leq 1$, soit
$150h\leq 2$, $h\leq \dfrac{1}{75}$. Bien que le probl\`eme soit tout \`a fait bien pos\'e, on voit
qu'il est n\'ecessaire de prendre un pas assez petit, et donc de faire des calculs plus co\^uteux
que d'ordinaire.

\bigskip
\majdeb
\state D\'efinition 3 -- |  {\it On dit qu'un probl\`eme est bien conditionn\'e si les
m\'ethodes num\'eriques usuelles peuvent en donner la solution en un temps raisonnable.}
\par\majfin\par

\bigskip
Un probl\`eme sera bien conditionn\'e si la constante de stabilit\'e $S$ n'est pas trop grande
(disons nettement $<10^{10}$  si la pr\'ecision des
calculs est $10^{-10}$). Sinon, on dit qu'on a affaire \`a un {\it  probl\`eme raide}.

On sait qu'en g\'en\'eral la constante de stabilit\'e $S$ est
major\'ee par $e^{\Lambda T}$. Dans un probl\`eme raide, on peut avoir
typiquement $\Lambda T=10^3$, $e^{\Lambda T}>10^{400}$. Il existe des
algorithmes permettant de traiter certains probl\`emes raides, mais
nous n'aborderons pas cette question. Le lecteur pourra consulter sur
ce point le livre de Crouzeix-Mignot.
\eject

\section 3. M\'ethodes de Runge-Kutta|

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

\noindent On consid\`ere un probl\`eme de Cauchy
$$
\cases{
y'=f(t,y),\qquad t\in[t_0,t_0+T]\cr
\noalign{\vskip3pt}
y(t_0)=y_0\hfill\cr
}
$$
et on cherche \`a discr\'etiser ce probl\`eme par rapport \`a une subdivision
$t_0<t_1<\dots<t_N=t_0+T$. L'id\'ee est de calculer par r\'ecurrence les points $(t_n,y_n)$ en
utilisant des points interm\'ediaires $(t_{n,i},y_{n,i})$ avec 
$$
t_{n,i}=t_n+c_i h_n,\quad 1\leq i\leq q,\quad c_i\in[0,1].
$$
A chacun de ces points on associe la pente correspondante
$$
p_{n,i}=f(t_{n,i},y_{n,i}).
$$
Soit $z$ une solution exacte de l'\'equation. On a
$$
\eqalign{
z(t_{n,i})&=z(t_n)+\int^{t_{n,i}}_{t_n} f(t,z(t))dt\cr
&=z(t_n)+h_n\int^{c_i}_0 f(t_n+uh_n, z(t_n+uh_n))du\cr
}
$$
gr\^ace au changement de variable $t=t_n+uh_n$. De m\^eme
$$
z(t_{n+1})=z(t_n)+h_n\int^1_0f(t_n+uh_n,z(t_n+uh_n)du.
$$
On se donne alors pour chaque $i=1,2,\dots,q$ une m\'ethode d'int\'egration approch\'ee
$$
\int^{c_i}_0g(t)dt\simeq\sum_{1\leq j<i}a_{ij}g(c_j),\leqno({\rm M}_i)
$$
ces m\'ethodes pouvant \^etre {\it a priori} diff\'erentes. On se donne \'egalement une m\'ethode
d'int\'egration approch\'ee sur $[0,1]$~:
$$
\int^1_0g(t)dt\simeq\sum_{1\leq j\leq q}b_jg(c_j).\leqno({\rm M})
$$
En appliquant ces m\'ethodes d'int\'egration \`a  $g(u)=f(t_n+uh_n,z(t_n+uh_n))$, il
vient
$$
\eqalign{
z(t_{n,i})&\simeq z(t_n)+h_n\sum_{1\leq j<i}a_{ij}f(t_{n,j},z(t_{n,j})),\cr
z(t_{n+1})&\simeq z(t_n)+h_n\sum_{1\leq j\leq q}b_{j}f(t_{n,j},z(t_{n,j})).\cr
}
$$
La m\'ethode Runge-Kutta correspondante est d\'efinie par l'algorithme
$$
\cases{
\left[\matrix{
t_{n,i}=t_n+c_ih_n\hfill\cr
\noalign{\vskip3pt}
\displaystyle y_{n,i}=y_n+h_n\sum_{1\leq j< i}a_{ij}p_{n,j}\hfill\cr
\noalign{\vskip3pt}
p_{n,i}=f(t_{n,i},y_{n,i})\hfill\cr
}\right]\qquad 1\leq i\leq q\cr
\noalign{\vskip6pt}
\enspace t_{n+1}=t_n+h_n\hfill\cr
\noalign{\vskip3pt}
\enspace\displaystyle y_{n+1}=y_n+h_n\sum_{1\leq j\leq q}b_j p_{n,j}.\hfill\cr
}
$$
On la repr\'esente conventionnellement par le tableau

$$
\vbox{\offinterlineskip\halign{
#\hfill\quad&
\hfill\quad#\hfill\quad&
\tv\quad\hfill#\hfill\quad&
\hfill#\quad\hfill&
\hfill#\quad\hfill&
\hfill#\quad\hfill&
\hfill#\quad\hfill\cr
(M$_1$)&$c_1$&0&0&$\dots$&0&0\cr
(M$_2$)&$c_2$&$a_{21}$&0&$\dots$&0&0\cr
&$\vdots$&$\vdots$&$\vdots$&$\ddots$&$\vdots$&$\vdots$\cr
\noalign{\vskip-8pt}
&$\vdots$&$\vdots$&$\vdots$&&0&0\cr
(M$_q$)&$c_q$&$a_{q1}$&$a_{q2}$&$\dots$&$a_{qq-1}$&0\cr
&\multispan6\hrulefill\cr
(M)&&$b_1$&$b_2$&$\dots$&$b_{q-1}$&$b_q$\cr
}}
$$

\vskip6pt
\noindent o\`u les m\'ethodes d'int\'egration approch\'ees correspondent aux lignes. On pose par
convention $a_{ij}=0$ pour $j\geq i$.


\bigskip
\majdeb
\state Hypoth\`ese -- |  {\it On supposera toujours que les m\'ethodes d'int\'egration (M$_i$)
et (M) sont d'ordre 0 au moins, c'est-\`a-dire
$$
c_i=\sum_{1\leq j<i}a_{ij},\qquad 1=\sum_{1\leq j\leq q}b_j.
$$} 
\vskip-9pt
\par\majfin\par

\bigskip
En particulier, on aura toujours
$$
c_1=0,\quad t_{n,1}=t_n,\quad y_{n,1}=y_n,\quad p_{n,1}=f(t_n,y_n).
$$

\subsection 3.2. Exemples|

\state Exemple 1 -- |  Pour  $q=1$, le seul choix possible est\qquad$\matrix{
0&\Big|&0\cr
\noalign{\vskip-6pt}
\noalign{\hbox to 3pc{\hrulefill}}
\noalign{\vskip-3pt}
&\Big|&1\cr
}$

On a ici $c_1=0$, $a_{11}=0$, $b_1=1$. L'algorithme est donn\'e par
$$
\cases{
p_{n,1}=f(t_n,y_n)\hfill\cr
\noalign{\vskip3pt}
t_{n+1}=t_n+h_n\hfill\cr
\noalign{\vskip3pt}
y_{n+1}=y_n+h_np_{n,1}\hfill\cr
}
$$
Il s'agit de la m\'ethode d'Euler.

\bigskip
\state Exemple 2 -- |  Pour $q=2$, on consid\`ere les tableaux de la forme
$$
\matrix{
0&\Big|&0&0\cr
\noalign{\vskip-2pt}
\alpha&\Big|&\alpha&0\cr
\noalign{\vskip-6pt}
\noalign{\hbox to 7pc{\hrulefill}}
\noalign{\vskip-2pt}
&\Bigg|&\displaystyle 1-\dfrac{1}{2\alpha}&\displaystyle\dfrac{1}{2\alpha}\cr
},\qquad\qquad\hbox{o\`u}\quad \alpha\in\ ]0,1].
$$
L'algorithme s'\'ecrit ici
$$
\cases{
p_{n,1}=f(t_n,y_n)\hfill\cr
\noalign{\vskip3pt}
t_{n,2}=t_n+\alpha h_n\hfill\cr
\noalign{\vskip3pt}
y_{n,2}=y_n+\alpha h_np_{n,1}\hfill\cr
\noalign{\vskip3pt}
p_{n,2}=f(t_{n,2},y_{n,2})\hfill\cr
\noalign{\vskip3pt}
t_{n+1}=t_n+h_n\hfill\cr
\noalign{\vskip3pt}
y_{n+1}=y_n+h_n\Big(\Big(1-\dfrac{1}{2\alpha}\Big)\ p_{n,1}+\dfrac{1}{2\alpha}\
p_{n,2}\Big),\hfill\cr 
}
$$
ou encore, sous forme condens\'ee~:
$$
y_{n+1}=y_n+h_n\Big(\Big(1-\dfrac{1}{2\alpha}\Big)\ f(t_n,y_n)+\dfrac{1}{2\alpha}\ f(t_n+\alpha
h_n,y_n+\alpha h_nf(t,y_n))\Big). 
$$
C'est n\'eanmoins la premi\`ere formulation qui est la plus efficace en pratique, puisqu'elle
requiert seulement deux \'evaluations de la fonction $f$ au lieu de 3 pour la forme condens\'ee. 

\medskip
{\parindent1pc
\item{$\bullet$} Pour $\alpha=\dfrac{1}{2}$, on retrouve la m\'ethode du point milieu
$$
y_{n+1}=y_n+h_nf\Big(t_n+\frac{h_n}{2},y_n+\frac{h_n}{2}\ f(t_n,y_n)\Big),
$$
\item{} qui est bas\'ee sur la m\'ethode d'int\'egration du point milieu~:
$$
\int^1_0 g(t)dt\simeq g\Big(\frac{1}{2}\Big).\leqno{\rm (M)}
$$

\item{$\bullet$} Pour $\alpha=1$, on obtient la {\it m\'ethode de Heun}~:
$$
y_{n+1}=y_n+h_n\Big(\frac{1}{2}\ f(t_n,y_n)+\frac{1}{2}\ f(t_{n+1},y_n+h_nf(t_n,y_n))\Big),
$$
\item{} qui repose sur la m\'ethode d'int\'egration des trap\`ezes~:
$$
\int^1_0 g(t)dt\simeq \frac{1}{2}\ (g(0)+g(1)).\leqno{\rm (M)}
$$}
\smallskip

\vbox{\state Exemple 3 -- |  {\it M\'ethode de Runge-Kutta \og classique\fg~:}

Il s'agit de la m\'ethode d\'efinie par le tableau}
$$
q=4,\qquad\qquad\matrix{
0&\Big|&0&0&0&0\cr
\noalign{\vskip-2pt}
\dfrac{1}{2}&\Bigg|&\dfrac{1}{2}&0&0&0\cr
\noalign{\vskip-2pt}
\dfrac{1}{2}&\Big|&0&\dfrac{1}{2}&0&0\cr
\noalign{\vskip-2pt}
1&\Big|&0&0&1&0\cr
\noalign{\vskip-6pt}
\noalign{\hbox to 9pc{\hrulefill}}
\noalign{\vskip-2pt}
&\Bigg|&\dfrac{1}{6}&\dfrac{2}{6}&\dfrac{2}{6}&\dfrac{1}{6}\cr
}
$$
L'algorithme correspondant s'\'ecrit
$$
\cases{
p_{n,1}=f(t_n,y_n)\hfill\cr
\noalign{\vskip3pt}
t_{n,2}=t_n+\dfrac{1}{2}\ h_n\hfill\cr
\noalign{\vskip3pt}
y_{n,2}=y_n+\dfrac{1}{2}\ h_np_{n,1}\hfill\cr
\noalign{\vskip3pt}
p_{n,2}=f(t_{n,2},y_{n,2})\hfill\cr
\noalign{\vskip3pt}
y_{n,3}=y_n+\dfrac{1}{2}\ h_np_{n,2}\hfill\cr
\noalign{\vskip3pt}
p_{n,3}=f(t_{n,2},y_{n,3})\qquad\quad\hbox{(noter que $t_{n,3}=t_{n,2}$)}\hfill\cr
\noalign{\vskip3pt}
t_{n+1}=t_n+h_n\qquad\qquad\kern6pt\hbox{(noter que $t_{n,4}=t_{n+1}$)}\hfill\cr
\noalign{\vskip3pt}
y_{n,4}=y_n+h_np_{n,3}\hfill\cr
\noalign{\vskip3pt}
p_{n,4}=f(t_{n+1}, y_{n,4})\hfill\cr
\noalign{\vskip3pt}
y_{n+1}=y_n+h_n\Big(\dfrac{1}{6}\ p_{n,1}+\dfrac{2}{6}\ p_{n,2}+\dfrac{2}{6}\ p_{n,3}+\dfrac{1}{6}\
p_{n,4}\Big)\hfill\cr 
}
$$

\noindent On verra plus loin que cette m\'ethode est d'ordre 4.
Dans ce cas les m\'ethodes d'int\'egration (M$_i$) et (M) utilis\'ees sont respectivement~:
$$
\matrix{
({\rm M}_2)&\displaystyle\int^{\frac{1}{2}}_0 g(t)dt\simeq \dfrac{1}{2}\
g(0):\qquad\hbox{rectangles \`a gauche,}\hfill\cr 
\noalign{\vskip3pt}
({\rm M}_3)&\displaystyle\int^{\frac{1}{2}}_0 g(t)dt\simeq \dfrac{1}{2}\
g\Big(\dfrac{1}{2}\Big):\quad\kern4pt\hbox{rectangles \`a droite,}\hfill\cr 
\noalign{\vskip3pt}
({\rm M}_4)&\displaystyle\int^1_0g(t)dt\simeq g\Big(\dfrac{1}{2}\Big):\qquad\kern8pt\hbox{point
milieu,}\hfill\cr 
\noalign{\vskip3pt}
({\rm M})\hfill&\displaystyle\int^1_0 g(t)dt\simeq \dfrac{1}{6}\ g(0)+\dfrac{2}{6}\
g\Big(\dfrac{1}{2}\Big)+ \dfrac{2}{6}\ g\Big(\dfrac{1}{2}\Big)+\dfrac{1}{6}\
g(1):\quad\hbox{Simpson.}\hfill\cr  }
$$
\vskip-3pt

\subsection 3.3. Stabilit\'e des m\'ethodes de Runge-Kutta|

\noindent Les m\'ethodes de Runge-Kutta sont des m\'ethodes \`a un pas
$$
y_{n+1}=y_n+h_n\Phi(t_n,y_n,h_n)
$$
avec $\displaystyle\Phi(t_n,y_n,h_n)=\sum_{1\leq j\leq q}b_jp_{n,j}$. La fonction $\Phi$ est
d\'efinie de mani\`ere explicite par
$$
\cases{
\displaystyle \Phi(t,y,h)=\sum_{1\leq j\leq q}b_jf(t+c_jh,y_j)\quad{\rm avec}\hfill\cr
\noalign{\vskip3pt}
\displaystyle y_i=y+h\sum_{1\leq j< i}a_{ij}f(t+c_jh,y_j),\quad 1\leq i\leq q.\hfill\cr
}\eqno(*)
$$
Supposons que $f$ soit $k$-lipschitzienne en $y$. On va montrer que $\Phi$ est alors \'egalement
lipschitzienne. Soit $z\in\bR$ et supposons $\Phi(t,z,h)$ et $z_i$ d\'efinis \`a partir de $z$ comme
dans la formule $(*)$.


\bigskip
\majdeb
\state Lemme -- |  {\it Soit $\displaystyle\alpha=\max_i\Big(\sum_{1\leq j\leq
i}|a_{ij}|\Big)$. Alors
$$
|y_i-z_i|\leq (1+(\alpha kh)+(\alpha kh)^2+\dots+(\alpha kh)^{i-1})|y-z|.
$$
\vskip-9pt
} 
\par\majfin\par
\bigskip

On d\'emontre le lemme par r\'ecurrence sur $i$. Pour $i=1$, on a $y_1=y$, $z_1=z$ et le r\'esultat
est \'evident. Supposons l'in\'egalit\'e vraie pour tout $j<i$. Alors
$$
\eqalign{
|y_i-z_i|&\leq |y-z|+h\sum_{j<i}|a_{ij}|\cdot k\cdot \max_{j<i}|y_j-z_j|,\cr
|y_i-z_i|&\leq |y-z|+\alpha kh\max_{j<i}|y_j-z_j|.\cr
}$$
Par hypoth\`ese de r\'ecurrence il vient
$$
\max_{j<i}|y_j-z_j|\leq(1+\alpha kh+\dots+(\alpha kh)^{i-2})|y-z|,
$$
et l'in\'egalit\'e s'ensuit \`a l'ordre $i$.\hfill\blacksquare

\medskip
La formule $(*)$ entra\^\i ne maintenant
$$
|\Phi(t,y,h)-\Phi(t,z,h)|\leq\sum_{1\leq j\leq q}|b_j|\ k\ |y_j-z_j|
\leq\Lambda|y-z|\qquad\hbox{\rm avec}
$$

\bigskip
\majdeb
\hfil$\displaystyle\Lambda=k\sum_{1\leq j\leq q}|b_j|(1+(\alpha kh_{\max})
+\dots+(\alpha kh_{\rm max})^{j-1})$.\hfill 
\par\majfin\par

\bigskip
\majdeb
\state Corollaire -- |  {\it Les m\'ethodes de Runge-Kutta sont stables, avec 
constante de stabilit\'e $S=e^{\Lambda T}$.} 
\par\majfin\par

\bigskip
\state Remarque -- |  Dans le cas fr\'equent o\`u les coefficients $b_j$ 
sont $\geq 0$, on a la relation 
$$
\Lambda\leq k(1+(\alpha kh_{\max})+\dots+(\alpha kh_{\max})^{q-1}).
$$
Si les coefficients $a_{ij}$ sont eux-m\^emes $\geq 0$, on a 
$\displaystyle\alpha=\max_i c_i$.

Lorsque $h_{\rm \max}$ est assez petit devant $1/\alpha k$, la
constante de stabilit\'e est donc de l'ordre de grandeur de $e^{kT}$.
Ces observations montrent que les m\'ethodes de Runge-Kutta d\'ecrites
dans les exemples 1, 2, 3 du \S$\,$3.3 poss\`edent une excellente
stabilit\'e (il est facile de voir que $e^{kT}$ est la borne
inf\'erieure possible pour $S$, quelle que soit la m\'ethode
utilis\'ee~: consid\'erer pour cela l'\'equation $y'=ky)$.

\subsection 3.4. Ordre des m\'ethodes de Runge-Kutta|

\noindent Pour d\'eterminer l'ordre, on peut appliquer le crit\`ere du 
\S\ 2.4 consistant \`a \'evaluer les d\'eriv\'ees
$\dfrac{\partial^l\Phi}{\partial h^l}\ (t,y,0)$~: l'ordre est au moins
\'egal \`a $p$ si et seulement si cette d\'eriv\'ee est \'egale \`a
$\dfrac{1}{l+1}\ f^{[l]}(t,y)$ pour $l\leq p-1$. Gr\^ace \`a la
formule $(*)$ du \S\ 3.3, on obtient facilement les d\'eriv\'ees
successives de $\Phi$~:

\medskip

\noindent$\bullet$ $\displaystyle\Phi(t,y,0)=\sum_{1\leq j\leq q}b_jf(t,y)=f(t,y)$.

Les m\'ethodes de Runge-Kutta sont donc toujours d'ordre $\geq 1$ (c'est-\`a-dire
consistantes).

\medskip
\noindent$\bullet$\kern7pt $\displaystyle\dfrac{\partial\Phi}{\partial h}\ (t,y,h)=\sum_j
b_j\Big(c_jf'_t(t+c_jh,y_j)+f'_y(t+c_jh,y_j)\ \frac{\partial y_j}{\partial h}\Big)$,

$\displaystyle\dfrac{\partial y_i}{\partial
h}=\sum_{j<i}a_{ij}f(t+c_jh,y_j)+h\sum_{j<i}a_{ij}\Big(c_jf'_t+f'_y\ \frac{\partial y_j}{\partial
h}\Big)$.

\noindent Pour $h=0$, on obtient donc
$$
\eqalign{
&\frac{\partial y_i}{\partial h}\Big|_{h=0}=\Big(\sum_{j<i}a_{ij}\Big)f(t,y)=c_if(t,y)\cr
&\frac{\partial\Phi}{\partial h}\ (t,y,0)=\sum_j b_jc_j(f'_t+f'_yf)(t,y)=\Big(\sum b_jc_j\Big)
f^{[1]}(t,y).\cr 
}
$$
D'apr\`es le \S\ 2.4, la m\'ethode est d'ordre $\geq 2$ si et seulement si 
$\sum b_jc_j=\dfrac{1}{2}$.

\medskip
\noindent$\bullet$\kern7pt $\displaystyle\dfrac{\partial^2\Phi}{\partial h^2}
\ (t,y,h)=\sum_j b_j\Big(c^2_jf''_{tt}+2c_jf''_{ty}\ 
\dfrac{\partial y_j}{\partial h}+f''_{yy}\ \Big(\dfrac{\partial
y_j}{\partial h}\Big)^2+f'_y\ \dfrac{\partial^2 y_j}{\partial h^2}\Big)$,

$\displaystyle\dfrac{\partial^2 y_i}{\partial h^2}
=2\sum_{j<i}a_{ij}\Big(c_jf'_t+f'_y\ \dfrac{\partial y_j}{\partial h}\Big)
+h\sum_{j<i}a_{ij}\Big(c^2_jf''_{tt}+\dots\Big)$.

\noindent Pour $h=0$, il vient
$$
\eqalign{
&\frac{\partial^2 y_i}{\partial h^2}\Big|_{h=0}
=2\sum a_{ij}c_j(f'_t+f'_yf)(t,y),\cr
&\frac{\partial^2\Phi}{\partial h^2}\ (t,y,0)
=\sum_j b_jc^2_j(f''_{tt}+2f''_{ty}f+f''_{yy}f^2)(t,y)
+2\sum_{i,j}b_ia_{ij}c_jf'_y(f'_t+f'_yf)(t,y).\cr}
$$
Or $f^{[2]}$ est donn\'e par
$$
\eqalign{
f^{[2]}(t,y)&=(f^{[1]})'_t+(f^{[1]})'_y f\cr
&=(f'_t+f'_yf)'_t+(f'_t+f'_yf)'_yf\cr
&=f''_{tt}+f''_{ty}f+f'_yf'_t+f''_{ty}f+f''_{yy}f^2+f'^2_y f.\cr
&=(f''_{tt}+2f''_{ty}f+f''_{yy}f^2)+f'_y(f'_t+f'_yf).\cr
}
$$
La condition $\dfrac{\partial^2\Phi}{\partial h^2}\ (t,y,0)
=\dfrac{1}{3}\ f^{[2]}(t,y)$ se traduit en g\'en\'eral par les conditions
$$
\sum_jb_jc^2_j=\dfrac{1}{3},\qquad\sum_{i,j}b_ia_{ij}c_j=\dfrac{1}{6}
$$
(prendre respectivement $f(t,y)=t^2$, puis $f(t,y)=t+y$ pour obtenir ces
deux conditions). Un calcul analogue (p\'enible~!) de 
$\dfrac{\partial^3\Phi}{\partial h^3}$ conduirait au r\'esultat suivant.


\bigskip
\majdeb
\state Th\'eor\`eme -- |  {\it La m\'ethode de Runge-Kutta d\'efinie par
le tableau des coefficients $c_i$, $a_{ij}$, $b_j$ est
\vskip3pt
\noindent$\bullet$ d'ordre $\geq 2$ ssi\ \  
$\displaystyle\sum_j b_jc_j=\dfrac{1}{2}$.
\vskip1pt
\noindent$\bullet$ d'ordre $\geq 3$ ssi\ \  $\displaystyle\sum_j b_jc_j
=\dfrac{1}{2}$~;\quad
$\displaystyle\sum_jb_jc^2_j=\frac{1}{3}$~;\quad 
$\displaystyle\sum_{i,j}b_ia_{ij}c_j=\dfrac{1}{6}$. 
\vskip1pt
\noindent$\bullet$ d'ordre $\geq$ 4 ssi 
$$
\eqalign{
&\sum_j b_jc_j=\frac{1}{2}~;\quad\sum_jb_jc^2_j=\frac{1}{3}~;\quad
\sum_j b_jc^3_j=\frac{1}{4}\cr
&\sum_{i,j}b_ia_{ij}c_j=\frac{1}{6}~;\quad
\sum_{i,j}b_ia_{ij}c^2_j=\frac{1}{12}~;\quad
\sum_{i,j}b_ic_ia_{ij}c_j =\frac{1}{8}~~;\cr
&\sum_{i,j,k}b_ia_{ij}a_{jk}c_k=\frac{1}{12}.\cr
}
$$} 
\vskip-\baselineskip
\majfin
\vskip\baselineskip

Pour la v\'erification pratique, on notera que certaines des
expressions pr\'ec\'edentes sont des produits des matrices
$C=\pmatrix{ c_1\cr \vdots\cr c_q\cr}$,\quad
$A=(a_{ij})$,\quad$B=(b_1b_2\dots b_q)$. Ainsi
$$
\sum_{i,j}b_ia_{ij}c_j=BAC,\qquad
\sum_{i,j,k}b_ia_{ij}a_{jk}c_k=BA^2C.
$$

Pour les exemples du \S\ 3.2, on voit ainsi que la m\'ethode d'Euler
est d'ordre 1, et que les m\'ethodes de l'exemple 2 sont d'ordre 2. De
plus, dans une m\'ethode avec $q=2$, il y a {\it a priori} un seul
coefficient $a_{ij}$ non nul, \`a savoir $\alpha=a_{21}$. On a alors
$c_2=\sum_{j<2}a_{2j}=\alpha$ et la m\'ethode est
d'ordre 2 au moins ssi $\displaystyle\sum b_jc_j=b_2\alpha=1/2$, soit
$b_2=1/2\alpha$ et $b_1=1-b_2=1-1/2\alpha$. On voit donc qu'il n'y
avait pas d'autres choix possibles pour une m\'ethode d'ordre~2 avec
$q=2$.

Enfin, la m\'ethode Runge-Kutta \og classique\fg\ pr\'esent\'ee dans
l'exemple 3 est d'ordre~4 (l'ordre n'est pas $\geq 5$ car
$\displaystyle\sum b_jc_j^4\not=1/5$). C'est si l'on peut dire la \og
m\'ethode reine\fg\ des m\'ethodes \`a un pas~: ordre \'elev\'e, grande
stabilit\'e (gr\^ace \`a la positivit\'e des coefficients, voir
remarque finale du \S$\,$3.3). Il existe des m\'ethodes d'ordre encore
plus \'elev\'e (voir exercice 5.4), mais leur plus grande 
complexit\'e les rend peut-\^etre un peu moins praticables.

\section 4. Contr\^ole du pas|

\noindent La mani\`ere la plus simple pour appliquer une m\'ethode de 
r\'esolution num\'erique consiste \`a utiliser un pas constant $h_n=h$.

La principale difficult\'e est alors de d\'eterminer $h_{\max}$ de
fa\c{c}on que l'erreur globale ne d\'epasse pas une certaine
tol\'erance $\varepsilon$ fix\'ee \`a l'avance~; on ne sait pas en
effet quelle sera l'\'evolution de la solution \'etudi\'ee, de sorte
qu'il est difficile de pr\'evoir {\it a priori} les erreurs de
consistance.

L'utilisation d'algorithmes \`a pas variables pr\'esente de ce point
de vue deux avantages majeurs~:

\medskip
{\parindent1pc
\item{$\bullet$} l'adaptation du pas \`a chaque \'etape permet
  d'optimiser l'erreur commise en fonction de la tol\'erance prescrite
  $\varepsilon$, sous r\'eserve qu'on dispose d'une estimation \og
  instantan\'ee\fg\ de l'erreur de consistance $e_n$.
  
\item{$\bullet$} l'approche d'une discontinuit\'e ou d'une
  singularit\'e de l'\'equation diff\'erentielle ne peut se faire
  g\'en\'eralement qu'avec une r\'eduction importante du pas. Dans
  cette circonstance, il convient d'arr\^eter l'algorithme avant de
  traverser la discontinuit\'e, faute de quoi les erreurs deviennent
  impr\'evisibles. Le calcul du pas sert alors de test d'arr\^et.
\par}


\subsection 4.1. Principe g\'en\'eral du contr\^ole|

\noindent Soit $[t_0,t_0+T]$ l'intervalle de temps consid\'er\'e. On suppose fix\'ee une tol\'erance
$\varepsilon$ pour l'erreur globale
$$
\max_{0\leq n\leq N}|y_n-z(t_n)|.
$$
Supposons \'egalement qu'on dispose d'un estimation $S$ de la constante de stabilit\'e. En
n\'egligeant l'erreur initiale $|y_0-z(t_0)|$ et les erreurs d'arrondi, on a alors
$$
\eqalign{
\max_{0\leq n\leq N}|y_n-z(t_n)|&\leq S\sum_{0\leq n<N}|e_n|=S\sum h_n\frac{|e_n|}{h_n}\cr
&\leq S\Big(\sum h_n\Big)\max\,\Big(\frac{|e_n|}{h_n}\Big)\leq ST\
\max\,\Big(\frac{|e_n|}{h_n}\Big).\cr 
}
$$
Il suffit donc de choisir les pas $h_n$ en sorte que
$$
\max\,\Big(\frac{|e_n|}{h_n}\Big)\leq\delta=\frac{\varepsilon}{ST}~;
$$
$\dfrac{|e_n|}{h_n}$ repr\'esente intuitivement l'erreur de consistance par unit\'e de temps, et
c'est ce rapport qu'il s'agit de contr\^oler.

Il est bien entendu impossible de d\'eterminer exactement $e_n$, sinon on conna\^\i trait du m\^eme
coup la solution exacte par la formule $z(t_{n+1})=y_{n+1}+e_n$~! On va supposer n\'eanmoins qu'on
dispose d'une estimation $e^*_n$ de $e_n$.

Dans la pratique, on se fixe un encadrement $[h_{\min},h_{\max}]$ du pas
($h_{\min}$ est impos\'e par les limitations du temps de calcul et par
l'accroissement des erreurs d'arrondi quand le pas diminue, cf. \S\ 
2.5). On essaie alors de choisir $h_n\in[h_{\min},h_{\max}]$ de fa\c{c}on
que $|e^*_n|/h_n\leq\delta$, et si l'erreur est notablement
inf\'erieure on se permet d'augmenter prudemment $h_n$. Par exemple~:

\medskip
\noindent$\bullet$ si\quad $\dfrac{1}{3}\,\delta\leq\dfrac{|e^*_n|}{h_n}
\leq\delta$, alors\quad $h_{n+1}:= h_n$~;

\noindent$\bullet$ si\quad
$\dfrac{|e^*_n|}{h_n}<\dfrac{1}{3}\,\delta$, alors\quad $h_{n+1}:=
\min (1.25\,h_n,h_{\max})$~;

\noindent$\bullet$ si\quad $\dfrac{|e^*_n|}{h_n}>\delta$, alors\quad
$h_{n+1}:= 0.8\,h_n$, avec arr\^et de l'algorithme si $h_{n+1}<h_{\min}$.

\medskip 

Ce dernier cas peut correspondre \`a l'approche d'une discontinuit\'e,
l'erreur augmentant continuellement malgr\'e la diminution du pas.

\noindent Pour l'initialisation du pas, on prend $h_0=h_{\min}$, \`a moins que 
l'on connaisse une valeur initiale plus appropri\'ee.

\subsection 4.2. Estimation du rapport\quad$\vert e_n\vert/h_n$|

\noindent Pour estimer $e_n$, il n'est pas question d'utiliser les expressions
analytiques des d\'eriv\'ees $\dfrac{\partial^l\Phi}{\partial h^l}$ et
$f^{[l]}$, beaucoup trop co\^uteuses en temps de calcul. On est donc
amen\'e \`a rechercher des estimations {\it ad hoc}, qui n'utilisent
si possible que des quantit\'es d\'ej\`a calcul\'ees par l'algorithme,
et qui ne r\'eclament pas de nouvelle \'evaluation de $f$.

\bigskip
\noindent{\bf $\bullet$ M\'ethode d'Euler}
\medskip

Si $p_n=f(t_n,y_n)$, alors
$$
\eqalign{
p_{n+1}-p_n&=f(t_n+h_n,y_n+h_nf(t_n,y_n))-f(t_n,y_n)\cr
&=h_nf'_t(t_n,y_n)+h_nf(t_n,y_n)f'_y(t_n,y_n)+o(h_n)\cr
&=h_nf^{[1]}(t_n,y_n)+o(h_n).\cr
}
$$
Comme $e_n=\dfrac{1}{2}\ h^2_nf^{[1]}(t_n,y_n)+o(h_n^2)$ on a une approximation de $e_n$ donn\'ee par
$$
\frac{e^*_n}{h_n}=\dfrac{1}{2}\ (p_{n+1}-p_n).
$$
Ceci ne n\'ecessite aucun calcul suppl\'ementaire puisque $p_{n+1}$ est de toute fa\c con
n\'ecessaire pour l'\'etape suivante.

\bigskip
\noindent{\bf $\bullet$ M\'ethode de Runge-Kutta d'ordre 2}
\medskip
$$
\matrix{
0&\Big|&0&0\cr
\noalign{\vskip-2pt}
\alpha&\Big|&\alpha&0\cr
\noalign{\vskip-6pt}
\noalign{\hbox to 7pc{\hrulefill}}
\noalign{\vskip-2pt}
&\Bigg|&\displaystyle 1-\dfrac{1}{2\alpha}&\displaystyle\dfrac{1}{2\alpha}\cr
}
$$
D'apr\`es le \S\ 2.4 on a ici
$$
e_n=h^3_n\Big(\frac{1}{3!}\,f^{[2]}(t_n,y_n)-\frac{1}{2!}\ \frac{\partial^2\Phi}{\partial
h^2}\,(t_n,y_n,0)\Big)+o(h^3_n), 
$$
et les calculs du \S\ 3.4 donnent
$$
\eqalign{
e_n&=h^3_n\Big(\frac{1}{6}\ f^{[2]}-\frac{\alpha}{4}\
(f''_{tt}+2f''_{ty}f+f''_{yy}f^2)\Big)(t_n,y_n)+o(h^3_n)\cr 
&=h^3_n\Big(\Big(\frac{1}{6}-\frac{\alpha}{4}\Big)(f''_{tt}+2f''_{ty}f+f''_{yy}f^2)+\frac{1}{6}\
f'_yf^{[1]}\Big)+o(h_n^3).\cr 
}
$$
On est par ailleurs amen\'e \`a calculer les quantit\'es
$$
\eqalign{
&p_{n,1}=f(t_n,y_n),\cr
&p_{n,2}=f(t_n+\alpha h_n,y_n+\alpha h_np_{n,1}),\cr
&p_{n+1,1}=f\Big(t_n+h_n,y_n+h_n\Big(\Big(1-\frac{1}{2\alpha}\Big)p_{n,1}+\frac{1}{2\alpha}\,p_{n,2}
\Big)\Big).\cr
}
$$
Des d\'eveloppements limit\'es d'ordre 2 donnent apr\`es calcul~:
$$
\eqalign{
p_{n,2}-p_{n,1}&=\alpha
h_nf^{[1]}+\alpha^2\,\frac{h_n^2}{2}\,(f''_{tt}+2f''_{ty}f+f''_{yy}f^2)+o(h^2_n),\cr
p_{n+1,1}-p_{n,1}&=h_nf^{[1]} +h_n\,\frac{1}{2\alpha}\ (p_{n,2}-p_{n,1})f'_y\cr
&\qquad\qquad+\frac{h_n^2}{2}\ (f''_{tt}+2f''_{ty}f+f''_{yy}f^2)+o(h_n^2),\cr
p_{n+1,1}-p_{n,1}&-\frac{1}{\alpha}\ (p_{n,2}-p_{n,1})=(1-\alpha)\ \frac{h_n^2}{2}\
(f''_{tt}+2f''_{ty}f+f''_{yy}f^2)\cr
&\hskip13.26pc+\frac{h^2_n}{2}\,f'_yf^{[1]}+o(h^2_n).\cr 
}
$$
On peut approximer tr\`es grossi\`erement $e_n/h_n$ par
$$
\frac{e^*_n}{h_n}=\frac{1}{3}\ \Big(p_{n+1,1}-p_{n,1}-\frac{1}{\alpha}\,(p_{n,2}-p_{n,1})\Big).
$$
Il n'y a pas de justification th\'eorique s\'erieuse pour cela, l'id\'ee est simplement que les
d\'eveloppements limit\'es se ressemblent formellement.

\bigskip
\noindent{\bf $\bullet$ M\'ethode de Runge-Kutta classique}
\medskip

Il n'y a pas ici de m\'ethode simple permettant d'\'evaluer $e_n$, m\^eme grossi\`erement. On peut
cependant observer que
$$
e_n=h^5_n\times\hbox{(d\'eriv\'ees d'ordre $\leq 4$ de $f$)}~;
$$
D'autre part, des calculs analogues \`a ceux ci-dessus montrent que les quantit\'es
$$
\lambda_n=p_{n,4}-2p_{n,2}+p_{n,1}\qquad{\rm et}\qquad \mu_n=p_{n,3}-p_{n,2}
$$
sont de la forme
$$
h^2_n\times\hbox{(d\'eriv\'ees d'ordre $\leq 2$ de $f$).}
$$
Au lieu de comparer $\dfrac{|e_n|}{h_n}$ \`a $\delta$, on peut essayer de comparer
$$
\lambda^2_n+\mu^2_n\qquad\hbox{\`a}\qquad\delta
$$
ou (plus rapide)
$$
|\lambda_n|+|\mu_n|\qquad\hbox{\`a}\qquad\delta'=\sqrt{\delta}=\sqrt{\frac{\varepsilon}{ST}},
$$
quitte \`a ajuster \'eventuellement les valeurs de $\delta$ et $\delta'$ par tatonnements.

Si l'on d\'esire une \'evaluation plus pr\'ecise de $e_n$, il est n\'ecessaire d'utiliser des
techniques plus \'elabor\'ees, telles que les m\'ethodes de Runge-Kutta embo\^\i t\'ees (voir par
exemple le livre de Crouzeix-Mignot, chapitre 5, \S\ 6).

\section 5. Probl\`emes|

{\parindent1.4pc
\noindent{\bf 5.1.} On \'etudie la m\'ethode num\'erique (M) de 
r\'esolution de l'\'equation
diff\'erentielle $y'=f(x,y)$ d\'efinie par
$$
\eqalign{
&y_{n+1}=y_n+h_n\Phi(t_n,y_n,h_n),\cr
&\Phi(t,y,h)=\alpha f(t,y)+\beta f(t+\dfrac{h}{2},y+\dfrac{h}{2}\ f(t,y))+\gamma f(t+h,y+hf(t,h))\cr
}
$$
o\`u $\alpha,\beta,\gamma$ sont des r\'eels compris entre 0 et 1.

\medskip
\item{(a)} Pour quelles valeurs du triplet $(\alpha,\beta,\gamma)$ 
retrouve-t-on

\smallskip
\item{} $\bullet$ la m\'ethode d'Euler~?

\item{} $\bullet$ la m\'ethode du point milieu~?

\item{} $\bullet$ la m\'ethode de Heun~?

\medskip

\item{(b)} Dans cette question et la suivante, on supposera que la fonction $f(t,y)$ est de classe
$C^{\infty}$ sur $[t_0,t_0+\tau]\times\bR$, et $k$-lipschitzienne en $y$. Pour quelles valeurs
de $(\alpha,\beta,\gamma)$ la m\'ethode propos\'ee est-elle stable~?

\medskip
\item{(c)} Quelles relations doivent satisfaire $(\alpha,\beta,\gamma)$ pour que la m\'ethode soit

\smallskip
\item{}  consistante~?\quad  convergente~?\quad  d'ordre $\geq 1$~?\quad d'ordre $\geq 2$~?

\medskip
La m\'ethode (M) peut-elle \^etre d'ordre sup\'erieur~?

\bigskip
\noindent{\bf 5.2.} On consid\`ere la m\'ethode de Runge-Kutta d\'efinie
par le tableau
$$
\vbox{\offinterlineskip\halign{
\hfill#\hfill\enspace&
\tv\enspace\hfill#\hfill\quad&
\hfill#\hfill\quad&
\hfill#\hfill\cr
0&0&0&0\cr
1/4&1&0&0\cr
3/4&$-$9/20&6/5&0\cr
\noalign{\hrule}
1&1/9&1/3&5/9\cr
}}
$$
On consid\`ere l'\'equation $y'=f(t,y)=t+y+1$. R\'esoudre cette \'equation et calculer
$f^{[n]}(0,0)$ pour tout $n$. Que peut-on dire de $\partial^n\Phi/\partial 
h^n(0,0;0)$~? D\'eterminer l'ordre de cette m\'ethode.

\bigskip
\noindent{\bf 5.3.} On se propose d'obtenir une borne explicite pour 
l'ordre $p$ d'une m\'ethode de
Runge-Kutta dont le nombre de points interm\'ediaires $q$ est fix\'e.
\medskip

\item{(a)} Dans l'espace vectoriel $\cP_n$ des polyn\^omes \`a coefficients r\'eels de degr\'e
inf\'erieur ou \'egal \`a $n$, on consid\`ere la forme bilin\'eaire 
$$
\langle P,Q\rangle=\int^1_0 P(x)Q(x)dx.
$$
\item{} D\'eterminer la matrice de cette forme bilin\'eaire dans la base canonique
$(1,x,\dots,x^n)$. En d\'eduire que la matrice sym\'etrique
$$
M=\pmatrix{
1&1/2&\dots&1/(n+1)\cr
1/2&1/3&&1/(n+2)\cr
&&\vdots&\cr
1/(n+1)&1/(n+2)&&1/(2n+1)\cr
}
$$
\item{} est d\'efinie positive et que det$\,M>0$.

\medskip
\item{(b)} On consid\`ere pour $q\in\bN$ le syst\`eme d'\'equations
$$
{\rm (S)}\quad\cases{
b_1+b_2+\dots b_q=1\hfill\cr
\noalign{\vskip3pt}
b_1c_1+b_2c_2+\dots+b_qc_q=1/2\hfill\cr
\quad \vdots\hfill\cr
b_1c_1^{2q}+\dots+b_qc_q^{2q}=1/(2q+1)\hfill\cr
}
$$
\item{} o\`u les $b_i$ et $c_i$ appartiennent \`a $\bR^*_+$. Dans
$\bR^{q+1}$ on note $F$ l'espace vectoriel engendr\'e par les
vecteurs $(1,c_j,c_j^2,\dots,c^q_j)$, $1\leq j\leq q$.

\smallskip
\itemitem{$(\alpha)$} Quelle est la dimension maximum de $F$~?

\smallskip
\itemitem{$(\beta)$} Montrer  que si (S) avait une solution, les vecteurs
$$
\eqalign{
V_1&=(1\qquad 1/2\dots1/(q+1))\cr
V_2&=(1/2\quad 1/3\dots1/(q+2))\cr
V_{q+1}&=(1/(q+1)\dots 1/(2q+1))\cr
}
$$
appartiendraient \`a $F$. En d\'eduire que det$(V_1,V_2,\dots,V_{q+1})=0$ 
puis que le syst\`eme (S) est impossible.
\medskip

\item{(c)} On consid\`ere une m\'ethode de Runge-Kutta
$$
\vbox{\offinterlineskip\halign{
#\quad&
\tv\quad#\hfill\cr
$c_1$&\cr
\vdots&$A=(a_{ij})$\cr
\noalign{\vskip-2pt}
\vdots&\qquad $1\leq i$, $j\leq q$\cr
$c_q$&\cr
\noalign{\hrule}
&$b_1\dots\ \ldots\ b_q$\cr
}}
$$
\item{} associ\'ee \`a la m\'ethode d'int\'egration (INT) $\displaystyle\int^1_0
f(x)dx\simeq\sum^q_{j=1}b_jf(c_j)$.

\item{} On suppose que cette m\'ethode est d'ordre $p$.

\smallskip
\itemitem{$(\alpha)$} Montrer que (INT) est d'ordre $p-1$.

\smallskip
\itemitem{$(\beta)$} En d\'eduire que $p\leq 2q$.

\bigskip
\noindent{\bf 5.4.} On consid\`ere une m\'ethode \`a un pas de la forme
$$
y_{n+1}=y_n+h_n\Phi(t_n,y_n,h_n)\leqno{\rm (M)}
$$
qu'on suppose \^etre d'ordre $p$. On se donne par ailleurs une m\'ethode d'int\'egration approch\'ee
$$
\int^1_0 g(u)du\simeq\sum_{1\leq j\leq q}b_jg(c_j)\leqno{\rm (I)}
$$
d'ordre $p$ au moins (il en existe pour $p$ quelconque).
\smallskip

\item{(a)} On consid\`ere la m\'ethode \`a un pas avec points interm\'ediaires
$t_{n,i}=t_n+c_ih_n$ d\'efinie comme suit
$$
\cases{
\left[\matrix{
t_{n,i}=t_n+c_ih_n\hfill\cr
y_{n,i}=y_n+c_ih_n\Phi(t_n,y_n,c_ih_n)\hfill\cr
p_{n,i}=f(t_{n,i},y_{n,i})\hfill\cr
}\right]\quad 1\leq i\leq q\cr
\noalign{\vskip2pt}
\hskip5pt t_{n+1}=t_n+h_n\hfill\cr
\noalign{\vskip2pt}
\displaystyle\hskip4pt y_{n+1}=y_n+h_n\sum_{1\leq j\leq q}b_j p_{n,j}.\hfill\cr
}\leqno{\rm (M')}
$$
Montrer que (M$'$) est d'ordre $\geq p+1$.
\smallskip
\item{(b)} Gr\^ace \`a un raisonnement par r\'ecurrence, montrer qu'il existe des m\'ethodes de
Runge-Kutta d'ordre $p$ arbitrairement \'elev\'e.
\par}

\bye

