Gravitation

Qu'est-ce que vous appelez un système périodiques ? Quand les points sont récurrents ? Quand l'orbite est vraiment périodique ? Parce que c'est une propriété vraiment très spéciale d'être périodique et en ce sens, quasiment tout serait chaotique …

L'apériodicité du comportement à long terme, c'est le fait que la solution du système ne tende à l'infini ni vers un point fixe ni vers un motif périodique ou quasi-périodique.

D'après ce lien, déterminisme, sensibilité aux conditions initiales et apériodicité du comportement à long terme (éventuellement complétées d'autres conditions plus marginales) seraient toutes trois nécessaires (et suffisantes).

Ok merci de la précision (histoire de débattre sur autre chose :p).

Après c'est une histoire de définitions. Pourquoi pas le définir comme ça mais est-ce que c'est intéressant ?

Généralement ce qui compte quand on parle de système chaotique c'est plutôt la dépendance aux conditions initiales, non ?

@Ryuuken : je ne sais pas pourquoi j'ai mis un "généralement" dans ma parenthèse, en fait.

Généralement ce qui compte quand on parle de système chaotique c'est plutôt la dépendance aux conditions initiales, non ?

Non. Voir l'exemple plus haut. Cette condition est faible, et finalement d'intérêt moindre pour l'étude d'un système chaotique. L'apériodicité est beaucoup plus problématique pour étudier ce genre de système. Dans la vraie vie, on n'a jamais les conditions initiales réelles, et comme le système y est fortement dépendant, s'y intéresser n'est pas forcément pertinent. Ce qui est beaucoup plus intéressant, c'est de chercher des tendances dans le comportement apériodique du système (éventuellement dans l'espace des conditions initiales), voir si on est capable de paramétriser son comportement moyen, etc.

Ok je comprends mieux ton point de vue …

M'enfin pour moi ton exemple plus haut n'a rien de problématique, ça tend vers l'infini et puis c'est tout :-).

M'enfin pour moi ton exemple plus haut n'a rien de problématique, ça tend vers l'infini et puis c'est tout :-).

Je n'ai pas dit qu'il était problématique. Il prouve simplement que la dépendance aux conditions aux limites n'est pas propre aux systèmes chaotiques.

J'avais fais un petit code dessus l'année dernière avec un schéma d'Euler explicite. Mon code est en GNU/Octave et est pas bien foutu, mais en bourrinant avec un pas de temps tout petit c'est presque stable ! Preuve en vidéo : hébergé ici

En vrai, les planètes s'écartent les unes des autres, mais j'ai pas su si c'était du aux conditions initiales trop imprécisent, à la methode de calcul ou à un éloignement réel (les 3 pouvant jouer de concert…).

Code : Schéma euleur explicite mais avec une bourde dans l'adimensionnement, il me semble que j'avais oublié d'adimensionner le temps.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108

% Distance Terre-Soleil = 151e6 km
% Vitesse de la Terre = 107217 km/h
% Vitesse de jupiter = 47006 km/h
% Masse du Soleil = 1.989e30 kg
% Masse de la Terre = 5.972e24 kg
% Masse de la Lune = 0.0123 mT
% Distance T-L = 0.00257 u.a.
% Vitesse de la lune = 3679.2 km/h + vT = 110896.2 km/h = 1.0343 vT
temps=10000001
pas=1e-5

X=zeros(2,2,3);
V=zeros(2,2,3);
XS=zeros(2,2);
VS=zeros(2,2);

% Adimensionnement des paramètres
d_ST=1 ; % Distance adim
mS=1 ; % Masse adim
mT=5.972e24/1.989e30; % masse Terre adim 
vT=1 ; % Vitesse adim
G=6.67e-11*1.989e30/(150e9 * (29783)^2) ;%G*mS/d_ST*vT² constante de gravitation adim

% Conditions initiales
X(:,1,1)=[d_ST ; 0]                   ; 
V(:,1,1)=[0 ; vT ]                        ;
mA= mT                                ;

