Bonjour tout le monde !

Je suis un peu nouveau ici donc c'est mon tout premier post, et histoire de bien commencer j'ai déjà un problème
Donc voilà j'explique la situation ! C'est un problème mélangeant à la fois informatique (Programmation) et physique ! J'espère que vous serez capable de m'aider puisque j'ai du mal concernant l'exploitation de mes résultats c'est-à-dire la dernière partie de ce devoir ! J'ai déjà fais mes programmes concernant l'informatique c'est juste l'exploitation qui me pose problème ! Donc je vous présente le problème :

L'objectif de ce microprojet est la résolution de l'équation différentielle du mouvement d'une sonde dans un champ gravitationnel planétaire par des méthodes numériques. On souhaite mettre en évidence l'impact du pas de discrétisation sur la qualité des résultats et sur le temps de calcul. On comparera les résultats obtenus avec ceux des fonctions déjà implémentées dans les bibliothèques numériques de Python.

1. Présentation du problème physique
On considère une sonde spatiale de masse m dans le champ gravitationnel de la planète Terre de masse M P =5,97 × 1024 kg et de rayon RT =6.6356∗106 m . On suppose cette sonde soumise à cette unique force centrale gravitationnelle. Le mouvement est donc plan et repéré par les coordonnées polaires ( r ,θ ) avec ⃗OM=ru⃗r
La mise en équation par la deuxième loi de Newton conduit à :
r ̇=dr dt
m⃗a=−GmMP u⃗ r2 r
  θ ̇=dθ
dt avec α=GMP
{r ̈ ( t ) = − α + r θ ̇ 2 r2
θ ̈ ( t ) = − 2 r ̇ θ ̇ r
On prendra comme conditions initiales r(t0)=R+h
θ(t0)=0
{r ̇(t0)=0 avec vi∈[0.8v0,0.9v0,v0,1.1v0,1.2v0] ,
θ ̇(t0)= vi R+h
(altitude de Landsat)
2. Résolution de l'équation du mouvement. 2.1. Par la méthode d'Euler
v0=7365m.s−1
et h=705.0km
  • Réaliser la discrétisation du problème ci-dessus qui permette la résolution approchée de l'équation différentielle par la méthode d'Euler.
• Implémenter en python un programme qui permette de résoudre l'équation différentielle du mouvement par la méthode d'Euler et ainsi de déterminer l'évolution temporelle de r(t) et θ(t) en fonction des conditions initiales.
Vous justifierez notamment le choix de la durée de l'étude
2.2. Grâce à des fonctions numériques de la bibliothèques scipy
En utilisant la fonction déjà implémentée « odeint » de la bibliothèque scipy.integrate de Python, écrire un programme qui permette de résoudre à nouveau l'équation différentielle du mouvement du pendule simple et de déterminer l'évolution temporelle de r(t) et θ(t) .
3. Analyse des résultats obtenus
On souhaite comparer les résultats obtenus par ces méthodes numériques avec une solution exacte qui est connue et vaut
p C2 (r20 θ ̇0)2 p r(θ)=1+ecosθ avec p=GM =GM e=r −1
PP0
• Dans un premier temps, la comparaison de ces deux méthodes avec la solution exacte se fera de manière graphique en superposant, sur un même graphique, les courbes obtenues pour les différentes conditions initiales. Vous pourrez ainsi illustrer l'influence du pas d'intégration sur la qualité de la résolution numérique. Vous préciserez bien pour chaque courbe présentée les conditions initiales et le pas d'intégration
• Dans un deuxième temps, afin de quantifier plus rigoureusement l'influence du pas d'intégration, réaliser une étude qui mette en évidence l'évolution du temps d'intégration et de l'erreur en fonction du pas d'intégration. Vous pourrez synthétiser vos résultats dans un tableau
4. Exploitation
A quelle condition la sonde rentre en collision avec la planète de rayon R ?

Concernant l'exploitation c'est la dernière phrase du haut

import matplotlib.pyplot as plt
import scipy.integrate as spi
import numpy as np
import time as time

def Euler(r0,rprime0,theta0,thetaprime0,a,b,N):
    Tinitial=time.time()    
    t=a
    r=r0
    rprime=rprime0
    theta=theta0
    thetaprime=thetaprime0
    h=(b-a)/N
    R=[r0]
    Rprime=[rprime0]
    T=[a]
    Theta=[theta0]
    Thetaprime=[thetaprime0]
    for i in range(N):
        r,rprime,theta,thetaprime,t=r+h*rprime,rprime+h*((-3.98428*10**14/r**2)+(r*thetaprime**2)),theta+h*thetaprime,thetaprime-(2*h*rprime*thetaprime)/r,t+h
        R.append(r)
        Rprime.append(rprime)
        Theta.append(theta)
        Thetaprime.append(thetaprime)
        T.append(t)
    plt.figure('méthode Euler')
    plt.subplot(211)
    plt.title('theta(t)')
    plt.plot(T,Theta)
    plt.subplot(212)
    plt.title('rayon(t)')
    plt.plot(T,R)
    plt.figure('r(theta)')
    plt.plot(Theta,R)
    plt.title('Méthode Euler')
    plt.ylabel('r(theta)')
    plt.xlabel('theta')
    Tfinal=time.time()    
    return Theta,R
R=np.array(Euler(7.34021*10**6,0,0,0.001003323978966297,0,86400,100000)[1])
Theta=np.array(Euler(7.34021*10**6,0,0,0.001003323978966297,0,86400,100000)[0])

def odep(y,t):
    r,theta,rprime,thetaprime=y
    dydt=[rprime,thetaprime,(-3.98428*10**14/r**2)+r*thetaprime**2,-2*rprime*thetaprime/r]
    return dydt
    
y0=[7.34021*10**6,0,0,0.001003323978966297]

T=np.linspace(0,86400,1000)
X=spi.odeint(odep,y0,T)

Rayon=X[:,0]
theta=X[:,1]
plt.figure('scipy')
plt.plot(theta,Rayon)
plt.show



#
Tethasol=[i/100 for i in range(10001)]

def solutionexacte(tetha):
    thetapoint=7365*1./(7.34021*10**6)
    p=(5.38787*10**13*thetapoint)**2/(3.98428*10**14)
    e=(p/(7.34021*10**6))-1
    return p/(1+e*np.cos(tetha))
Rsol=[solutionexacte(j) for j in Tethasol]
plt.figure('solution exacte v=v0')
plt.plot(Tethasol,Rsol)
plt.xlabel('tetha en rad')
plt.ylabel('rayon en metre')
#
#
#print('euler')
#print((R-Rsol)/Rsol)
#print('scipy')
#print((Rayon-Rsol)/Rsol)
#
#
#def moyennelistes(R):
#    m=0
#    for i in range(len(R)):
#        m=m+R[i]
#    return(m/len(R))
#    



plt.title("en vert scipy, en bleu euler, en rouge solution exacte")
plt.plot(theta,Rayon,"gs",Euler(7.34021*10**6,0,0,0.001003323978966297,0,86400,1000)[0],Euler(7.34021*10**6,0,0,0.001003323978966297,0,86400,1000)[1],"b",Tethasol,Rsol,"r--")

Voila ! Je ne peux malheureusement pas joindre les graphiques obtenus qui sont beaucoup trop nombreux en espérant que vous comprendrez !

Merci d'avance pour votre aide !