Intégration numérique d'une EDO:
schéma d'Euler implicite
Ce TD est composé de deux parties: la programmation d'un intégrateur implicite et la programmation d'un intégrateur Runge-Kutta d'ordre 4.
Le problème
Dans la première partie, nous allons traiter le mouvement d'un pendule
rigide, comme au TD précédent: une masse m est reliée à
un pivot par une tige de longueur l, supposée rigide et de
masse nulle. Le pendule se trouve dans un champ de pesanteur homogène
g.
Dans la deuxième partie, nous allons considérer le mouvement d'une planète autour du soleil.
Pas à pas
Partie I
-
Les équations d'Euler implicites sont
thetak+1 = thetak + dt*omegak+1 omegak+1 = omegak - dt*sin(thetak+1)
Il n'est pas possible d'écrire la solution explicitement. Mais on peut observer qu'en considérantf=(f1,f2)avecf1(x1,x2) = thetak + dt*x2 f2(x1,x2) = omegak - dt*sin(x1)
on a une fonction vectorielle dontx∞=(thetak+1,omegak+1)est un point fixe. Calculez la Jacobienne(dfi/dxj)i,j=1..2de la fonction vectorielle(f1,f2)en(thetak+1,omegak+1)et montrez que ses valeurs propres ont une norme inférieure à un ( |lambda1/2|<1 ) pour un choix convenable dedt. -
Cela signifie que
f(x)est contractante dans un voisinage du point fixex∞, et qu'on peut trouver le point fixe en itérant la fonction:f(f(f(f(...f(f(...))))))) → x∞
À quelle vitesse s'approche-t-on du point fixe, autrement dit de combien diminue l'écart au point fixe à chaque application supplémentaire def? En pratique, combien de fois doit-on itérer au moins ? -
Écrivez un programme C
qui utilise ce fait pour trouver à chaque pas de temps le nouvel état
du système
(thetak+1,omegak+1)par solution itérative de l'équation implicite. -
Calculez une trajectoire pour les conditions initiales
theta0=0,omega0=1.9sur 50 unités de temps en pas de 0.01 et 0.05. Comparez le comportement de l'intégrateur par rapport au schéma d'Euler simple.
Partie II
- Écrivez les équations du mouvement d'une planète autour du soleil, qu'on supposera fixe à l'origine (le rapport m/M des masse de la terre et du soleil par exemple est d'ordre 10-7, ce qui justifie l'approximation consistant à dire que le soleil n'est pas influencé par le mouvement de la terre). Posez GM = 1, où G est la constante de gravitation et M la masse du soleil.
- Écrivez un programme utilisant le schéma d'intégration de Runge et Kutta d'ordre 4 pour simuler le mouvement de la planète.
- Simulez au moins deux cas, planète en orbite circulaire et planète en orbite elliptique, en choisissant bien les conditions initiales.