oa112023 14.41 TO N*S : Equations difrenteies ordinares TD N°5: Equations différentielles ordinaires Ahmed Ammar ([Link] at [Link]) Institut Préparatoire aux Etudes Scientifiques et Techniques, Université de Carthage. Mar 10, 2020 Table des matiéres Exercice 1: Pendule simple Exercice 2: Comparaison des schémas d’Buler explicite et implicite Exercice 3: Atterrissage d'un vaisseau spatial Exercice 1: Pendule simple On considére un pendule simple de masse m = 1 kg , de longueur ! avant a cause du champ de gravité de la Terre g = 9.8 m/s? [Link] iofFormationPythonHAPZ0 19-2020/pub/TOSITDS. html 1 m qui va osciller d’arriére enoa112023 14.41 TO N'S : Equations diférentlies ordinaires Y on Le pendule a 'équation du mouvement : 6 Fsin(a) () Pour les petites amplitudes d’oscillation, @ < 1, on peut faire l'approximation sin(@) ~ 0, on retrouve alors l’équation différentielle d'un oscillateur harmonique: (2) La solution exacte de cette équation est simplement: 4(t) =O cos(wot) (3) od wy = /g/l et nous avons supposé que le pendule partait du repos avec un déplacement initial 0 = 0.2 rad. ‘Nous allons transformer I'équation différentielle d’ordre 2 (Eq. (2)) en deux équations différentielles ordre 1 afin de pouvoir utiliser simplement la méthode d’Euler. En posant w(t) = 6 (t) la vitesse angulaire du pendule, on obtient le systéme de deux fonctions inconnues suivant : 4(t) = w(t) (4) a(t) = —wp A(t) 3) Pour résoudre ce systéme nous devons connaitre les deux conditions initiales suivantes (0) = 6 w(0) =0 a) Définir une fonction sol_exacte(t) qui renvoie la solution exacte de I'oscillateur harmonique donnée par I'équation (3). Tracer cette solution pour t € [0, 10] et pour un pas de At = 0.01 s. 1G Indication. * Utiliser la fonction [Link]() pour créer le vecteur temps t. * Utiliser la fonction [Link]() pour tracer sol_exacte(t) - '@ Solution. bitpsastrax. github iF ormationPython#AP20 19-2020/pub/TDS/[Link]| ansoart2029 14.64 TONS: Equations difrentels ordinates illateur harmonique est le suivant Le programme Python qui renvoie et trace la solution exacte de Io import numpy as np import [Link] as plt # SYSTEME: PENDULE STMPLE 8 = 9.8 # accélération de pesanteur [m/s*2] 1 = 1 # Longeur du pendule [m} dt = 0.01 # pas du temps [s] Tf = 10 # temps finale de La simulation [s] thetad = 0.2 # angle initiale [rad] omega® = [Link](g/1) # SOLUTION EXACTE def pendule_exacte(t): return thetad * [Link](omega® * t) t = [Link](@, Tf, dt) [Link](t, pendule_exacte(t), linewidth=2, ol exacte") pit. legend() [Link] ("Amplitude d'oscillation [rad]") [Link]("Tenps [s]") pit. title("Oscillateur Harmonique") [Link]() L'exécution de ce programme donne la figure suivante: Osciltaur Harmanique isa = un + AEA su, = (I+ ALA) ui (6) oau, = (%),a4= (9 — 1) cerestia matrice identité. —g/l 0 a(t) ¢) Calculer u (vo ) avec la méthode d’Euler explieite pour ¢ € (0, 10] et pour un pas d’intégration At = 001s. bitpsastrax. github iF ormationPython#AP20 19-2020/pub/TDS/[Link]| 33oa112023 14.41 TO N*S : Equations diférentoies ordinares ‘Tracer + Dans un méme graphique, la variation de l'amplitude d'oscillation 8 en fonction du temps t et le diagramme des phases (vitesse angulaire w en fonction de 6). + Dans un graphique 3D, la vitesse angulaire w et I'amplitude d'oscillation en fonction du temps ¢. Que remarquez-vous pour le résultat trouvé? 1G Indication. 1@ Solution. 4) Rappeler lexpression de la méthode d'Ewler implicite pour ce systéme. '@ Solution. up = (I— AtA)? - uy (7) 4 0 1 ovw, ( A ( ) et Festa matic ident on g/t 0 e) Calculer u = ( ote w(t At=0.01s. , ) avec la méthode d’Euler implicite pour t < [0, 10) et pour un pas d'integration ‘Tracer: + Dans un méme graphique, la variation de 'amplitude d'oscillation 9 en fonction du temps t et le diagramme des phases (vitesse angulaire w en fonction de 6). + Dans un graphique 3D, la vitesse angulaire w et 'amplitude d'oscillation @ en fonction du temps t. Que remarquez-vous pour le résultat trouvé? '@ Solution. bitpsastrax. github iF ormationPython#AP20 19-2020/pub/TDS/[Link]| ansoar112023 14.41 TO N*S : Equations diférentlles ordinaires #85 EULER IMPLICITE from [Link] import inv Timp = [Link](nsteps) Uimp = [Link]((2, nsteps)) Timp[@] = 0.0 UimpL:,0] = ue # ITERATION for k in range(nsteps-1): Tamptiea] = Tinptkl + a ] = [Link](inv([Link](2) - dt * A), Uimp[:,k]) pit. Figure(Figsize=(10,5)) # PLOT POSITION vs TEMPS [Link](1,2,1) [Link](Timp,Uimp[2,:], linewidth=2) [Link]("Temps [s]") [Link] (“Amplitude d'oscillation [rad]") [Link]("Oscillateur Harmonique (Euler implicite)") # DIAGRAMME DE PHASE [Link](1,2,2) [Link](Uimp[®,:],Uimp[1,:], linewidth=2) [Link] (“Amplitude d'oscillation (rad]") [Link] ("Vitesse angulaire [rad/s]" [Link]("Espace des phases (Euler implicite)") pit. savefig("Pendule_culerimp1D.png"); [Link]("Pendule_culerimp10. pdf") pit. show() # DIAGRAMME DE PHASE 30 [Link]() ax = [Link](projection="34") [Link](Timp, Uimp[@, :],Uimp(1, ax.set_xlabel("Tenps [s]") ax.set_ylabel("Anplitude d'oscillation [rad]") ax.set_zlabel("Vitesse angulaire [rad/s]") ax.set_title("Espace des phases (Euler implicite)") [Link]("Pendule_Eulerimp3D.png"); [Link]("Pendule_Eulerimp3D. pdf") [Link]() » linewidth=2) Liexécution de ce programme donne les figures suivantes: bitpsastrax. github iF ormationPython#AP20 19-2020/pub/TDS/[Link]| 53TO N'S : Equations diférentolies oinaires oart2029 14.64 Oscilateur Harmonigue (Euler implicite) Espace des phases (Euler implicte) 020 os ois ith % ow s : 3 £ cos £ a £ ay | io 2-008 $2 &-010 nos -oas 08 -0.20 of 4 6 8 2. 21 00 01 @2 Temps] Amplitude doscilation fad] et la figure en 3D: Espace des phases (Euler implicite) 4 ze ov o2 & oo & alias -oa8 ~067 Remarque. Dans le cas d'intégration avec la méthode d'Euler implicite nous avons un probléme de diminution d’amplitude dans le cas d’un oscillateur libre non amorti, Plus le temps de simulation est long, plus 'amplitude diminue, ce qui n'est pas ce que nous attendons de 'évolution du systéme dans le temps. nipsifastrax github jo/FormationPythonHAP2019-2020/pub/[Link] ensoastti2028 14:44 TONS : Equations diferentes ordnaires f) Tracer dans un méme graphique pour t € [0, 10] et avec un pas At = 0.01 s: * sol_exacte(t) calculée dans a). + 6(¢) calculée dans ¢) par la méthode d'Euler explicite. * 6(¢) calculée dans e) par la méthode d’Euler implicite. Que remarquez-vous si nous modifions la valeur du pas d'intégration par At = 0.001 s? Expliquer le résultat trouvé, 1@ Solution. ‘#85 TLLUSTRATION [Link]() [Link](t, pendule_exacte(t), linewidth=2, label: [Link](t, Uexp[@,:], linewidth=2, linestyl [Link](t, Uimp[@,:], linewidth=2, linestyl pit. legend() [Link] ("Temps [5] [Link] (“Amplitude d'oscillation [rad]") [Link]("Oscillateur Harmonique avec "+ r"$\Delta t =$"+str(dt)) pit. savefig("Pendule_illustration.png"); [Link]("Pendule_: [Link]() ol exact") ", label="Euler explicite") . "Euler implicite") Pour At = 0.01, l'exécution du code donne la figure suivante: Oscilateur Harmonique avec At = 0.02 wien 23) —— euler expen i == fulerimpicte 7, Aa 5 Temps (5) Pour At = 0.001, l'exécution du code donne la figure suivante: bitpsastrax. github iF ormationPython#AP20 19-2020/pub/TDS/[Link]| [Link]") 73oa112023 14.41 TO N*S : Equations diférentoies ordinares Oscilateur Harmonique avec At = 0.001 oa Eo : silent % 00 tule eve g fuerte Bo 02 ¥ Temps (5) Remarque. Les deux méthodes d'Euler, explicite et implicite, posent un probléme fondamental avec ses amplitudes eroissantes et décroissantes pour le cas d'oscillateur libre non amorti. Un trés petit At est nécessaire pour obtenir des résultats satisfaisants. Plus la simulation est longue, plus At doit étre petit. Exercice 2: Comparaison des schémas d’Euler explicite et implicite On considére le probléme de Cauchy: #0 =1- e tER, 2(0) = 2 (8) On rappelle que la solution exaete de ce probleme est donnée par: At) =" (uae * (9) a) Définir une fonetion sol_exacte(t, mu, 28) qui renvoie la solution exacte donnée par I'équation (9). Tracer sur un méme graphique pour yu = 1 et zo € {0, 1, 2} ces solutions. Soit ¢ < [0, 2] et pour un pas de At = 0.15. 1@ Solution. 1b) Méme questions pour pr = 0.05 et zo € {0, 1,2} 1@ Solution. On suppose dans cette question que 1 = 0.05 et que zy = 2. ©) Rappeler lexpression de la méthode dEuler explicite pour ce probleme. Calculer 2(t) avec la méthode d’Euler explicite pour t < [0,2] et pour un pas d'intégration At = 0.1. 1@ Solution. 4) Montrer que l'expression de la méthode d'Zuler implicite est: Itpsiastraxthubo/FormationyhonAPZ019-2020/pub/[Link] ansoa112023 14.41 TO N*S : Equations diférentlles ordinaires nt At sna = En 01,2001 1448 Calculer 2(t) avec la méthode d’Euler implicite pour t € [0, 2] et pour un pas dintégration At = 0.1 s. ‘Solution, €) Tracer dans un méme graphique pour t € [0, 2] et avec des pas d'intégration At = 0.5, 0.1, 0.05, 0.01, 0.005 s * La solution exacte: sol_exacte(t, 0.05, 2) + 2(t) calculée par la méthode d’Euler expli + 2(t) caleulée par la méthode d’Euler explicite. te. Que remarquez-vous pour les résultats trouvés? Quelle est la méthode la plus proche de la solution exacte? '@ Solution. Exercice 3: Atterrissage d'un vaisseau spatial Un vaisseau spatial est lancé & 'altitude H = 772 km au-dessus du niveau de la mer avec la vitesse v9 = 6700 m/s dans la direction indiquée sur la figure ci-dessus. Les équations différenticlles décrivant Je mouvement du vaisseau spatial sont: oti ret 6 sont les coordonnées polaires du vaisseau spatial. Les constantes impliquées dans le mouvement sont: * G = 6.672 x 10°"! m’kg-1s~? = constante gravitationnelle universelle. * Mr = 5.9742 x 1074 kg = masse de la terre. + Rr = 6378.14 km = rayon de la terre au niveau de la mer. bitpsastrax. github iF ormationPython#AP20 19-2020/pub/TDS/[Link]| onsoart2029 14.64 TON'S : Equations diférentolis ordinakos a) Dériver les équations différentielles du premier ordre et les conditions initiales de la forme ¥ = F(t,y),y(0) = 1@ Solution, ‘Nous avons: GMr = (6.672 x 10°!) (5.9742 x 104) = 3.9860 x 10'4 m/s? yo r vw r y we a ys 6 Les équations équivalentes du premier ordre deviennent: vo n 2 14) iy yoy? — 3.9860 x 10" /y’ Fity) = |% $ ° Yn Ys vs —2yrys/yo et les conditions initiales sont: (0) = Rr + H = (6378.14 + 772) x 10? = 7.15014 x 10° m (0) (0) =0 (0) = up/r(0) = (6700)/ (7.15014 x 10°) = 0.937045 x 10~* rad/s Ainsi, 7.15014 x 10° 0 y 0 0.937045 x 10-* b) Utiliser la méthode Runge-Kutta du quatriéme ordre (RK4) pour intégrer les équations depuis le lancement jusqu’a ce que le vaisseau spatial touche la terre. Déterminez 0 au site d'impact. 1G Solution. Le code Python est le suivant: bitpsastrax. github iF ormationPython#AP20 19-2020/pub/TDS/[Link]| sonsoar112023 14.41 TO N*S : Equations diférentlles ordinaires ## NOM DU PROGRAMME: #% IMPORTATION import numpy as np import [Link] as plt def RKa(F,t,y,h): ke = hFF(t,y) Ki = htF(t + h/2.0, y + K0/2.0) RKAEX3. py K2 = h*F(t + h/2.0, y + K1/2.0) K3 = hF(t +h, y + K2) return (KO + 2.0%K1 + 2.0%K2 + K3)/6.0 def F(t,y): F = [Link](4) Fle] = y{2) F[4] = y[o}*(y[3]**2) - 3.9860e24/(y[o]**2) Fl2] = yf3] F[3] = -2.0+y[1}*y[3]/yLe] return F t-0 tstop = 1200. h = 50.0 y = [Link]((7.15014e6, 0.0, 8.0, .937045e-3]) T=] y= [Link](t) [Link](y) while t < tStop: y = y + RKA(F,tyy,h) teteh [Link](t) [Link](y) T,Y = [Link](T),[Link](Y) print ("T Yre] y[a] Y2] Y[3]") for i in range(len(T)): print(T[i], "{:2.4e}".format(V[i,0]), "{:2.4e}".format(¥[i,1]), "{:2.4e}".format(¥[i,2]), "{:2.4e}".format (Yi, 3])) [Link](figsize=(8,5)) [Link](T, Y[:,@], '-0',1w = 2) [Link] ("temps [s]") [Link]("r [m]") [Link](y=6.3781466, color= # plot point d’impact [Link] (1034. 184, 6.3781466, 'rX") pt. savefig("[Link]"); [Link]("[Link]") [Link]() La sortie du programme est : bitpsastrax. github iF ormationPython#AP20 19-2020/pub/TDS/[Link]| 103oar112023 14.41 TO N*S : Equations diférentlles ordinaires T Y[2] you) y[2] YL3] 2.0 7.1501e+06 0.0000¢+00 0.0000c+00 9.37040-04 50.0 7.1482e106 -7.5908e101 4.6861e-02 9.3754e-04 100.0 7.1426e+06 -1.5173e+02 9.3771¢-02 9.3904e-04 150.0 7.1331e+06 -2.2738e+02 1.4078e-01 9.4153e-04 200.0 7.1198e+@6 -3.0276e+02 1.8794e-01 9.4504e-04 250.0 7,1028e+06 -3.7779e+02 2.353@e-01 9.4958e-04 300.0 7.0820e+06 -4.5236e+02 2.8292e-01 9.5515e-04 350.0 7.0576e+06 -5.2638e+02 3.3084e-01 9.6179e-04 400.0 7.0294e+06 -5.9973e+2 3.7911e-01 9.6951e-04 450.0 6.9976e+06 -6.7229e+02 4.2780e-01 9.7834e-04 500.0 6.9622e+06 -7.4393¢+02 4.7697e-01 9.8832e-04 550.0 6.9232e+06 -8.1452e+02 5.2666e-01 9.9947e-04 600.0 6.8808e+06 -8.8389e+02 5.7693e-01 1.0118e-03 650.0 6.8349e+06 -9.5189e+02 6.2786e-21 1.0255e-03 700.0 6.7856e+06 -1.0183e+03 6.7950e-01 1.0404e-03 750.0 6.7331e+06 -1.0830e+03 7.3193e-01 1.0567e-03 800.0 6.6773e+06 -1.1456e+03 7.8520e-01 1.0744e-03 850.0 6.6185e+06 -1.2060c+03 8.39400-01 1.09360-03 900.0 6.5568e+06 -1.2639e+03 8.9459e-01 1.1143¢-03 950.0 6.4922e+06 -1.3189e+03 9.5085e-01 1.1366e-03 1000.0 6.4250e+06 -1.3708e+03 1.0083e+00 1.1605e-03 1050.0 6.3552e+86 -1.4191e+03 1.0669e+00 1.1861e-03 1100.0 6.2831e+06 -1.4634e+03 1.1269e+00 1.2135e-03 1150.0 6.2089e+06 -1.5034e+03 1.1883e+00 1.2427e-03 1200.0 6.1329e+06 -1.5384e+03 1.2512e+00 1.2737e-03 Le vaisseau spatial frappe la terre lorsque r est égal A Rr = 6.37814 x 10° m. Cela se produit entre + = 1000 et 1050 s. Laissant 1000 + At étre le moment de l'impact, nous pouvons écrire: (1000 + At) = Rr En développant r dans une série Taylor a deux termes, nous obtenons: 7(1000) + #(1000)At = Rr 6.4250 x 10° — 1.3708 x 10°At = 6378.14 x 10* a partir duquel: At = 34.184 s La coordonnée @ du site d'impact peut étre estimée de maniére similaire. En utilisant 4 nouveau deux termes de la série Taylor, nous avons: (1000 + At) = 0(1000) + 6(1000)At = .0083 + 1.1605 x 107% x (34.184) bitpsastrax. github iF ormationPython#AP20 19-2020/pub/TDS/[Link]| 1.0480 rad = 60.00° sansoa112023 14.41 TO N‘S : Equations diférentolles ordinares nipsifastrax github jo/FormationPythonHAP2019-2020/pub/[Link] 133