MENUCours D'introduction à l'analyse numérique

Les méthodes classiques utilisées pour approcher les solutions des équations différentielles (Euler, Runge-Kutta ou Verlet) présentent l'inconvénient d'avoir un pas fixe. Or, dans de nombreuses situations, une grandeur physique peut varier rapidement un moment, puis plus lentement ensuite. Il est alors judicieux d'adapter le pas en fonction des besoins, c'est-à-dire de le diminuer lorsque les variations sont rapides et de l'augmenter dans l'autre cas, de façon à économiser du temps de calcul. C'est la philosophie des algorithmes à pas adaptatif dont la méthode décrite ici fait partie.

Schéma d'Euler-Richardson

Généralités

Tout d'abord, précisons que l'on cherche à approcher la solution d'une équation différentielle avec condition initiale qui se met sous la forme \[\left\{\begin{array}{lccc} \dot{\mathbf{y}} & = & f(t,\mathbf{y}(t)), & 0\leq t \leq T \\ \mathbf{y}(0) & = & \mathbf{y}_{0} & \\ \end{array}\right.\] avec \(\mathbf{y}(t)\) un scalaire ou un vecteur et \(f(t,\mathbf{y}(t))\) une fonction suffisamment régulière pour que l'existence et l'unicité de la solution ne pose pas de problème.

À l'instar de la méthode d'Euler, la méthode d'Euler-Richardson est basée sur la discrétisation de la variable \(t\) ainsi que sur l'utilisation de différences finies pour approcher les dérivées. La différence réside ici dans l'utilisation d'un pas variable, dont la valeur est liée à la précision exigée par l'utilisateur. Ainsi, l'intervalle \([0,T]\) est divisé en subdivisions de longueur variable \(h_n\). On a la relation \[ t_n=t_{n-1}+h_n\] Le nombre de points n'est alors connu qu'a posteriori.

Rappelons enfin, qu'avant d'effectuer la résolution numérique, il est conseillé de procéder à un travail de simplification qui passe par la définition d'un nouveau système d'unités adapté au problème. On commence donc par ré-exprimer l'équation différentielle à l'aide de grandeurs sans dimension (voir La méthode d'Euler).

Principe de l'algorithme

Commençons par rappeler la méthode d'Euler. La valeur \(\mathbf{y}_{n+1}^{e1}\) correspondant à l'instant \(t_{n+1}\) est donnée par \[\mathbf{y}_{n+1}^{e1}=\mathbf{y}_{n}+h\,f(t_{n},\mathbf{y}_{n})\] où \(h\) est le pas \(t_{n+1}-t_n\). Nous savons que la méthode d'Euler est d'ordre un. Elle produit une erreur de troncature \(\varepsilon_{1}\simeq kh^{2}\) où \(k\) est une constante liée aux valeurs de la dérivée de \(f(t,\mathbf{y}(t))\).

Envisageons maintenant une évolution en deux étapes de pas \(h/2\), ce qui nous donne une seconde évaluation que l'on note \(\mathbf{y}_{n+1}^{e2}\) : \[\mathbf{y}_{n+1}^{e2}=\mathbf{y}_{n}+\frac{h}{2}f(t_{n},\mathbf{y}_{n})+\frac{h}{2}f(t_{n+1/2},\mathbf{y}_{n+1/2})\] ce qui produit une erreur \[\varepsilon_{2}\simeq 2\times k\left(\frac{h}{2}\right)^{2}=k\frac{h^{2}}{2}=\frac12 \varepsilon_{1}\] Ainsi, on peut espérer que les erreurs se compenseront (tout du moins à l'ordre 2) si l'on calcule \(\mathbf{y}_{n+1}\) de la façon suivante : \[\mathbf{y}_{n+1}=2\,\mathbf{y}_{n+1}^{e2}-\mathbf{y}_{n+1}^{e1}= \mathbf{y}_{n}+h\,f(t_{n+1/2},\mathbf{y}_{n+1/2})\] d'où le schéma numérique suivant : \[ \left\{ \begin{array}{lcl} \mathbf{k}_{n} & = & f(t_{n},\mathbf{y}_{n}) \\ \mathbf{k}'_{n} & = & f(t_{n}+\frac{h}{2},\mathbf{y}_{n}+\frac{h}{2}\mathbf{k}_{n}) \\ \mathbf{y}_{0} & = & \mathbf{y}(0) \\ \mathbf{y}_{n+1} & = & \mathbf{y}_{n}+h\,\mathbf{k'}_{n}, \\ \end{array}\right. \]

