Partage
  • Partager sur Facebook
  • Partager sur Twitter

Orbite terrestre Python

Simulation d'orbite

    28 juin 2020 à 17:28:30

    Bonjour,
    Je souhaite modéliser l'orbite terrestre sur python. J'ai fais une simulation sur python; le problème est que la terre dévie un peu de son orbite et au bout de quelque période la terre, l'orbite se décale d'une 100 de millions de kilomètres (ça fais bcp :) )
    Pour la simulation je n'ai utilisé que la mécanique newtonienne. Je n'ai pas pris en compte les autres planètes. Le soleil est fixe. J'ai utilisé la méthode d'Euler pour intégrer.
    Je précise que je suis d'un niveau d'un bon élève de prépa math sup en python.
    Voici mon code, si vous pouviez jeter un œil pour me dire comment le corriger merci 
    # initialisation des constantes
    #ra=rayon aphélie
    #rb=rayon périphélie
    
    ra=152093407000
    rb=147093407000
    a=(ra+rb)/2
    
    #masse de la terre
    mt=5.92*10**24
    
    #masse du soleil
    ms=1.9884*10**30
    
    #constante gravitationnelle
    G=6.6742*10**(-11)
    
    va=28851 #vitesse aphélie
    
    ax=[]
    ay=[]
    az=[]
    
    vx=[0]
    vy=[va]
    
    
    x=[ra]
    y=[0]
    
    
    
    d=1
    
    F=G*ms/(d**2)
    
    theta=0
    
    
    
    import numpy as np
    import matplotlib.pyplot as plt
    
    t=np.linspace(0,365*24*3600,365*24)
    dt=1/10
    for k in range(300000000):
        #print(k)
    
    
    
        d=(x[k]**2+y[k]**2)**0.5
    
        F=G*ms/(d**2)
    
        
    
        theta=np.abs(np.arctan(y[k]/x[k]))
     
    
        if x[k]<=0:
            if y[k]>=0:
                ax.append(F*np.cos(theta))
                ay.append(-F*np.sin(theta))
    
            if y[k]<0:
                ax.append(F*np.cos(theta))
                ay.append(F*np.sin(theta))
        if x[k]>0:
            if y[k]>=0:
                ax.append(-F*np.cos(theta))
                ay.append(-F*np.sin(theta))
            if y[k]<0:
                ax.append(-F*np.cos(theta))
                ay.append(F*np.sin(theta))
    
    
        
    
    
        vx.append(vx[k]+dt*ax[k])
        vy.append(vy[k]+dt*ay[k])
    
    
    
        x.append(x[k]+dt*vx[k])
        y.append(y[k]+dt*vy[k])
    
    
    
        plt.scatter(0,0,s=50,c='black',marker='+')
        plt.axis([-2*10**11,2*10**11,-2*10**11,2*10**11])
        plt.plot(x,y)
    
        plt.pause(0.01)
    plt.show()
    
    
    
    
    Les membres obtiennent plus de réponses que les utilisateurs anonymes Les membres obtiennent plus de réponses que les utilisateurs anonymes
    Gratuit depuis 1999Gratuit depuis 1999
    Mettez en avant votre expertise en aidant les autres&nbsp;! Mettez en avant votre expertise en aidant les autres !
    Navigation plus rapide sans aucune publicité&nbsp;!

    Les informations recueillies sont destinées à CCM BENCHMARK GROUP pour vous assurer l'envoi de votre newsletter.

    Elles seront également utilisées sous réserve des options souscrites, à des fins de ciblage publicitaire.

    Vous bénéficiez d’un droit d’accès et de rectification de vos données personnelles, ainsi que celui d’en demander l’effacement dans les limites prévues par la loi.

    Vous pouvez également à tout moment revoir vos options en matière de ciblage. En savoir plus sur notre politique de confidentialité.

    <iframe id="utif_ba_right_cdf53271-0d20-4485-a28a-6cd97f6b0f05" style="border: 0px; border-image: none;" name="utif_ba_right_cdf53271-0d20-4485-a28a-6cd97f6b0f05" frameborder="0" marginwidth="0" marginheight="0" scrolling="no" width="300" height="600" />

    -
    Edité par HadrienLevechin 28 juin 2020 à 17:29:14

    • Partager sur Facebook
    • Partager sur Twitter
      29 juin 2020 à 1:51:36

      Salut,
      Tu veux la réponse d'un astronome? Je n'en suis pas un ...
      À première vue (je n'ai pas vraiment analysé) tu fais tes calculs au mètre près.
      Ce que je pense est que tu calcules la position suivante en partant de la position courante (ce qui devrait être logique).
      Il y a ta variable dt qui vaut 1/10. C'est un dixième de quoi?
      Je pense qu'il faut toujours calculer la position par rapport à une position de référence, pour éviter la propagation des erreurs.
      La ligne suivante me semble une cause de propagation des erreurs. Je me trompe peut-être.
          theta=np.abs(np.arctan(y[k]/x[k]))
      Après deux tours, tu auras tourné de 4*PI radians. Tu pourrais soustraire 2*PI, mais ce n'est pas la solution.
      Tu dois te rendre compte que tu es revenu au point de départ et ramener ton angle à 0.
      Ne fonctionne pas en fonction des distances linéaires, mais des angles qui seront des fraction de tour.
      Par exemple, tu peux diviser un tour en 1000 ou 10000 parties, et tu comptes les parties.
      Je pense que tu devrais avoir de meilleurs résultats.
      Bonne chance.

      -
      Edité par PierrotLeFou 29 juin 2020 à 6:59:55

      • Partager sur Facebook
      • Partager sur Twitter

      Le Tout est souvent plus grand que la somme de ses parties.

        29 juin 2020 à 10:50:00

        Un sujet intéressant !

        Tu dis que, au bout de quelques périodes, ça donne un très grand écart. Qu'entends-tu par « quelques » ? Deux ou trois ? Un millier ?

        Pour l'instant je ne vois pas d'erreur visible (j'ai quand même un petit doute sur le calcul, parce que je ne vois où il y a une intégration, mais je ne me souviens pas de la formule donc ça ne veut rien dire), mais j'ai quelques petites remarques :

        1)

        ra=152093407000
        rb=147093407000

        Ces valeurs semblent calculées ainsi :

        149 593 407 km +/- 2 500 000 km

        donc à partir d'un demi-grand axe précis au km près, et d'une valeur connue à un demi millions de km près ! Tu devrais utiliser des valeurs plus précises, en tout cas plus cohérentes.

        2)

        d=(x[k]**2+y[k]**2)**0.5 # on calcule la racine carrée de x² + y²
        F=G*ms/(d**2)            # et on la remet au carré ! gaspillage...

        Si le programme doit effectuer beaucoup de calculs (et je crois que c'est le cas), voici deux lignes qui devraient à mon avis être optimisées. Calculer une racine carrée est toujours bien plus long que juste additionner ou multiplier, et là on pourrait très bien faire :

        d2=x[k]**2+y[k]**2       # on ne calcule pas la racine carrée de x² + y²
        F=G*ms/d2                # vu que c'est le carré de d qui est utilisé (ça tombe bien)

        3) Au fait, ce que tu nommes F n'est pas la force mais l'accélération. C'est trompeur de l'appeler F ! (Mais tu l'as utilisé correctement dans les calculs plus loin, dont tu le savais). (Du coup la masse de la Terre ne sert jamais ?)

        4)

        theta=np.abs(np.arctan(y[k]/x[k]))

        Tu devrais utiliser la fonction 'atan2' qui gère les quadrants.

        5) Pourquoi stocker l'accélération, la vitesse et la position dans une liste ? Il me semble qu'on peut utiliser une seule variable pour chaque donnée, par exemple :

        vx += dt*ax
        x  += dt*vx
        

        6) À quoi sert t (ligne 44) ?

        -
        Edité par robun 29 juin 2020 à 11:17:58

        • Partager sur Facebook
        • Partager sur Twitter

        Orbite terrestre Python

        × Après avoir cliqué sur "Répondre" vous serez invité à vous connecter pour que votre message soit publié.
        × Attention, ce sujet est très ancien. Le déterrer n'est pas forcément approprié. Nous te conseillons de créer un nouveau sujet pour poser ta question.
        • Editeur
        • Markdown