Méthodes numériques pour quelques EDO

Un problème de tir

À t=0, on lance un projectile (assimilé à un point de masse m) depuis un point repéré par son altitude y=0 et sa position x=0, z=0. La vitesse initiale du projectile est v = (vx0, vy0, 0). En négligeant les effets de rotation de la terre (force de Coriolis), la trajectoire du projectile s'effectue dans le plan défini par z=0. En négligeant en outre les effets de frottements dûs à l'air, on trouve l'équation du mouvement suivante (donnée par les lois de la dynamique de Newton):

où g est la constante de gravité. La solution exacte de ce problème vaut:

ce qui montre que le projectile suit une parabole et retombe au sol au temps .

Question 1

Écrivez un programme contenant tous les paramètres physiques du problème et qui calcule la solution exacte pour un nombre N (au début on prendra N=100) de valeurs de t entre 0 et tmax. On prendra les valeurs suivantes:

On créera un module comme suit:

module tir

  implicit none

  ! les parametres physiques en parametre (par exemple g)
  real*8 , parameter :: g = 10.0D0

  ! etc
  
end module

program prog_tir

  use tir

  implicit none

  integer :: N

  N = 100
  ! etc
end program

On utilisera la double précision (real*8). On écrira la solution exacte dans un fichier avec le format suivant:

x0 y0
x1 y1
...
xN yN

La première colonne du fichier contient les abscisses x, la seconde colonne les ordonnées y aux différents temps. On lancera ensuite ipython pour visualiser les résultats:

ipython --pylab

Si le fichier s'appelle data.txt, on écrira dans la console ipython:

u = loadtxt('data.txt')
plot(u[:, 0], u[:, 1])

On devrait voir la trajectoire du projectile.

Question 2

On considère le schéma d'Euler explicite pour l'EDO autonome qui s'écrit:

où Δt est le pas de temps. Programmer le schéma d'Euler explicite en écrivant une sous-routine AdvanceEuler (toutes les sous-routines seront implémentées dans le module tir) sous la forme suivante (pour traiter une EDO générique) :

// entree :
// m : taille du vecteur X 
// Xn : X au temps tn
// dt : pas de temps Delta t
// F : fonction F sous la forme d'une sous-routine

// sortie :
// Xnext : X au temps tn+1
subroutine AdvanceEuler(m, dt, Xn, F, Xnext)

  implicit none

  integer :: m
  real*8 :: dt, Xn(*), Xnext(*)

  interface
    subroutine F(X, Fx)
      real*8 :: X(*)
      real*8 :: Fx(*)
    end subroutine
  end interface

  ! variables locales
  real*8 :: Fx(m)

  ! si vous voulez evaluer F
  call F(Xn, Fx)

  ! ensuite il faut calculer Xnext

end subroutine

// example de fonction F
// en entree X, en sortie Fx
subroutine chute_libre(X, Fx)

    implicit none

    real*8 :: X(*), Fx(*)

    // il faut remplir Fx
end subroutine
 

Dans le programme principal, on appellera AdvanceEuler comme suit:

dt = tmax / (N-1)
Xn = (/ 0.0D0, 0.0D0, Vx0, Vy0 /)
do i = 1, N
  call AdvanceEuler(4, dt, Xn, chute_libre, Xnext)
  Xn = Xnext
end do

On remarquera qu'on ne stocke pas tous les instantanés, ils sont écrasés au fur et à mesure. On a ici appelé la sous-routine calculant F(X) (sous-routine à créer) chute_libre. Écrire la solution numérique dans un fichier. On pourra afficher la trajectoire de la solution numérique sous ipython (plot superpose les différentes courbes par défaut) en même temps que la solution exacte. On pourra ici calculer l'erreur relative en tapant:

norm(u - u_numerique) / norm(u)

Vérifier que si vous multipliez N par 2, l'erreur relative est divisépar deux.

Question 3

Programmer à présent (avec une routine AdvanceRK2 par exemple) un schéma de Runge-Kutta d'ordre 2, par exemple celui du point-milieu, donné par

Calculer l'erreur relative obtenue avec ce schéma, pourquoi obtient-on une telle précision ?

Question 4

Programmer à présent (avec une routine AdvanceRK4 par exemple) un schéma de Runge-Kutta d'ordre 4

Calculer l'erreur relative obtenue avec ce schéma.

Question 5