En prime, on peut accéder à la quantité \[\varepsilon_{2}=\left|\mathbf{y}_{n+1}^{e2}-\mathbf{y}_{n+1}^{e1}\right|= \frac{h}{2}\left|f(t_{n+1/2},\mathbf{y}_{n+1/2})-f(t_{n},\mathbf{y}_{n})\right |= \frac{h}{2}\left|\mathbf{k}'_n-\mathbf{k}_n\right|\] qui représente l'erreur que l'on fait avec la méthode d'Euler lorsque l'on multiplie par deux le pas (\(h/2\to h\)). Cette erreur va nous servir d'indicateur pour savoir si l'on peut augmenter le pas ou s'il faut le diminuer. Il suffit alors de fixer un seuil de précision (notons le \(\varepsilon_{\text{seuil}}\)). Si l'erreur est plus grande que \(\varepsilon_{\text{seuil}}\) d'un facteur \(\alpha\), il suffit de réduire le pas d'un facteur \(\sqrt{\alpha}\) puisque \(\varepsilon_2\) varie comme \(h^2\). Symétriquement, si l'erreur est plus petite que le seuil d'un facteur \(\alpha\), alors on s'autorise à augmenter le pas d'un facteur \(\sqrt{\alpha}\), l'idée étant de ne pas gaspiller du temps de calcul. Résumons maintenant la méthode :

Exemple : Le mouvement képlerien

Illustrons cette méthode sur l'exemple du mouvement d'un astre A en orbite elliptique autour du Soleil S. Dans le cas où l'ellipse est de forte excentricité, la distance \(r=AS\) varie de façon importante entre le périhélie (point le plus proche du Soleil) et l'aphélie (point le plus éloigné) ce qui signifie que la force de gravitation -inversement proportionnelle au carré de la distance \(AS\)- varie sur de grandes échelles. On comprend dès lors qu'une méthode à pas variable convienne mieux qu'une méthode classique.

Travail préliminaire

Si l'on note \(m_{\star}\) la masse du Soleil, et \(m\ll m_{\star}\) celle de l'astre, le mouvement de ce dernier dans le référentiel héliocentrique Référentiel que nous supposons galiléen pour notre propos. En réalité, le Soleil oscille légèrement autour du centre de masse du système solaire. est régi par l'équation du mouvement : \begin{equation} m\frac{\mathrm{d}^{2} \overrightarrow{\rm SA}}{\mathrm{d}t^{2}}=-\frac{\mathcal{G}mm_{\star}}{r^{3}} \overrightarrow{\rm SA} \end{equation} Adoptons le système d'unités astronomiques (u.a) dans lequel l'unité de longueur est le demi-grand axe de l'orbite terrestre \(a\simeq 150.10^{6}\) km, l'unité de temps est la période orbitale terrestre \(T_{0}=1\) an et l'unité de masse, celle du Soleil \(m_{\star}\). La troisième loi de Kepler nous donne la valeur de \(\mathcal{G}\) dans ce système d'unités \[\frac{a^{3}}{{T_{0}}^{2}}=\frac{\mathcal{G}m_{\star}}{4\pi^{2}}\quad\Longrightarrow\quad \mathcal{G}=4\pi^{2}\] ce qui permet de reformuler l'équation (1) dans le plan \((x,y)\) : \[ \ddot{x} = -\dfrac{4\pi^{2}x}{(x^{2}+y^{2})^{3/2}} \quad\text{et}\quad \ddot{y} = -\dfrac{4\pi^{2}y}{(x^{2}+y^{2})^{3/2}} \] Il s'agit donc de deux équations scalaires du second ordre que l'on peut transformer en une équation vectorielle du premier ordre, de la forme \(\dot{\mathbf{y}} = f(\mathbf{y})\) si l'on pose : \[\mathbf{y}=\left(\begin{array}{c}x\\y\\v_{x}\\v_{y}\end{array}\right) \quad\text{et}\quad f(\mathbf{y})=\left(\begin{array}{c}v_{x}\\v_{y}\\[1mm]-\dfrac{4\pi^{2}x}{(x^{2}+y^{2})^{3/2}}\\[3mm] -\dfrac{4\pi^{2}y}{(x^{2}+y^{2})^{3/2}} \end{array}\right) \] Notons qu'ici la fonction d'évolution \(f(\mathbf{y})\) ne dépend pas explicitement du temps (on dit que l'équation différentielle est autonome).