X(:,1,2)=[(1+0.00257)*d_ST ; 0]       ;
V(:,1,2)=[0 ; 1.0343*vT]              ;
mB= 0.0123*mT                             ;

X(:,1,3)=[5.3*d_ST ; 0]               ; 
V(:,1,3)=[0 ; 0.4384*vT ]                 ;
mC= 313*mT                            ;

XS(:,1)=[0 ; 0]                       ;
VS(:,1)=[0 ; 0.01*vT]                     ;
mS= mS                                ;

dSA=((X(1,1,1)-XS(1,1))^2+(X(2,1,1)-XS(2,1))^2)^(1/2);
dSB=((X(1,1,2)-XS(1,1))^2+(X(2,1,2)-XS(2,1))^2)^(1/2);
dSC=((X(1,1,3)-XS(1,1))^2+(X(2,1,3)-XS(2,1))^2)^(1/2);
dAB=((X(1,1,1)-X(1,1,2))^2+(X(2,1,1)-X(2,1,2))^2)^(1/2);
dAC=((X(1,1,1)-X(1,1,3))^2+(X(2,1,1)-X(2,1,3))^2)^(1/2);
dBC=((X(1,1,2)-X(1,1,3))^2+(X(2,1,2)-X(2,1,3))^2)^(1/2);

tic
fdata=fopen('3planetes.dat', 'w+')

fprintf(fdata, '%f ', X(:,1,:)); % 6 colonnes
fprintf(fdata, '%f ', XS(:,1));  % 2 colonnes
fprintf(fdata, '%f ', dSA);    % 1 colonnes
fprintf(fdata, '%f ', dSB);        % 1 colonnes
fprintf(fdata, '%f ', dSC);        % 1 colonnes
fprintf(fdata, '%f ', dAB);        % 1 colonnes
fprintf(fdata, '%f ', dAC);    % 1 colonnes
fprintf(fdata, '%f ', dBC);    % 1 colonnes
fprintf(fdata, '\n');

for t=2:temps

    %% Déplacement de A
    V(:,2,1)=V(:,1,1) - G * pas * ( mS*(X(:,1,1)-XS(:,1))/(dSA)^3 + mB*(X(:,1,1)-X(:,1,2))/(dAB)^3 + mC*(X(:,1,1)-X(:,1,3))/(dAC)^3 ) ;
    X(:,2,1)=X(:,1,1) + V(:,1,1)*pas;
    
    %% Déplacement de B
    V(:,2,2)=V(:,1,2) - G * pas * ( mS*(X(:,1,2)-XS(:,1))/(dSB)^3 + mA*(X(:,1,2)-X(:,1,1))/(dAB)^3 + mC*(X(:,1,2)-X(:,1,3))/(dBC)^3 ) ;
    X(:,2,2)=X(:,1,2) + V(:,1,2)*pas;
 
    %% Déplacement de C
    V(:,2,3)=V(:,1,3) - G * pas * ( mS*(X(:,1,3)-XS(:,1))/(dSC)^3 + mA*(X(:,1,3)-X(:,1,1))/(dAC)^3 + mB*(X(:,1,3)-X(:,1,2))/(dBC)^3 ) ;
    X(:,2,3)=X(:,1,3) + V(:,1,3)*pas;
  
    %% Déplacement de S
    VS(:,2)=VS(:,1) - G * pas * ( mA*(XS(:,1)-X(:,1,1))/(dSA)^3 + mB*(XS(:,1)-X(:,1,2))/(dSB)^3 + mC*(XS(:,1)-X(:,1,3))/(dSC)^3 ) ;
    XS(:,2)=XS(:,1) + VS(:,1)*pas;

    %% Distance entre les astres
  dSA=((X(1,2,1)-XS(1,2))^2+(X(2,2,1)-XS(2,2))^2)^(1/2);
  dSB=((X(1,2,2)-XS(1,2))^2+(X(2,2,2)-XS(2,2))^2)^(1/2);
  dSC=((X(1,2,3)-XS(1,2))^2+(X(2,2,3)-XS(2,2))^2)^(1/2);
  dAB=((X(1,2,1)-X(1,2,2))^2+(X(2,2,1)-X(2,2,2))^2)^(1/2);
  dAC=((X(1,2,1)-X(1,2,3))^2+(X(2,2,1)-X(2,2,3))^2)^(1/2);
  dBC=((X(1,2,2)-X(1,2,3))^2+(X(2,2,2)-X(2,2,3))^2)^(1/2);

    %% Ecriture tous les 500 pas de temps uniquement
  if mod(t, 500)==1
        fprintf(fdata, '%f ', X(:,2,:));
      fprintf(fdata, '%f ', XS(:,2));
        fprintf(fdata, '%f ', dSA);
        fprintf(fdata, '%f ', dSB);
        fprintf(fdata, '%f ', dSC);
        fprintf(fdata, '%f ', dAB);
        fprintf(fdata, '%f ', dAC);
        fprintf(fdata, '%f ', dBC);
      fprintf(fdata, '\n');
    end
    
    % Décalage des posisitions des matrices. (n devient n-1)
  X(:,1,:)=X(:,2,:);
  V(:,1,:)=V(:,2,:);
  XS(:,1)=XS(:,2);
  VS(:,1)=VS(:,2);
  