On tient compte à présent des frottements subis par le projectile. On suppose que l'air exerce sur le solide une force f = - k v qui s'oppose donc au mouvement. Si le coefficient de frottement k est constant, on peut encore résoudre le problème analytiquement. En réalité, le coefficient de frottement est plutôt non-linéaire, par exemple . Il est alors beaucoup plus difficile de résoudre le problème analytiquement. Écrire le nouveau système différentiel et écrire une sous-routine chute_amortie par exemple implémentant ce système. On prendra le paramètre et on arrêtera le calcul (ou l'écriture) quand le projectile touche le sol. Calculer la solution avec le schéma d'Euler et de RK2, afficher la solution obtenue sous ipython. Pour comparer la performance des deux schémas, on pourra calculer une solution de référence avec RK4 et calculer l'erreur relative en tapant :

norm(u_rk4 - u_euler) / norm(u_rk4)

Si on veut atteindre une erreur relative de 1%, quel N faut-il prendre pour le schéma d'Euler ? Quel N faut-il pour le schéma du point-milieu ? Quel schéma est donc le plus rapide pour une précision fixée ?

Question 6

On augmente maintenant les frottements en prenant . Constater que le schéma d'Euler explicite donne une solution non-physique pour N=101 (le projectile repart en arrière). On pourra regarder ce qui se passe pour N=181, puis N=101.

Mécanique céleste : orbites planétaires

Le but de cet exercice est de calculer numériquement les trajectoires autour du soleil de 5 planètes du système solaire : Jupiter, Saturne, Uranus, Neptune et Pluton. Il est tiré du polycopié d'Analyse Numérique, de Ernst Hairer.

On considère que chacune des planètes n'est soumise qu'à l'attraction gravitationnelle exercée par toutes les autres, plus celle du soleil. On rappelle que la force gravitationnelle exercée par un corps B de masse mB situé en sur un corps A de masse mA situé en est

où G est la constante de gravitation universelle. La loi fondamentale de la dynamique nous permet alors d'écrire les équations suivantes :

pour i=1 à 5, où xi, mi sont les positions et masses de Jupiter (i=1), Saturne(i=2), Uranus (i=3), Neptune (i=4) et Pluton (i=5). Le mouvement de ces planètes est décrit dans le référentiel lié au soleil qui a donc la position fixe x0=0. Les inconnues du problème sont donc les positions xi(t) pour i=1 à 5. En notant vi = xi'< la vitesse de la planète numéro i, on peut réécrire les relations précédentes sous forme d'un système d'équations différentielles du premier ordre :

pour i = 1 à 5.

La masse du soleil est m0 = 1.00000597682 (modifiée pour tenir compte de la masse des planètes proches) et la constante de gravitation universelle est G = 2.95912208286e-4. Les unités de masse sont relatives à celle du soleil, les distances en unités astronomiques (1 A.U. = 149 597 870 km) et le temps est compté en jours terrestres. Pour les 5 planètes, les positions et vitesses initiales calculées le 5 septembre 1994 sont données dans la table ci-dessous:

x1 -3.5023653 -3.8169847 -1.5507963
x2 9.0755314 -3.0458353 -1.6483708
x3 8.3101420 -16.2901086 -7.2521278
x4 11.4707666 -25.7294829 -10.816945
x5 -15.5387357 -25.2225594 -3.1902382
v1 0.00565429 -0.00412490 -0.00190589
v2 0.00168318 0.00483525 0.00192462
v3 0.00354178 0.00055029 0.00055029
v4 0.00288930 0.00039677 0.00039677
v5 0.00276725 -0.00136504 -0.00136504

Les masses sont :

Le temps maximum est tmax = 100 000, temps au bout duquel toutes les planètes ont fait au moins une révolution complète autour du soleil.

On demande de résoudre numériquement ce système avec les méthodes d'Euler explicite, de RK2 et RK4 (codées dans le premier exercice). On écrira dans un fichier les positions des 5 planètes (on aura donc 15 colonnes). Pour visualiser la trajectoire des planètes, on tapera dans ipython :

u = loadtxt("trajectoire.txt")

from mpl_toolkits.mplot3d import Axes3D
fig = figure(); ax = Axes3D(fig)
ax.plot3D(u[:,0], u[:,1], u[:,2])
ax.plot3D(u[:,3], u[:,4], u[:,5])
ax.plot3D(u[:,6], u[:,7], u[:,8])
ax.plot3D(u[:,9], u[:,10], u[:,11])
ax.plot3D(u[:,12], u[:,13], u[:,14])

Question 7

Tester le schéma d'Euler explicite, observez qu'il est instable quel que soit le pas de temps choisi (essayez Δ t = 300, puis Δ t = 10) : toutes les planètes quittent leur orbite (et leur trajectoire part plus ou moins vite vers l'infini). Cela est d'autant plus visible que les planètes font beaucoup de tours autour du soleil (donc celles qui en sont le plus proche comme Jupiter et Saturne).

Question 8

Tester RK2 et RK4, commenter.