Résultat

Supposons qu'à l'instant initial \(t=0\), l'astre soit lancé depuis le point \((x=0,5,y=0)\) avec un vecteur vitesse \((v_x=0,v_y=12)\) (dans les unités astronomiques). Observons le résultat correspondant à une durée \(T=4,768\) ans, un pas initial \(h=0,04\) (envirion 15 jours) et une précision \(\varepsilon_{\text{seuil}}=0,018\). On applique la méthode d'Euler-Richardson :

1. Initialisation : \(h=0,04\), \(T=4,768\) ans et \(\varepsilon_{\text{seuil}}=0,018\).
2. Conditions initiales : \(t=0\quad \text{et}\quad\mathbf{y}=\left(\begin{array}{l}y_1=0,5\\y_2=0\\y_3=0\\y_4=11,5\end{array}\right)\)
3. Tant que \(t\leq T\) faire :
   1. Calcul de \(\mathbf{k}=f(\mathbf{y})\) : \[\mathbf{k}=\left(\begin{array}{c} y_3\\ y_4\\ -\dfrac{4\pi^{2}y_1}{({y_1}^2+{y_2}^2)^{3/2}}\\ -\dfrac{4\pi^{2}y_2}{({y_1}^2+{y_2}^2)^{3/2}}\end{array}\right)= \left(\begin{array}{c} k_1\\ k_2\\ k_3\\ k_4 \end{array}\right)\]
   2. Calcul de \(\mathbf{k'}=f(\mathbf{y}+\frac{h}{2}\mathbf{k})\) : \[\mathbf{k'}=\left(\begin{array}{c} y_3+(h/2)k_3\\ y_4+(h/2)k_4\\ -\dfrac{4\pi^2[y_1+(h/2)k_1]}{\left[(y_1+(h/2)k_1)^2+(y_2+(h/2)k_2)^2\right]^{3/2}}\\ -\dfrac{4\pi^{2}[y_2+(h/2)k_2]}{\left[(y_1+(h/2)k_1)^2+(y_2+(h/2)k_2)^2\right]^{3/2}} \end{array}\right)= \left(\begin{array}{c} k'_1\\ k'_2\\ k'_3\\ k'_4 \end{array}\right)\]
   3. Calcul de l'erreur \(\varepsilon=|\mathbf{k}-\mathbf{k'}|\frac{h}{2}\) et \(\alpha\) : \[\varepsilon=\frac{h}{2}\left(|k'_1-k_1|+|k'_2-k_2|+|k'_3-k_3|+|k'_4-k_4|\right) \quad\text{et}\quad \alpha=\varepsilon/\varepsilon_{\rm seuil}\]
   4. Si \(\alpha>1\) faire \(h=0,9\,h/\sqrt{\alpha}\).
   5. Si \(\alpha<1\) faire \(\mathbf{y}=\mathbf{y}+\mathbf{k'}h\) ; \(h=0,9\,h/\sqrt{\alpha}\) ; \(t=t+h\) ; enregistrement des données.

La Fig. 1 montre le résultat et le compare à celui obtenu avec la méthode de Runge-Kutta. Dans les deux cas, le nombre de pas est le même : \(N=119\). Avec la méthode d'Euler Richardson, l'astre fait un tour en 4,768 ans alors qu'avec Runge-Kutta il ne met que 4,218 ans.

Théoriquement, on prévoit une trajectoire elliptique de grand-axe \(2a=5,675\) ua et de période \(T_{A}=4,779\) ans. Malgré le faible nombre de points, la méthode d'Euler-Richardson donne un résultat tout à fait acceptable (0,2% d'erreur sur la période). Les points sont d'autant plus denses que l'astre est près du périhélie car l'attraction et donc l'accélération, est grande. En revanche, la méthode de Runge-Kutta est décevante (11% d'erreur sur la période). Le pas étant trop grand près du périhélie, la trajectoire est instable et l'astre a tendance à spiraler en convergeant sur le Soleil pour quitter le système solaire au bout d'un certain temps.

En conclusion, le calcul numérique est un outil intéressant à condition d'avoir un regard critique sur le problème et les résultats. Bien que les méthodes classiques soient très utilisées, il faut savoir qu'il existe de nombreuses méthodes numériques plus ou moins sophistiquées, les plus précises étant les méthodes à pas adaptatifs comme la méthode que l'on vient de présenter.