end
fclose(fdata)
toc

@Gwend@l : tiens, c'est rigolo, je n'avais pas fais gaffe quand tu me l'avais montré, mais pourquoi tu as mis une vitesse verticale pour le Soleil ?

Deux raisons :

1) au début c'était pour avoir un 3 ème corp en déplacement et voir si mon modèle réussisait quand même à faire tourner la Terre.

2) plus physiquement, le soleil se déplace vraiment (comme tout dans l'univers). C'est une visualisation basique de ce mouvement.

plus physiquement, le soleil se déplace vraiment (comme tout dans l'univers).

Je pense que tu as un peu forcé la dose par rapport aux autres vitesses, ça explique peut être que la Lune parte en couille.

Comme j'ai du temps, je vais essayer une modélisation à 9 corps (8 planètes + Soleil) mobiles, ça peut être marrant. Dès que c'est prêt, je poste ça !

Je suis arrivé à quelque chose !

Je vais essayé d’être clair sur la méthode utilisée. C’est la première fois que je m’attaque à un système fortement non-linéaire, il y a peut être des erreurs de débutants dedans, hésitez pas à me le faire remarquer !

Je suis parti du PFD appliqué au corps $i$ dont la position $\mathbf x_i$ varie avec le temps selon : $\ddot{\mathbf x_i}=\sum_{j\neq i}G\dfrac{M_j}{d_{ij}^3}(\mathbf x_j-\mathbf x_i)$ Où $d_{ij}=||\mathbf x_j-\mathbf x_i||_2$ est la distance entre deux planètes.

J’ai ensuite utilisé $G$ comme échelle pour les $\mathrm{m\cdot s^{-2}\cdot kg^{-1}\cdot m^2}$, une masse $M$ quelconque pour les $\mathrm{kg}$ (en pratique celle du Soleil par exemple), et l’unité astronomique comme unité de longueur.

Ceci mène à l’échelle de temps $\sqrt{\dfrac{L^3}{GM}}$ soit environ $5\cdot 10^6\mathrm{s}$ dans le cas où on prend la masse du Soleil comme référence.

Les équations adimensionnées s’écrivent alors $\ddot{\mathbf x_i}=\sum_{j\neq i}\dfrac{M_j}{d_{ij}^3}(\mathbf x_j-\mathbf x_i)$

Pour résoudre ce système, je suis parti sur un schéma implicite. En utilisant les deux développements en série de Taylor : $x_i^n=x_i^{n+1}-\Delta t\partial_tx_i^{n+1}+\dfrac{\Delta t^2}{2}\partial^2_tx_i^{n+1}$ $x_i^{n-1}=x_i^{n+1}-2\Delta t\partial_tx_i^{n+1}+2\Delta t^2\partial^2_tx_i^{n+1}$

On obtient le schéma implicite suivant : $(1-\Delta t^2\partial^2_t)x_i^{n+1}=2x_i^n-x_i^{n-1}$

Pour calculer $x_i^{n+1}$ il faut donc inverser l’opérateur $\mathbf F:(1-\Delta t^2\partial^2_t)\bullet -2x_i^n+x_i^{n-1}$, autrement dit résoudre l’équation $\mathbf F(\mathbf x^{n+1})=0$

Comme cette équation est non-linéaire, je me sers de la méthode itérative de Newton suivante : $\mathbf J_{\mathbf F}(\mathbf x_k)(\mathbf x_{k+1}-\mathbf x_k)=-\mathbf F(\mathbf x_k)$ où $\mathbf J_{\mathbf F}(\mathbf x_k)$ est la matrice Jacobienne $J_{ij}=\dfrac{\partial F_i}{\partial x_j}$ évaluée en $\mathbf x_k$. L’expression analytique de la matrice Jacobienne est tout aussi connue que dégueulasse (j’avoue avoir la flemme de l’écrire ici ), il suffit d’appliquer bêtement la définition et ne pas se gourer en dérivant.

D’un point de vue pratique, j’ai codé en fortran le bordel dans le cas général (en utilisant LAPACK pour l’inversion de la matrice Jacobienne). On met autant de corps qu’on veut, le problème est résolu avec 3 dimensions spatiales.

Il suffit d’entrer le nombre de corps, leurs masses, et deux positions initiales séparées d’un pas de temps dans le fichier d’entrée.

Pour l’instant je n’ai testé qu’avec le Soleil et la Terre, et ça marche très bien. C’est même très rapide puisque 100000 pas de temps (donc en gros 16 ans avec un pas de temps adimensionné de $10^{-3}$) sont réalisés en grossièrement une demi-seconde. J’hébergerai le code quelque part (sur GitHub probablement), et j’essayerais de faire tourner ça avec 9 planètes en prenant une situation initiale réaliste.

Hésitez pas à me dire si c’est pas clair (ou si il y a une boulette).

Je veux bien voir le code !

Je veux bien voir le code !

Voici le repo git annoncé (le nbody.in étant celui utilisé pour un simulacre de système Terre+Soleil, la Terre s'approche un peu du Soleil au cours du temps). Hésitez pas à jouer avec le code.

Je vais me pencher sur le cas avec 8 planètes + Soleil (EDIT: solar.in sur le dépôt), et il faudra que je fasse des tests de performances (peut être que ça vaut le coup de paralléliser les calculs de $\mathbf F(\mathbf x_k)$ et de $\mathbf J_{\mathbf F}(\mathbf x_k)$).

Très jolie ! Spontanément, je ne vois pas d'erreur.

Avis très personnel, je trouve que c'est un peu compliqué pour un système physique aussi simple. Avec un pas de temps similaire (1e5 pas de 5e3 secondes, soit ~16 ans), la résolution des équations du mouvement non matriciel prend 6 secondes en python chez moi.¹

Si les performances ne suivent pas avec 8 planètes, plutôt que de tenter une parallélisation qui va te donner, au mieux, un rapport 10 sur une machine personnelle, je te recommande de négliger l’interaction des planètes légères entre elles. Le Soleil, c'est quelque chose comme 99% de la masse du système solaire, et Jupiter représente 2/3 de ce qu'il reste, donc l'influence de Neptune sur Venus ou Mercure est négligeable. Probablement même physiquement négligeable devant les corrections relativistes (je rappelle que la théorie de Newton ne décrit par correctement l'orbite de Mercure).

En tout cas, c'est très intéressant de voir deux approches assez différente de la résolution d'un même problème. Je zieuterai le code d'ici quelques jours.

Si les performances ne suivent pas avec 8 planètes, plutôt que de tenter une parallélisation qui va te donner, au mieux, un rapport 10 sur une machine personnelle, je te recommande de négliger l’interaction des planètes légères entre elles.

C'est pas bête, mais vu comme le machin est codé, ça risque d'être compliqué à mettre en place. Si on néglige les interactions entre les deux dernières planètes et les quatre premières, on met juste quelques termes à 0 dans la matrice jacobienne et on fait sauter quelques termes dans les sommes de $\mathbf F$. Mais je pense pas que ça change grand chose pour ce qui est de l'inversion de $\mathbf J$ ou du calcul de $\mathbf F$ (remplacer un petit calcul par un saut conditionnel, ça doit même coûter plus que ça ne rapporte je pense).

EDIT : d'ailleurs, en y repensant, si on vire les 4 premières planètes vis-à-vis des deux dernières, on fait 2 fois moins de calculs pour les deux dernières planètes et 1.5 fois moins pour les 4 premières. En idéalisant, on gratte un facteur 2 maximum. Parallelization wins again!

Avis très personnel, je trouve que c'est un peu compliqué pour un système physique aussi simple.

Pas si simple que ça. Le problème à n-corps fait partie des problèmes les plus difficiles à résoudre de la physique actuelle. Je placerais ça au même niveau que la mécanique des fluides.

Tu ne l'explicite pas, donc pour en être sûr : pour ton équation discrétisée, tu considères la quantité $(\mathbf x_j^n-\mathbf x_i^n)$, i.e. à l'instant précédent, n'est-ce pas ?

Je ne comprends pas ta notation pour l'opérateur $\mathbf F$. La méthode de Newton de base n'est pas garantie de converger, dans le cas général. As-tu vérifié qu'elle converge sans faute dans le cas de $\mathbf F$ ?

J'avoue ne pas être totalement convaincu par ton schéma d'intégration numérique ; intuitivement, il me semble plutôt pas mal, mais peut-être un petit peu lent et je verrais bien les planètes lentement s'approcher du Soleil et s'y crasher. Mais bon, l'intuition dans le domaine de l'analyse numérique… Pour moi, c'est à confronter expérimentalement à Euler semi-implicite.

Si tu prends uniquement Jupiter et le Soleil, tu passes de 9 corps en interaction (88 = 64 calculs) à 9 corps interagissant avec 2 (72 + 1 = 15 car pas d'auto-interaction). Bon OK, un rapport 4, soit peu ou prou ce que tu peux obtenir en parallélisant un quadri-cœur. Et en perdant un peu de précision (ce qui n'est pas le cas avec le parallélisme).

J'ai voulu vérifier ce que ça donnais avec le module linalg de python, qui utilise Lapack derrière. Je gagne 10%. Donc, avec ta méthode, cette « simplification physique » n'apporte pas de gain de calcul. Ce qui la rend vraiment meilleure que la mienne. Je n'ai pas le réflexe de faire des calculs matriciel, mais quand je vois ça, je me dit que c'est un tort !

Je dirais que ça signifie : l'opérateur qui à la fonction $X(t)$ associe $ (1−\Delta t^2 \partial_t^2) X(t) −2x^n_i + x^{n−1}_i$.

Tu ne l'explicite pas, donc pour en être sûr : pour ton équation discrétisée, tu considères la quantité $(\mathbf x_j^n-\mathbf x_i^n)$, i.e. à l'instant précédent, n'est-ce pas ?

C'est-à-dire ? Je résous l'équation modifiée directement pour $\mathbf x$ qui est un vecteur colonne qui contient toutes les positions des planètes. Lors du calcul de $\mathbf x^{n+1}$, $\mathbf x_j^{n+1}-\mathbf x_i^{n+1}$ n'intervient que dans le calcul de l'interaction de la planète $j$ sur la planète $i$. Comme c'est un schéma implicite, les opérateurs sont évidemment évalués en $n+1$, d'où la nécessité d'inverser un opérateur non-linéaire.

Je ne comprends pas ta notation pour l'opérateur $\mathbf F$.

Lorsque je veux calculer $\mathbf x$ à l'instant $n+1$, je résous en fait : $$\mathbf F(\mathbf x^{n+1})=\mathbf x^{n+1}-\Delta t^2\partial _t^2\mathbf x^{n+1}-2\mathbf x^n+\mathbf x^{n-1}=0$$

Introduire l'opérateur $\mathbf F$ ne sert qu'à alléger les notations suivantes.

La méthode de Newton de base n'est pas garantie de converger, dans le cas général. As-tu vérifié qu'elle converge sans faute dans le cas de $\mathbf F$ ?

Je pars de $\mathbf x^n$ lors de l'itération servant à calculer $\mathbf x^{n+1}$. Comme le pas de temps est suffisamment petit, $\mathbf x^n$ et $\mathbf x^{n+1}$ sont assez proches pour que Newton converge en quelques itérations (3 ou 4 en pratique avec un résidu de $10^{-30}$).

Oui, il est évident que rien ne garanti à priori que le schéma que je propose est pertinent. Le défi pour moi, c'était surtout d'arriver à faire un schéma implicite avec un opérateur non-linéaire à inverser.

Si tu prends uniquement Jupiter et le Soleil, tu passes de 9 corps en interaction (88 = 64 calculs) à 9 corps interagissant avec 2 (72 + 1 = 15 car pas d'auto-interaction). Bon OK, un rapport 4, soit peu ou prou ce que tu peux obtenir en parallélisant un quadri-cœur. Et en perdant un peu de précision (ce qui n'est pas le cas avec le parallélisme).

Ça me semble un peu osé de dire que Venus n’interagit pas sur la Terre par exemple. Tout comme ça me parait un peu osé de prendre en compte l'action de Jupiter sur Mars, mais pas celle de Mars sur Jupiter.

Vu que l'on fait des problèmes à N corps, pourquoi ne pas utiliser des algorithmes de dynamique moléculaire ? Genre l'algorithme de Verlet-Vitesse ou de Beeman ?

Ça me semble un peu osé de dire que Venus n’interagit pas sur la Terre par exemple.

Au pifomètre, le Soleil détermine le mouvement et Jupiter le 1er ordre des perturbations. Après, on peut toujours raffiner. En gros, Jupiter est 400 fois plus massive que Vénus, et situé à 5 UA de la Terre, contre (en gros) 1 UA pour Vénus. L'attraction gravitationnelle de Jupiter est donc environ 15 fois plus forte.

Tout comme ça me parait un peu osé de prendre en compte l'action de Jupiter sur Mars, mais pas celle de Mars sur Jupiter.

On ne peut pas négliger l'influence de la Terre sur toi, mais le contraire est vrai. L'influence de Mars sur Jupiter (en terme d'accélération) est de l'ordre de la masse de Mars sur la distance carrée, celle de Jupiter sur Mars est de l'ordre de la masse de Jupiter sur la distance carrée.

Vu que l'on fait des problèmes à N corps, pourquoi ne pas utiliser des algorithmes de dynamique moléculaire ?

C'est, en très gros, ce que fait l'algorithme proposé en page 1. C'est une version pas optimisée de Verlet (Euler). Après, ça marche !

L'attraction gravitationnelle de Jupiter est donc environ 15 fois plus forte.

Si on fait sauter grossièrement 10% des interactions, ça fait quand même une grosse erreur (surtout sur un système chaotique…).

L'influence de Mars sur Jupiter (en terme d'accélération) est de l'ordre de la masse de Mars sur la distance carrée, celle de Jupiter sur Mars est de l'ordre de la masse de Jupiter sur la distance carrée.

Mars fait en gros $3\cdot 10^{-4}$ fois la masse de Jupiter. Négliger un terme de l'ordre du pour mille (en prenant en compte les 7 autres planètes), c'est se tirer une balle dans le pied et ruiner tous les efforts fait sur la discrétisation du problème (qui entraine ici une erreur de l'ordre de $\Delta t^3\approx 10^{-9}$ avec un pas $10^{-3}$).