Bonjour à tous,
On parlait (là) de faire une simulation numérique de la gravitation pour simuler le mouvement de Rosetta. Je proposais (ici) une méthode. Je l'ai fait.
Caractéristiques techniques
- Langage : Python2/3.
- Dépendances : Scipy, Numpy et Matplotlib.
- Limitation : Résolution des EDO par la méthode d'Euler, qui est non stable !
- Licence : zlib, donc libre sans viralité. (ce n'est pas précisé dans les fichier d'ailleurs !)
Ça marche, c'est-à-dire qu'en donnant des valeurs réels pour le système Terre-Soleil, ça donne le bon résultat. C'est très limité, car le système est très sensible aux conditions initiales (et j'en ai bavé à cause de cela !), le tout avec une méthode pas stable… La version avec scipy.integrate.odeint devrait venir sous peu, mais comme je maitrise mal, je voulais avoir un truc qui marche avant (et comme ça marche, je le présente).
Utilisation
Le script s'utilise comme suit :
1 | python(3|2) gravite.py test.json |
Commentaires
Le fonctionnement interne est le suivant : la classe ObjetCeleste (fichier ObCe) contient un objet du système solaire avec son nom, sa masse, ses conditions initiales de position et vitesse et deux booléens attracteur/attiré (avec un seul "t" dans le code ! Yeah l'orthographe ! ) pour limiter les calculs (le soleil est fixe, les comètes n'attirent rien). Le tout est initialisé avec un fichier json lut par parseur (fichier json_parser). Le fichier gravite fait l'affichage et la résolution du système. La dérivation du système se fait dans la classe ObjetCeleste MAIS devra se faire dans une fonction annexe quand on passera à scipy.integrate.odeint. Le code est relativement bien commenté (sans vérification orthographique cependant). Aucune gestion des exception.
Du code de physicien. Mais du bon code de physiscien : l'entrée est via la console, et le nombre d'objet n'est pas défini dans le corps du programme. J'ai déjà fait bien pire…
Le code
Je rappelle : la méthode n'est pas stable et le système est très sensible aux conditions initiales. Les 4 fichiers sont ci-dessous :
ObCe.py
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 | #!/usr/bin/env python3 # -*-coding:Utf-8 -* """ Définit la classe ObjetCeleste qui represente un objet dans le système solaire. """ #Python2 special import. from __future__ import (division, absolute_import, with_statement, print_function, unicode_literals) #Scientific import. import numpy as np import scipy as sp class ObjetCeleste: """ Définit un objet céléste par : - ses coordonnées en espace ; (couple de flottant) - ses coordonnées en vitesse ; (couple de flottant) - son nom ; (chaine de caractère) - sa masse ; (flottant) - comment il se comporte vis-à-vis des autres. (couple de booléen) Possède les fonctions résolvant les équations différentiels du mouvement. """ def __init__(self, nom, masse, pos, vit, comportement): """Initialisation simple de l'objet.""" self.nom = nom self.pos = np.array(pos) self.vit = np.array(vit) self.masse = masse self.attracteur, self.atire = comportement #Atiré def gravite(self, objet): """Renvoie l'accélération dû à l'attraction d'un autre corps.""" G = 6.67e-11 module_acc = -G*objet.masse / sum((self.pos-objet.pos)**2) #array : carré terme à terme. #On projete dans le plan. Pour cela, on calcule l'angle theta entre objet.pos/self.pos et l'axe des X pris en self.pos. #Passage par l'arctan est le quart de cercle. arctan = sp.arctan( (objet.pos[1]-self.pos[1]) / (objet.pos[0]-self.pos[0]) ) if (arctan > 0 and self.pos[1] > 0) or (arctan < 0 and self.pos[1] < 0): #Cas facile : artan donne la bonne valeur de theta entre -pi/2 et pi/2 dans le demi-cercle droit. theta = arctan else: #Demi-cercle gauche ; atan va de -pi/2 (|y|>>|x|, y>0) à pi/2 (|y|>>|x|, y<0), #avec x<0 sur tout le demi-cercle et atan~0 pour |x|>>|y|. #D'où theta = atan + pi. Theta ira de -pi/2 à 3pi/2 sur tout le cercle. theta = arctan + sp.pi ##DEBUG print(sum((self.pos-objet.pos)**2)) acc_x = module_acc * sp.cos(theta) acc_y = module_acc * sp.sin(theta) return np.array([acc_x, acc_y]) def __str__(self): """Affichage.""" return self.nom+", actuellement en " + str(self.pos) + " va à la vitesse " + str(self.vit) + \ ". Sa masse est " + str(self.masse) + \ ". Est-il attracteur-attiré ? " + str(bool(self.attracteur)) + "," + str(bool(self.atire)) + "." def resout(self, objets_celestes, dt): """Résout les équations différentiels du mouvement et mets à jour position et vitesse. objets_celestes est la liste de tout les objets céléstes présents.""" acc = np.array([0.,0.]) if self.atire: for oc in objets_celestes: if oc.attracteur and oc is not self: res = self.gravite(oc) acc[0] += res[0] acc[1] += res[1] ##DEBUG print(res) self.pos += self.vit*dt self.vit += acc*dt |
gravite.py
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 | #!/usr/bin/env python3 # -*-coding:Utf-8 -* """ Entree-sortie. """ #Python2 special import. from __future__ import (division, absolute_import, with_statement, print_function, unicode_literals) import sys import numpy as np import scipy as sp import matplotlib.pyplot as plt from json_parser import parseur #TRY fichier = sys.argv[1] #Objets célèste est la liste de tous les objets, avec leur position initiales. objets_celestes = parseur(fichier) ##DEBUG print(objets_celestes[0], objets_celestes[1]) T = 10000 ; dt = 1e4 #Memoire stoque les position et vitesses des objets. memoire = np.zeros((T, 4*len(objets_celestes))) for i, objet in enumerate(objets_celestes): memoire[0, 4*i] = objet.pos[0] ; memoire[0, 4*i+1] = objet.pos[1] memoire[0, 4*i+2] = objet.vit[0] ; memoire[0, 4*i+3] = objet.vit[1] for t in range(1,T): for i,objet in enumerate(objets_celestes): objet.resout(objets_celestes, dt) for i, objet in enumerate(objets_celestes): memoire[t, 4*i] = objet.pos[0] ; memoire[t, 4*i+1] = objet.pos[1] memoire[t, 4*i+2] = objet.vit[0] ; memoire[t, 4*i+3] = objet.vit[1] ##DEBUG print(memoire) plt.plot(memoire[:,0],memoire[:,1]) plt.show() |
json_parser.py
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 | #!/usr/bin/env python3 # -*-coding:Utf-8 -* """ Contient la fonction parseur qui lit un fichier json et le convertit en objet python. """ #Python2 special import. from __future__ import (division, absolute_import, with_statement, print_function, unicode_literals) import json from ObCe import ObjetCeleste def parseur(fichier): """ Renvoie une liste d'objet céléstes correctement formaté en lisant le fichier json 'fichier' (chaine de caractère). """ with open(fichier) as fich: json_decode = json.load(fich) #Un petit try serai bien ici. objets_celestes = [] for objet in json_decode: objet_celeste = ObjetCeleste(nom=objet['nom'], masse=objet['masse'],\ pos=objet['pos'], vit=objet['vit'],\ comportement=objet['bool']) objets_celestes.append(objet_celeste) return objets_celestes |
test.json
1 2 3 4 | [ {"nom":"planete", "pos":[0.00001,1e11], "vit":[30e3,0.0], "masse":6e24, "bool":[0, 1]}, {"nom":"etoile", "pos":[0,0], "vit":[0,0], "masse":1e30, "bool":[1, 0]} ] |
Je reste dans le coin (pas ce soir, dodo) pour les questions, et proposez une version stable d'ici quelques temps.
Bonne soirée !
edit Arius : ajout du tag CdS