Méthodes informatiques pour physiciens

introduction à C++ et
résolution de problèmes de physique par ordinateur
Corrigé 10

Professeur : Alessandro Bravar

[Link]@[Link]

Université de Genève
Section de Physique

Semestre de printemps 2015

Références :

M-Y. Bachmann, H. Catin, P. Epiney et al.

Méthodes numériques

A. Quateroni, F. Saleri et P. Gervasio

Calcul scientifique

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery

Numerical Recipes

[Link] :
pour les notes du cours, les exercices et les corrigés
10.1 Problèmes
1. Pendule Euler
Voir le programme Pendule [Link] developpé pendant le cours.
Soit L la longueur du pendule, ω sa vitesse angulaire et ϑ l’angle d’oscillation. L’équation
exacte du mouvement s’écrit :
d2 ϑ(t) g
= − sin ϑ(t) (1)
dt2 L
Si l’on considère les oscillations suffisement petite (ϑ < 5◦ ), l’on peut faire l’approxima-
tion :
d2 ϑ(t) g
2
≈ − ϑ(t) (2)
dt L
qui a comme solution :
ϑ(t) = Acos(ω0 t + φ0 ) (3)
avec : s
2π L
ω0 = et T0 = 2π (4)
T0 g
A et φ0 sont déterminés à partir des conditions initiales.
Dans ce programme on calcule les oscillations dans l’approximation des petites oscillations
en utilisant la méthode d’Euler (avec dans les deux cas : ϑ0 = 5◦ et ω0 = 0). Les valeurs de
l’angle et de la vitesse angulaire (en considérant l’approximation des petites oscillations)
sont données par :

ϑn+1 = ϑn + ωn ∆t
g
ωn+1 = ωn + an avec ici an = − ϑn = ϑn
L
La méthode de Runge étant d’ordre supérieur, elle donnera une meilleure approximation
de la solution exacte que ne le fait la méthode d’Euler.
2. Pendule Runge
Voir le programme Pendule [Link].
Ce programme fait également le calcul de la trajectoire d’un pendule avec l’approxima-
tion des petites oscillations mais en utilisant la méthode de Runge. La méthode de Runge
consiste à calculer la vitesse au milieu de l’intervalle ∆t et la position au début de l’in-
tervalle. Pour démarrer l’itération il faut redéfinir la vitesse initiale ω0 dans t − ∆t/2
selon :
g
ω0 = ω0 − a0 ∆t/2 a0 étant l’accélération a0 = − ϑ0
L
qui nous permet d’obtenir

ω1 = ω(∆t/2) = ω0 + a0 ∆t

à partir de la nouvelle vitesse initiale ω0 . Et on calcule ensuite l’accélération, la vitesse et

la position à chaque itération selon l’algorithme :
g g
an = − ϑn avec ici =1
L L
ωn+1 = ωn + an ∆t
ϑn+1 = ϑn + ωn+1 ∆t

2
 Figure 1 – Comparaison entre les méthodes d’Euler et de Runge pour un pendule.

La figure 1 montre en rouge la solution obtenue avec la méthode d’Euler, en bleu avec
la méthode de Runge et en vert la solution exacte, qui se superpose parfaitement à la
solution de Runge. La méthode de Runge est d’ordre supérieure à celle d’Euler et donne
donc un bien meilleure approximation de la solution exacte.

3. Oscillateur anharmonique
L’oscillateur anharmonique est défini comme un oscillateur avec un force non linéaire, son
équation de mouvement est la suivante :
x α
ẍ = −k f (x) = −k |x|
|x|

L’utilisateur doit en premier lieu entrer la valeurs des constantes k et α.

On considére l’élongation initiale x0 = 10 et la vitesse initiale v0 = 0.
Dans la méthode de Runge pour démarrer l’itération il faut d’abord, calculer v1 = v(∆t/2)
x0
v1 = v0 + a0 ∆t/2 avec a0 = −k |x0 |α
|x0 |

ou redéfinir la vitesse initiale comme ci-dessous :

v0 = v0 − a0 ∆t/2

Ensuite on calcule l’accélération, la vitesse et la position à chaque itération :

xn
an = −k |xn |α
|x|
vn+1 = xn + an ∆t
xn+1 = xn + vn+1 ∆t

La figure 2 montre les solutions pour valeurs de α = 0 en bleu, = 1 en rouge, et = 2

en vert. Le coefficient k/m = 1 dans cet exemple. On constate qu’il y a une certaine
periodicité et que l’amplitude maximale ne change que peu dans le temps d’intégration
consideré.
Osc [Link]

3
 Figure 2 – Solutions avec différents valeurs du coefficient anharmonique.

1 // o s c i l l a t e u r anharmonique
2 #include <i o s t r e a m >
3 #include <cmath>
4 #include ” d i s l i n . h”
5
6 using namespace s t d ;
7
8 int main ( ) {
9 c o u t << ” Etude de l ’ o s c i l l a t e u r anharmonique ” << e n d l ;
10 c o u t << ” En trez l a c o n s t a n t e k/m de l ’ o s c i l l a t e u r : ” ;
11 double k = 1 . ;
12 c i n >> k ;
13 c o u t << ” En trez l ’ e x p o s a n t a l p h a ( a l p h a = 1 −> o s c . harmonique ) : ” ;
14 double a l p h a = 1 . ;
15 c i n >> a l p h a ;
16 c o u t << ” En trez l ’ e l o n g a t i o n i n i t i a l e : ” ;
17 double a n g l e 0 = 1 0 . ;
18 c i n >> a n g l e 0 ;
19
20 // i n i t i a l i s a t i o n de DISLIN
21 m e t a f l ( ”XWIN” ) ; //XWIN ou PDF
22 disini () ;
23 name ( ” temps ( s ) ” , ”X” ) ;
24 name ( ” e l o n g a t i o n (m) ” , ”Y” ) ;
25 t i t l i n ( ” O s c i l l a t e u r anharmonique ” , 1 ) ;
26 g r a f ( 0 . , 1 0 0 . , 0 . , 2 0 . , − 2 ∗ a n g l e 0 , 2 ∗ a n g l e 0 , −2∗ a n g l e 0 , a n g l e 0 ) ;
27 t i t l e () ;
28
29 const int STEPS = 1 0 0 0 0 ;
30 double dt = 0 . 0 1 ;
31 double a n g l e [ STEPS ] , omega [ STEPS ] , time [ STEPS ] ;
32 angle [ 0 ] = angle0 ; // e l o n g a t i o n i n i t i a l
33 omega [ 0 ] = 0 . ; // v i t e s s e i n i t i a l e
34 time [ 0 ] = 0 . ;
35
36 // methode de Runge
37 // c o n d i t i o n s i n i t i a l e s
38 double a c c ;
39 a c c = − a n g l e [ 0 ] / abs ( a n g l e [ 0 ] ) ∗k ∗ pow ( abs ( a n g l e [ 0 ] ) , a l p h a ) ;

4
Figure 3 – Solutions du pendule physique avec longueur L = 40 m, γ = 0.04 s−1 , F = 0.2 s−2 ,
Ω = 1.0 rad/s et angle initial 5◦ .

40 omega [ 0 ] = omega [ 0 ] − a c c ∗ dt / 2 . ; // methode de Runge !

41
42 f o r ( int i =1; i <STEPS ; i ++) {
43 a c c = − a n g l e [ i −1]/ abs ( a n g l e [ i −1]) ∗k ∗ pow ( abs ( a n g l e [ i −1]) , a l p h a ) ;
44 omega [ i ] = omega [ i −1] + a c c ∗ dt ;
45 a n g l e [ i ] = a n g l e [ i −1] + omega [ i ] ∗ dt ;
46 time [ i ] = i ∗ dt ;
47 }
48
49 // d e s s i n de l a t r a j e c t o i r e
50 c o l o r ( ”RED” ) ;
51 c u r v e ( time , a n g l e , STEPS) ;
52 disfin () ; // f i n de DISLIN
53
54 return 0 ;
55 }

4. Pendule physique
Le pendule physique subit en plus de la gravité une force d’amortissement et une force
d’entrainement. Dans le programme la solution présentée utilise la méthode de Runge.
On calcule d’abord v1 = v(∆t/2). Dans l’algorithme de Runge on calcule l’accélération en
fonction de l’ancienne vitesse et position, ensuite la nouvelle vitesse en fonction de cette
accélération, et enfin la position en fonction de la nouvelle vitesse :
g
a = − sin(ϑi−1 ) − γ ωi−1 + F sin(Ω ti−1 )
L
ωi = ωi−1 + a ∆t
ϑi = ϑi−1 + ωi−1 ∆t

Les equations ci-dessus sont formulées pour l’angle en radians (p. ex. l’accélération est
donnée en rad/s2 ). Les fonctions trigonométriques de la bibliothèque C++ prennent comme
arguments les angles en radians. La trajectoire du pendule physique est montrée dans la
figure 3.
Pendule [Link]

5
 1 // p e n d u l e p h y s i q u e a v e c l a methode de Runge
2 #include <i o s t r e a m >
3 #include <cmath>
4 #include ” d i s l i n . h”
5
6 using namespace s t d ;
7
8 int main ( ) {
9 const int STEPS = 5 0 0 0 0 ;
10 double dt = 0 . 0 0 2 ;
11 double a n g l e [ STEPS ] , omega [ STEPS ] , time [ STEPS ] ;
12 a n g l e [ 0 ] = 5 . ∗ M PI / 1 8 0 . ; // a n g l e i n i t i a l en r a d i a n s
13 omega [ 0 ] = 0 . ; // v i t e s s e a n g u l a i r e i n i t i a l e
14 time [ 0 ] = 0 . ;
15
16 // c o n s t a n t e s du p e n d u l e p h y s i q u e
17 double g = 9 . 8 1 ;
18 double L = 4 0 . ;
19 double gamma = 0 . 0 4 ;
20 double f o r c e = 0 . 2 ;
21 double f r e q = 2 . ;
22
23 double accG , accAm , accEn , a c c ;
24 // v i t e s s e i n i t i a l e pour l a methode de Runge
25 accG = −g /L ∗ s i n ( a n g l e [ 0 ] ) ;
26 accAm = −gamma ∗ omega [ 0 ] ;
27 accEn = f o r c e ∗ s i n ( f r e q ∗ time [ 0 ] ) ;
28 a c c = accG + accAm + accEn ;
29 omega [ 0 ] = omega [ 0 ] − dt / 2 . ∗ a c c ;
30 f o r ( int i =1; i <STEPS ; i ++) {
31 accG = − g /L ∗ s i n ( a n g l e [ i −1]) ;
32 accAm = − gamma ∗ omega [ i − 1 ] ;
33 accEn = f o r c e ∗ s i n ( f r e q ∗ time [ i −1]) ;
34 a c c = accG + accAm + accEn ;
35 omega [ i ] = omega [ i −1] + a c c ∗ dt ;
36 a n g l e [ i ] = a n g l e [ i −1] + omega [ i ] ∗ dt ;
37 time [ i ] = i ∗ dt ;
38 }
39 f o r ( int i =0; i <STEPS ; i ++) // c o n v e r s i o n de l ’ a n g l e en d e g r e e s
40 a n g l e [ i ] ∗= 1 8 0 . / M PI ;
41
42 // i n i t i a l i s a t i o n de d i s l i n
43 m e t a f l ( ”XWIN” ) ; // XWIN ou PDF
44 disini () ;
45 name ( ” temps ( s ) ” , ”X” ) ;
46 name ( ” a n g l e ( d e g r e e s ) ” , ”Y” ) ;
47 graf (0. ,100. ,0. ,10. , −20. ,20. , −20. ,5.) ;
48 // d e s i n de l a t r a j e c t o i r e
49 thkcrv (10) ;
50 c o l o r ( ”RED” ) ;
51 c u r v e ( time , a n g l e , STEPS) ;
52 // f i n de DISLIN
53 disfin () ;
54
55 return 0 ;
56 }

5. Pendule élastique
La décomposition des forces agissant sur la balle attachée au ressort en coordonnèes

6
Figure 4 – Trajectoire du pendule élastique. La courbe verte décrit l’énergie en fonction du
temps décalé de −100 unités, la courbe blue l’angle et la courbe rouge l’élongation décalé de
+100 unités.

polaires est la suivante :

F~ = M~a = M ~r¨
M~g − k(r − L)r̂ = M (r̈ − rϑ̇2 )r̂ + M (rϑ̈ + 2ṙϑ̇)ϑ̂
M (g cos ϑ r̂ − g sin ϑ ϑ̂) − k(r − L)r̂ = M (r̈ − rϑ̇2 )r̂ + M (rϑ̈ + 2ṙϑ̇)ϑ̂

où r est l’élongation du ressort et ϑ l’angle entre le pendule et la verticale. L’on obtient
deux équations du mouvement pour l’élongation du ressort et l’angle, qui sont liées :
k
dans la composante r̂ : r̈ = g cos ϑ − (r − L) + rϑ̇2
M
g ṙϑ̇
dans la composante ϑ̂ : ϑ̈ = − sin ϑ − 2
r r
On doit résoudre les deux équations au simultanément. A chaque itération on calcule
l’accélération angulaire et l’accélération radiale, la vitesse angulaire et la vitesse radiale,
et la position angulaire et radiale avec la méthode de Runge.
La figure 4 montre l’élongation (en rouge), l’angle (en blue) et l’énergie (en vert) du
pendule élastique pour une masse de 0.1 kg attachée au ressort de constante k = 2 N/m
avec une élongation initiale de 2 m et d’angle initial de 15◦ .
Pendule [Link]
1 // p e n d u l e e l a s t i q u e
2 #include <i o s t r e a m >
3 #include <cmath>
4 #include <ctime>
5 #include ” d i s l i n . h” // b i l b i o t h e q u e DISLIN
6
7 using namespace s t d ;
8
9 int main ( ) {
10 const int STEPS = 1 0 0 0 0 ;
11 double dT = 0 . 0 0 5 ;
12 double a n g l e [ STEPS ] , l o n g u e u r [ STEPS ] , temps [ STEPS ] ;
13 double vitA [ STEPS ] , v i t L [ STEPS ] , accA [ STEPS ] , accL [ STEPS ] ;

7
14 double e n e r g i e [ STEPS ] ;
15
16 // p a r a m e t r e s du p e n d u l e
17 double L0 = 1 . ; // l o n g u e u r a r e p o s du r e s s o r t en m e t r e s
18 double k = 2 . ; // c o n s t a n t e e l a s t i q u e du r e s s o r t en N/m
19 double g = 9 . 8 1 ; // a c c e l e r a t i o n g r a v i t a t i n e l l e
20 double m = 0 . 1 ; // masse
21
22 c o u t <<” ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗\ n” ;
23 c o u t << ”Nous avons un p e n d u l e avec un r e s s o r t au l i e u d ’ une c o r d e . \ n”
24 << ”La l o n g u e u r a r e p o s du r e s s o r t e s t ” << L0 << ” m e t \n”
25 << ” l a c o n s t a n t e e l a s t i q u e e s t ” << k << ” N/m. \ n”
26 << ”La masse e s t ” << m << ” kg . \ n”
27 << ”On assume que l ’ o s c i l l a t i o n demarre d e p u i s l ’ e t a t de r e p o s \n”
28 << ” a un c e r t a i n a n g l e e t e l o n g a t i o n i n i t i a l e s . \ n\n” ;
29 c o u t << ” Q u e l l e e s t l ’ a n g l e i n i t i a l du p e n d u l e ? ” ;
30 c i n >> a n g l e [ 0 ] ;
31 c o u t << ”Et l ’ e l o n g a t i o n i n i t i a l e du r e s s o r t en m e t r e s ? ” ;
32 c i n >> l o n g u e u r [ 0 ] ;
33
34 int opt ;
35 do {
36 c o u t << ”Que type de g r a p h i q u e s o u h a i t e z −vous ? ” << e n d l ;
37 c o u t << ” 1 pour d e s s i n e r un g r a p h i q u e ” << e n d l ;
38 c o u t << ” 2 pour une a n i m a t i o n du mouvement du p e n d u l e ” << e n d l ;
39 c i n >> opt ; c o u t << e n d l ;
40 } while ( opt !=1 && opt !=2) ;
41
42 temps [ 0 ] = 0 . ;
43 a n g l e [ 0 ] ∗= M PI / 1 8 0 . ;
44
45 // c e s v e c t e u r s v o n t e t r e u t i l i s e s a v e c DISLIN
46 double x [ STEPS ] , y [ STEPS ] ;
47 x [ 0 ] = sin ( angle [ 0 ] ) ∗ longueur [ 0 ] ;
48 y [ 0 ] = −c o s ( a n g l e [ 0 ] ) ∗ l o n g u e u r [ 0 ] ;
49
50 vitA [ 0 ] = 0 . ;
51 vitL [ 0 ] = 0 . ;
52 accA [ 0 ] = −1.∗ s i n ( a n g l e [ 0 ] ) ∗ g / l o n g u e u r [ 0 ] ;
53 accL [ 0 ] = c o s ( a n g l e [ 0 ] ) ∗ g − k ∗ ( l o n g u e u r [0] − L0 ) /m;
54
55 e n e r g i e [ 0 ] = m∗g ∗y [ 0 ] + k ∗ ( l o n g u e u r [0] − L0 ) ∗ ( l o n g u e u r [0] − L0 ) / 2 . ;
56 c o u t << ”L ’ e n e r g i e i n i t i a l e du systeme e s t : ” << e n e r g i e [ 0 ] << e n d l ;
57
58 // i n i t i a l i s a t i o n de DISLIN
59 m e t a f l ( ”XWIN” ) ; //XWIN ou PDF
60 disini () ;
61
62 // methode de Runge
63 f o r ( int i =0; i <STEPS−1; i ++) {
64 vitA [ i +1] = vitA [ i ] + accA [ i ] ∗ dT ;
65 v i t L [ i +1] = v i t L [ i ] + accL [ i ] ∗ dT ;
66 a n g l e [ i +1]= a n g l e [ i ] + vitA [ i +1]∗dT ;
67 l o n g u e u r [ i +1]= l o n g u e u r [ i ] + v i t L [ i +1]∗dT ;
68 temps [ i +1]=( i +1)∗dT ;
69
70 accA [ i +1] = −1.∗ s i n ( a n g l e [ i +1]) ∗g / l o n g u e u r [ i +1]
71 − 2 . ∗ ( vitA [ i +1]∗ v i t L [ i +1]) / l o n g u e u r [ i + 1 ] ;
72 accL [ i +1] = c o s ( a n g l e [ i +1]) ∗ g − k ∗ ( l o n g u e u r [ i +1]−L0 ) /m
73 + l o n g u e u r [ i +1]∗ vitA [ i +1]∗ vitA [ i + 1 ] ;

8
 74
75 e n e r g i e [ i +1] = m∗ v i t L [ i +1]∗ v i t L [ i + 1 ] / 2 .
76 + m∗ l o n g u e u r [ i +1]∗ l o n g u e u r [ i +1]∗ vitA [ i +1]∗ vitA [ i + 1 ] / 2 .
77 + m∗ g∗y [ i +1] + k ∗ ( l o n g u e u r [ i +1]−L0 ) ∗ ( l o n g u e u r [ i +1]−L0 ) / 2 . ;
78
79 i f ( abs ( l o n g u e u r [ i +1]<=1.0E−6) ) {
80 l o n g u e u r [ i +1]=1.0E−6;
81 c o u t << ” A t t e n t i o n au c o n d i t i o n s i n i t i a l e s e t aux p a r a m e t r e s du
p e n d u l e ! ” << e n d l ;
82 }
83
84 // t a b l e a u x de 1 e l e m e n t pour DISLIN
85 double x p l o t [ 1 ] , y p l o t [ 1 ] ;
86 x p l o t [ 0 ] = s i n ( a n g l e [ i +1]) ∗ l o n g u e u r [ i + 1 ] ;
87 y p l o t [ 0 ] = −c o s ( a n g l e [ i +1]) ∗ l o n g u e u r [ i + 1 ] ;
88 x [ i +1] = x p l o t [ 0 ] ;
89 y [ i +1] = y p l o t [ 0 ] ;
90
91 double l i g n e x [ 2 ] = { 0 . , x [ i + 1 ] } ;
92 double l i g n e y [ 2 ] = { 0 . , y [ i + 1 ] } ;
93 // d e s s i n du p e n d u l e
94 i f ( opt==2) {
95 erase () ; // e f f a c e l a f e n e t r e du g r a p h i q u e
96 c o l o r ( ”WHITE” ) ;
97 name ( ” axe−X” , ”X” ) ;
98 name ( ” axe−Y” , ”Y” ) ;
99 t i t l i n ( ” Animation du p e n d u l e e l a s t i q u e ” , 1 ) ;
100 g r a f ( −4.∗L0 , 4 . ∗ L0 , −4.∗L0 , 1 . ∗ L0 , −5.∗L0 , 1 . ∗ L0 , −5.∗L0 , 1 . ∗ L0 ) ;
101 t i t l e () ;
102
103 incmrk ( −1) ;
104 c o l o r ( ”RED” ) ;
105 marker ( 2 1 ) ;
106 curve ( xplot , yplot , 1 ) ;
107 incmrk ( 0 ) ;
108 thkcrv (1) ;
109 c o l o r ( ”BLUE” ) ;
110 curve (x , y , i ) ;
111 thkcrv (5) ;
112 c o l o r ( ”YELLOW” ) ;
113 curve ( lignex , ligney , 2) ;
114 sendbf () ;
115 endgrf () ;
116 }
117
118 } // f i n de l a b o u c l e
119
120 // d e s s i n de l a t r a j e c t o i r e
121 i f ( opt==1) {
122 name ( ” temps ( s ) ” , ”X” ) ;
123 name ( ” e n e r g i e , a n g l e , e l o n g a t i o n ” , ”Y” ) ;
124 t i t l i n ( ” pendule e l a s t i q u e ” , 1) ;
125 g r a f ( 0 . , STEPS∗dT , 0 . , 1000∗dT , −200.∗L0 , 3 5 0 . ∗ L0 , −200.∗L0 , 1 0 0 . ∗ L0 ) ;
126 t i t l e () ;
127
128 f o r ( int i =0; i < STEPS ; i ++) {
129 a n g l e [ i ] ∗= 180/M PI ;
130 l o n g u e u r [ i ] ∗= 1 0 0 . ;
131 e n e r g i e [ i ] = e n e r g i e [ i ] ∗ 1 0 . −100.;
132 }

9
133
134 thkcrv (10) ;
135 c o l o r ( ”RED” ) ;
136 c u r v e ( temps , l o n g u e u r , STEPS) ;
137 c o l o r ( ”BLUE” ) ;
138 c u r v e ( temps , a n g l e , STEPS) ;
139 c o l o r ( ”GREEN” ) ;
140 c u r v e ( temps , e n e r g i e , STEPS) ;
141 }
142
143 disfin () ; // f i n de DISLIN
144
145 return 0 ;
146 }

6. Orbite elliptique
L’intégration de l’orbite est effectuée avec la méthode de Runge. On redéfinit d’abord la
vitesse initiale à t − ∆t/2 :
∆t
xvel[0] = xvel0 − xacc × ;
2
∆t
yvel[0] = yvel0 − yacc × ;
2
ou ∆t est la pas d’intégration pour la trajectoire. Pour calculer l’accélération à t = 0 on
utilise la distance au centre gravitationnel :
p
dist = xpos[0] × xpos[0] + ypos[0] × ypos[0];
xacc = −GM × xpos[0] / dist3 ;
yacc = −GM × ypos[0] / dist3 ;

GM = GM m. Pour l’algorithme d’intégration, voir le programme Orbite [Link]

développé pendant la leçon. On calcule l’énergie potentielle et cinétique à chaque pas
avec les équations Epot = −GM m/r et Ecin = 21 mv 2 . La figure 5 montre l’évolution de
l’énergie (en unités AU, M , années). La somme de l’énergie potentielle et cinétique reste
constante.
[Link]
1 // e t u d e de l ’ o r b i t e a v e c l a methode de Runge
2 #include <i o s t r e a m >
3 #include <c s t d l i b >
4 #include <cmath>
5 #include ” d i s l i n . h”
6
7 using namespace s t d ;
8
9 int main ( ) {
10 const int s t e p s = 2 0 0 0 0 ;
11 const double GM = 1 . ; //G N ∗ M S o l e i l
12 const double m = 1 . ; // masse p l a n e t e
13 double xpos [ s t e p s ] , ypos [ s t e p s ] , x v e l [ s t e p s ] , y v e l [ s t e p s ] , time [ s t e p s ] ;
14 double ePot [ s t e p s ] , eCin [ s t e p s ] , eTot [ s t e p s ] , l [ s t e p s ] ;
15 // c o n d i t i o n s i n i t i a l e s
16 xpos [ 0 ] = 1 . 0 ; ypos [ 0 ] = 0 . 0 ;
17 xvel [ 0 ] = 0 . 0 ; yvel [ 0 ] = 1 . 3 ;
18 time [ 0 ] = 0 . ;
19

10
Figure 5 – Evolution de l’énergie cinétique (bleu) et potentielle (rouge) de l’orbite. Le mo-
ment cinétique (jaune) et l’énergie totale (vert) sont constants. A noter que l’énergie totale est
négative, comme elle doit être pour une orbite fermé.

20 double dT = 0 . 0 1 ; // pas d ’ i n t e g r a t i o n
21
22 // methode de Runge
23 double d i s t , v e l , xacc , yacc ; // v a r i a l b e s a u x i l i a r e s
24 // r e d e f i n i t i o n de l a v i t e s s e i n i t i a l e
25 d i s t = s q r t ( xpos [ 0 ] ∗ xpos [ 0 ] + ypos [ 0 ] ∗ ypos [ 0 ] ) ;
26 vel = sqrt ( xvel [ 0 ] ∗ xvel [ 0 ] + yvel [ 0 ] ∗ yvel [ 0 ] ) ;
27 xacc = −GM ∗ xpos [ 0 ] / ( d i s t ∗ d i s t ∗ d i s t ) ; // f o r c e ou a c c e l e r a t i o n , GM = 1
28 yacc = −GM ∗ ypos [ 0 ] / ( d i s t ∗ d i s t ∗ d i s t ) ;
29 x v e l [ 0 ] = x v e l [ 0 ] − xacc ∗dT / 2 . ; // v i t e s s e
30 y v e l [ 0 ] = y v e l [ 0 ] − yacc ∗dT / 2 . ;
31 f o r ( int i =1; i <s t e p s ; i ++) {
32 d i s t = s q r t ( xpos [ i −1]∗ xpos [ i −1] + ypos [ i −1]∗ ypos [ i −1]) ;
33 v e l = s q r t ( x v e l [ i −1]∗ x v e l [ i −1] + y v e l [ i −1]∗ y v e l [ i −1]) ;
34 xacc = −GM ∗ xpos [ i −1] / ( d i s t ∗ d i s t ∗ d i s t ) ;
35 yacc = −GM ∗ ypos [ i −1] / ( d i s t ∗ d i s t ∗ d i s t ) ;
36 x v e l [ i ] = x v e l [ i −1] + xacc ∗dT ; // v i t e s s e
37 y v e l [ i ] = y v e l [ i −1] + yacc ∗dT ;
38 xpos [ i ] = xpos [ i −1] + x v e l [ i ] ∗ dT ; // p o s i t i o n
39 ypos [ i ] = ypos [ i −1] + y v e l [ i ] ∗ dT ;
40 time [ i ] = i ∗dT ;
41 ePot [ i ] = −GM ∗ m / d i s t ;
42 eCin [ i ] = 0 . 5 ∗ m ∗ v e l ∗ v e l ;
43 eTot [ i ] = ePot [ i ] + eCin [ i ] ;

11
44 //L = m ∗ ( r x v ) ! ! ! r e t v ne s o n t pas t o u j o u r s o r t h o g o n a u x
45 l [ i ] = m ∗ ( xpos [ i ] ∗ y v e l [ i ] − ypos [ i ] ∗ x v e l [ i ] ) ;
46 }
47
48 // p a r t i e g r a p h i q u e a v e c DISLIN
49 m e t a f l ( ”XWIN” ) ; //XWIN ou PDF
50 page ( 3 0 0 0 , 3 0 0 0 ) ; // pour c r e e r une image p r o p o r t i o n e e
51 disini () ; // i n i t i a l i s a t i o n de DISLIN
52
53 // d e s s i n de l ’ o r b i t e
54 name ( ” axe−X [UA] ” , ”X” ) ;
55 name ( ” axe−Y [UA] ” , ”Y” ) ;
56 t i t l i n ( ” Orbite e l l i p t i q u e ” ,1) ;
57 graf ( −7. ,3. , −7. ,1. , −5. ,5. , −5. ,1.) ;
58 t i t l e () ;
59 thkcrv (5) ;
60 c o l o r ( ”RED” ) ;
61 c u r v e ( xpos , ypos , s t e p s ) ;
62 // d e s s i n du s o l e i l
63 incmrk ( −1) ;
64 c o l o r ( ”GREEN” ) ;
65 marker ( 2 1 ) ;
66 double xx , yy ;
67 xx = 0 . ;
68 yy = 0 . ;
69 c u r v e (&xx ,&yy , 1 ) ; // l a f o n c t i o n s ’ a t t e n d un a d r e s s e memoire
70 endgrf () ;
71 erase () ;
72 system ( ”PAUSE” ) ; // a v a n t de d e s s i n e r l e deuxieme graph
73
74 // d e s s i n de l ’ e n e r g i e
75 name ( ” temps ” , ”X” ) ;
76 name ( ” e n e r g i e ” , ”Y” ) ;
77 t i t l i n ( ” Orbite e l l i p t i q u e ” ,1) ;
78 t i t l i n ( ” j a u n e = moment c i n e t i q u e ” , 2 ) ;
79 t i t l i n (” vert = energie total ” ,3) ;
80 t i t l i n ( ” rouge / blue = e n e r g i e p o t e n t i e l e / c i n e t i q u e ” ,4) ;
81 c o l o r ( ”FORE” ) ;
82 g r a f ( 0 . , dT∗ s t e p s , 0 . , dT∗ s t e p s / 5 . , −2. , 2 . , −2. , 0 . 5 ) ;
83 t i t l e () ;
84 incmrk ( 0 ) ;
85 thkcrv (2) ;
86 c o l o r ( ”RED” ) ;
87 c u r v e ( time , ePot , s t e p s ) ;
88 c o l o r ( ”BLUE” ) ;
89 c u r v e ( time , eCin , s t e p s ) ;
90 c o l o r ( ”GREEN” ) ;
91 c u r v e ( time , eTot , s t e p s ) ;
92 c o l o r ( ”YELLOW” ) ;
93 c u r v e ( time , l , s t e p s ) ;
94
95 disfin () ; // f i n de DISLIN
96
97 return 0 ;
98 }

7. Orbite avec frottement

Un satellite est sur une orbite circulaire proche de la Terre, à 100 km de la surface. La force
de frottement due à l’atmosphère représente 0.01% de la force gravitationnelle (situation

12
 hypothétique). La force de frottement est toujours parallèle à la vitesse instantanée. Es-
timez le temps qu’il faut pour que le satellite descende à 50 km.
On redéfinit la vitesse a t − ∆T /2 pour utiliser la méthode d’intégration de Runge :

∆t
xvel[0] = xvel0 − xacc × ;
2
∆t
yvel[0] = yvel0 − yacc × ;
2
Pour calculer l’accélération à t = 0 on calcule la distance au centre gravitationnel et on
utilise les formules développées pendant le cours (voir [Link]) pour calculer les
forces avec le frottement :
GM m  x y
Fx = − 3 − CX
r r r
GM m  y x
Fy = − 3 + CX
r r r
avec les instructions :
p
dist = xpos[0] × xpos[0] + ypos[0] × ypos[0];
xacc = −GM/dist3 × (xpos[i] − CX × ypos[i]);
yacc = −GM/dist3 × (ypos[i] + CX × xpos[i]);

Ensuite on étudie l’orbite avec la méthode de Runge : la position à t = (i + 1) × ∆T

est calculée en utilisant la vitesse au pas i + 1, qui correspond à la vitesse au temps
t = (i + 1) × ∆T /2. A chaque étape la distance du satellite à la Terre est calculée et
comparée à la hauteur finale limite (altfin). Quand la distance à la Terre est plus petite
que altfin le programme sort de la boucle et affiche le nombre de tours autours de la
Terre.

Figure 6 – L’orbite du satellite avec frottement, à gauche le mouvement complet et à droite

un zoom qui montre la diminution du rayon de l’orbite à chaque rotation.

[Link]
1 // e t u d e d ’ une o r b i t e a v e c f r o t t e m e n t dans l a h a u t atmosphere
2 #include <i o s t r e a m >
3 #include <cmath>

13
 4 #include ” d i s l i n . h”
5
6 using namespace s t d ;
7
8 int main ( ) {
9 const double RT = 6 3 7 8 . ; // en km
10 const double GM = 3 9 8 6 4 5 . ; //G N ∗ m a s s e t e r r e
11 double a l t I n i = 1 0 0 . ;
12 double a l t F i n = 5 0 . ;
13 const double R0 = RT + a l t I n i ;
14 const double V0 = s q r t (GM/R0) ;
15 const double CX = 0 . 0 0 0 1 ; // c o e f f i c i e n t de f r o t t e m e n t
16
17 const int s t e p s = 5 0 0 0 0 ;
18 double xpos [ s t e p s ] , ypos [ s t e p s ] , x v e l [ s t e p s ] , y v e l [ s t e p s ] ;
19 // c o n d i t i o n s i n i t i a l e s
20 xpos [ 0 ] = R0 ; ypos [ 0 ] = 0 . ;
21 x v e l [ 0 ] = 0 . ; y v e l [ 0 ] = V0 ;
22
23 double dT = 2 . ; // pas d ’ i n t e g r a t i o n
24
25 // methode de Runge
26 double d i s t , xacc , yacc ; // v a r i a b l e s a u x i l i a r e s
27 // r e d i f i n i t i o n de l a v i t e s s e i n i t i a l e
28 d i s t = s q r t ( xpos [ 0 ] ∗ xpos [ 0 ] + ypos [ 0 ] ∗ ypos [ 0 ] ) ;
29 xacc = −GM / pow ( d i s t , 3 ) ∗ ( xpos [ 0 ] − CX∗ ypos [ 0 ] ) ; // a c c e l e r a t i o n
30 yacc = −GM / pow ( d i s t , 3 ) ∗ ( ypos [ 0 ] + CX∗ xpos [ 0 ] ) ;
31 x v e l [ 0 ] = x v e l [ 0 ] −xacc ∗dT / 2 . ; // v i t e s s e
32 y v e l [ 0 ] = y v e l [ 0 ] −yacc ∗dT / 2 . ;
33 int n s t e p s = 0 ;
34 int n o r b i t s = 0 ;
35 // b o u c l e
36 f o r ( int i =1; i <s t e p s ; i ++) {
37 d i s t = s q r t ( xpos [ i −1]∗ xpos [ i −1] + ypos [ i −1]∗ ypos [ i −1]) ;
38 xacc = −GM / pow ( d i s t , 3 ) ∗ ( xpos [ i −1] − CX∗ ypos [ i −1]) ;
39 yacc = −GM / pow ( d i s t , 3 ) ∗ ( ypos [ i −1] + CX∗ xpos [ i −1]) ;
40
41 x v e l [ i ] = x v e l [ i −1] + xacc ∗dT ; // v i t e s s e
42 y v e l [ i ] = y v e l [ i −1] + yacc ∗dT ;
43
44 xpos [ i ] = xpos [ i −1] + x v e l [ i ] ∗ dT ; // p o s i t i o n
45 ypos [ i ] = ypos [ i −1] + y v e l [ i ] ∗ dT ;
46
47 n s t e p s ++;
48 i f ( ypos [ i −1]∗ ypos [ i ] < 0 . ) n o r b i t s = n o r b i t s + 1 ;
49 i f ( d i s t < RT) { // l e s a t e l l i t e e s t tombe s u r l a Terre
50 c o u t << ”Le s a t e l l i t e e s t tombe s u r l a T e r r e ! ” << e n d l ;
51 break ;
52 }
53 }
54
55 c o u t << ” a p r e s ” << int ( n s t e p s ∗dT) << ” s e c o n d s ” << e n d l ;
56 c o u t << ” e t ” << n o r b i t s /2 << ” o r b i t e s ” ;
57 i f ( n o r b i t s%2 == 1 ) c o u t << ” e t demi ” ;
58 c o u t << e n d l ;
59
60 // p a r t i e g r a p h i q u e v e c DISLIN
61 m e t a f l ( ”XWIN” ) ; //XWIN ou PDF
62 page ( 3 0 0 0 , 3 0 0 0 ) ;
63 disini () ; // i n i t i a l i s a t i o n de DISLIN

14
64 name ( ” axe−X” , ”X” ) ;
65 name ( ” axe−Y” , ”Y” ) ;
66 graf ( −7500. ,7500. , −7500. ,2500. , −7500. ,7500. , −7500. ,2500.) ;
67 thkcrv (2) ;
68 color (” yellow ”) ;
69 c u r v e ( xpos , ypos , n s t e p s ) ;
70
71 disfin () ; // f i n de DISLIN
72
73 return 0 ;
74 }

8. E cross B
Une particule de masse M et charge Q est soumise un champ électrique E = (0, Ey , 0)
et magnétique B = (0, 0, Bz ). La particule est au repos au point P = (0, 0, 0) quand le
champ électrique est allumé. Etudiez la trajectoire de la particule et dessinez la. Utilisez
la méthode de Runge pour intégrer l’équation du mouvement. Quelle est la vélocité de
dérive de la particule ? Et dans quelle direction ? Estimez la à partir de la trajectoire
(paramètres : M = 0.001 kg, Q = 0.001 C, E = 100 V, B = 1 T).
D’abord on analyse toutes les forces qui agissent sur la particule :

F = ma = Q(E + v × B)
Fx = max = Qvy B
Fy = may = −Qvx B + QE
Fz = 0

Vu que vz,0 = 0 et Fz = 0 la particule se déplace dans le plan x − y.

On redéfinit les conditions initiales pour intgrer l’équation du mouvement avec la méthode
de Runge :
Q
ax,0 = vy,0 B
m
Q Q
ay,0 = − vx,0 B+ E
m m
et
∆t
vx,0 = vx0 − ax,0
2
∆t
vy,0 = vy0 − ay,0
2
Puis on utilise le même algorithme que l’on a développé pour étudier les orbites, mais on
remplace la force gravitationnelle par la force définie précedemment. La trajectoire de la
particule est montrée dans la figure 7.
Pour estimer la vitesse de dérive on mesure le temps entre deux changements successives
de signe de la composante y de la vitesse.
[Link]
1 // e t u d e de l a t r a j e c t o r i e d ’ une p a r t i c u l e dans
2 //un champ e l e c t r o m a g n e t i q u e
3 #include <i o s t r e a m >
4 #include <cmath>
5 #include ” d i s l i n . h”
6

15
 Figure 7 – Trajectoire d’une particule charge dans un champ électromagnétique.

7 using namespace s t d ;
8
9 int main ( ) {
10 // p a r a m e t r e s
11 const int s t e p s = 6 0 0 0 0 ;
12 const double M = 0 . 0 0 1 ; // masse
13 const double Q = 0 . 0 0 1 ; // c h a r g e
14 const double E = 1 0 0 . ; //E
15 const double B = 1 . ; //B
16 double xpos [ s t e p s ] , ypos [ s t e p s ] , x v e l [ s t e p s ] , y v e l [ s t e p s ] ;
17 // c o n d i t i o n s i n i t i a l e s
18 xpos [ 0 ] = 0 . 0 ; ypos [ 0 ] = 0 . 0 ;
19 xvel [ 0 ] = 0 . 0 ; yvel [ 0 ] = 0 . 0 ;
20
21 double dT = 0 . 0 0 0 5 ; // pas d ’ i n t e g r a t i o n
22
23 // methode de Runge
24 double xacc , yacc ; // v a r i a l b e s a u x i l i a r e s
25 double v d r i f t ;
26 // r e d e f i n i t i o n de c o n d i t i o n s i n i t i a l e s
27 xacc = Q/M∗ y v e l [ 0 ] ∗ B ; // a c c e l e r a t i o n
28 yacc = −Q/M∗ x v e l [ 0 ] ∗ B + Q/M∗E ;
29 x v e l [ 0 ] = x v e l [ 0 ] − xacc ∗dT / 2 . ; // v i t e s s e
30 y v e l [ 0 ] = y v e l [ 0 ] − yacc ∗dT / 2 . ;
31
32 f o r ( int i =1; i <s t e p s ; i ++) {
33 xacc = Q/M∗ y v e l [ i −1]∗B ; // a c c e l e r a t i o n
34 yacc = −Q/M∗ x v e l [ i −1]∗B + Q/M∗E ;
35 x v e l [ i ] = x v e l [ i −1] + xacc ∗dT ; // v i t e s s e
36 y v e l [ i ] = y v e l [ i −1] + yacc ∗dT ;
37 xpos [ i ] = xpos [ i −1] + x v e l [ i ] ∗ dT ; // p o s i t i o n
38 ypos [ i ] = ypos [ i −1] + y v e l [ i ] ∗ dT ;
39
40 // v i t e s s e de d e r i v e

16
41 i f ( y v e l [ i −1] >0. && y v e l [ i ] < 0 . ) {
42 v d r i f t = xpos [ i ] / ( i ∗dT) ;
43 c o u t << ”La v i t e s s e de d e r i v e e s t : ” << v d r i f t << ” m/ s . ” << e n d l ;
44 }
45 }
46
47 // p a r t i e g r a p h i q u e a v e c DISLIN
48 m e t a f l ( ”XWIN” ) ; //XWIN ou PDF
49 disini () ; // i n i t i a l i s a t i o n de DISLIN
50
51 // d e s s i n de l a t r a j e c t o i r e
52 name ( ” axe−X [m] ” , ”X” ) ;
53 name ( ” axe−Y [m] ” , ”Y” ) ;
54 t i t l i n ( ”Champ E l e c t r o M a g n e t i q u e ” , 1 ) ;
55 graf (0. ,3100. ,0. ,500. , −20. ,220. , −20. ,40.) ;
56 t i t l e () ;
57 thkcrv (5) ;
58 c o l o r ( ”RED” ) ;
59 c u r v e ( xpos , ypos , s t e p s ) ;
60 disfin () ; // f i n de DISLIN
61
62 return 0 ;
63 }

9. Système binaire
Etudiez un système binaire pour deux objets de même masse. Il vous faut écrire 4
équations différentielles de premier ordre pour chaque objet.
Les étoiles ont une masse mA et mB respectivement, et la force qui agit sur chacune est
donnée par les équations :
GN mB mA
mA~aA = − (~rB − ~rA ) (5)
|~rB − ~rA |3
GN mA mB
mB~aB = − (~rA − ~rB )
|~rB − ~rA |3

Les 4 équations différentielles pour l’étoile A sont :

GN mB
v̇x,A = − (xB − xA ) (6)
|~rB − ~rA |3
GN mB
v̇y,A = − (yB − yA )
|~rB − ~rA |3
ẋA = vx,A
ẏA = vy,A

Par la conservation du moment cinétique,p la composante z peut être negligée. La distance

entre les deux étoiles est |~rB − ~rA | = (xA − xB )2 + (yA − yB )2 . Pour l’étoile B les
équations sont similaires, mais en remplaçant A par B et vice versa.
On a aussi ajouté une animation du système binaire. Pour ralentir l’exécution du pro-
gramme on utilise la fonction sleep, qui malheureusement est différente sur Windows
et Unix. Pour pouvoir utiliser le même programme sur Windows et Unix, on utilise les
directives de précompilation pour choisir la bonne fonction.
[Link]
1 // e t u d e d e s o r b i t e s de deux e t o i l e s b i n a i r e s
2 // methode de Runge

17
 3 #include <i o s t r e a m >
4 #include <cmath>
5 #include ” d i s l i n . h”
6
7 /∗
8 f o n c t i o n SLEEP(ms)
9 malheureusement l a f o n c t i o n SLEEP e s t d i f f e r e n t e s u r Windows e t Unix : −(
10 pour p o u v o i r u t i l i s e r l e meme programme s u r Windows e t Unix ,
11 on u t i l i s e l e d i r e c t i v e s de p r e c o m p i l a t i o n pour c h o i s i r l a bonne f o n c t i o n
12 ∗/
13 #i f d e f UNIX
14 #include <u n i s t d . h>
15 #define SLEEP(ms) u s l e e p (ms∗ 1 0 0 0 ) ;
16 #e l s e
17 #include <windows . h>
18 #define SLEEP(ms) S l e e p (ms) ;
19 #endif
20
21 using namespace s t d ;
22
23 int main ( ) {
24 const int s t e p s = 3 0 0 0 0 ;
25 const double GN = 1 . ;
26 // a t t e n t i o n aux c o n d i t i o n s i n i t i a l e s
27 // e t o i l e A
28 const double MA = 1 . ;
29 double xposA [ s t e p s ] , yposA [ s t e p s ] , xvelA [ s t e p s ] , yvelA [ s t e p s ] ;
30 xposA [ 0 ] = 0 . 3 ; yposA [ 0 ] = 0 . ;
31 xvelA [ 0 ] = 0 . ; yvelA [ 0 ] = 1 . 0 ;
32 // e t o i l e B
33 const double MB = MA;
34 double xposB [ s t e p s ] , yposB [ s t e p s ] , xvelB [ s t e p s ] , yvelB [ s t e p s ] ;
35 xposB [ 0 ] = −xposA [ 0 ] ; yposB [ 0 ] = 0 . ;
36 xvelB [ 0 ] = 0 . ; yvelB [ 0 ] = −yvelA [ 0 ] ;
37
38 double dT = 0 . 0 0 5 ;
39
40 // methode de Runge
41 double d i s t , xaccA , yaccA , xaccB , yaccB ; // v a r i a b l e s a u x i l i a r e s
42 // r e d e f i n i t i o n d e s v i t e s s e s i n i t i a l e s
43 d i s t = s q r t ( pow ( ( xposA [0] − xposB [ 0 ] ) , 2 ) + pow ( ( yposA [0] − yposB [ 0 ] ) , 2 ) ) ;
44 // e t o i l e A
45 xaccA = −GN∗MB / pow ( d i s t , 3 ) ∗ ( xposA [0] − xposB [ 0 ] ) ;
46 yaccA = −GN∗MB / pow ( d i s t , 3 ) ∗ ( yposA [0] − yposB [ 0 ] ) ;
47 xvelA [ 0 ] = xvelA [ 0 ] −xaccA ∗dT / 2 . ; // v i t e s s e
48 yvelA [ 0 ] = yvelA [ 0 ] −yaccA ∗dT / 2 . ;
49 // e t o i l e B
50 xaccB = −GN∗MA / pow ( d i s t , 3 ) ∗ ( xposB [0] − xposA [ 0 ] ) ;
51 yaccB = −GN∗MA / pow ( d i s t , 3 ) ∗ ( yposB [0] − yposA [ 0 ] ) ;
52 xvelB [ 0 ] = xvelB [ 0 ] −xaccB ∗dT / 2 . ; // v i t e s s e
53 yvelB [ 0 ] = yvelB [ 0 ] −yaccB ∗dT / 2 . ;
54 f o r ( int i =1; i <s t e p s ; i ++) {
55 // d i s t a n c e e n t r e l e s e t o i l e s
56 dist =
57 s q r t ( pow ( ( xposA [ i −1]−xposB [ i −1]) , 2 )+pow ( ( yposA [ i −1]−yposB [ i −1]) , 2 ) ) ;
58
59 // e t o i l e A
60 xaccA = −GN∗MB / pow ( d i s t , 3 ) ∗ ( xposA [ i −1]−xposB [ i −1]) ;
61 yaccA = −GN∗MB / pow ( d i s t , 3 ) ∗ ( yposA [ i −1]−yposB [ i −1]) ;
62

18
 63 xvelA [ i ] = xvelA [ i −1] + xaccA ∗dT ; // v i t e s s e
64 yvelA [ i ] = yvelA [ i −1] + yaccA ∗dT ;
65
66 xposA [ i ] = xposA [ i −1] + xvelA [ i ] ∗ dT ; // p o s i t i o n
67 yposA [ i ] = yposA [ i −1] + yvelA [ i ] ∗ dT ;
68
69 // e t o i l e B
70 //A e t B s o n t i n v e r s e dans l e c a l c u l de l ’ a c c e l e r a t i o n ! ! !
71 xaccB = −GN∗MA / pow ( d i s t , 3 ) ∗ ( xposB [ i −1]−xposA [ i −1]) ;
72 yaccB = −GN∗MA / pow ( d i s t , 3 ) ∗ ( yposB [ i −1]−yposA [ i −1]) ;
73
74 xvelB [ i ] = xvelB [ i −1] + xaccB ∗dT ; // v i t e s s e
75 yvelB [ i ] = yvelB [ i −1] + yaccB ∗dT ;
76
77 xposB [ i ] = xposB [ i −1] + xvelB [ i ] ∗ dT ; // p o s i t i o n
78 yposB [ i ] = yposB [ i −1] + yvelB [ i ] ∗ dT ;
79 }
80
81 // g r a p h i q u e DISLIN
82 m e t a f l ( ”XWIN” ) ;
83 page ( 3 0 0 0 , 3 0 0 0 ) ;
84 disini () ;
85 name ( ” axe−X” , ”X” ) ;
86 name ( ” axe−Y” , ”Y” ) ;
87 t i t l i n ( ” Animation d ’ un systeme b i n a i r e ” , 1 ) ;
88
89 /∗
90 // d e s s i n d e s o r b i t e s
91 graf ( −1. ,1. , −1. ,1. , −1. ,1. , −1. ,1.) ;
92 thkcrv (2) ;
93 c o l o r (” r e d ”) ;
94 c u r v e ( xposA , yposB , s t e p s ) ;
95 c o l o r (” y e l l o w ”) ;
96 c u r v e ( xposB , yposB , s t e p s ) ;
97 ∗/
98
99 // animation DISLIN
100 double xx [ 1 ] , yy [ 1 ] ; // v a r i a b l e s a u x i l i e r e s
101 f o r ( int j =1 ; j <s t e p s ; j ++) {
102 // pour a c c e l e r e r l ’ animation on ne d e s s i n e que t o u t e s l e s 100
iterations
103 SLEEP( 5 ) ; // pour a r r e t e r l ’ e x e c u t i o n pendant 10 ms
104 i f ( j %10==0) {
105 erase () ; // on e f f a c e l e t a b l e a u
106 c o l o r ( ”FORE” ) ; // on remet un axe d e c r i t en d e s s o u s en b l a n c
107 graf ( −0.5 ,0.5 , −0.5 ,0.2 , −0.5 ,0.5 , −0.5 ,0.2) ;
108 t i t l e () ;
109
110 // d e s s i n e r l e s e t o i l e s
111 incmrk ( −1) ;
112 c o l o r ( ”GREEN” ) ;
113 marker ( 2 1 ) ;
114 xx [ 0 ] = xposA [ j ] ;
115 yy [ 0 ] = yposA [ j ] ;
116 c u r v e ( xx , yy , 1 ) ;
117
118 incmrk ( −1) ;
119 c o l o r ( ”RED” ) ;
120 marker ( 2 1 ) ;
121 xx [ 0 ] = xposB [ j ] ;

19
122 yy [ 0 ] = yposB [ j ] ;
123 c u r v e ( xx , yy , 1 ) ;
124
125 sendbf () ;
126 endgrf () ;
127 }
128 }
129
130 disfin () ; // f i n de DISLIN
131
132 return 0 ;
133 }

10. Classe Planete

Voir le programme [Link] développé pendant le cours. Le programme est aussi
décomposé en plusieurs fichiers (voir le projet Planete) : Planete.h qui contient la
définition de la classe Planete, [Link] qui contient la définition des fonctions
membres de la classe Planete, et [Link] le programme principal qui utilise la classe
Planete. Chaque planète est définie par son nom, son rayon orbital, sa période orbitale,
l’excentricité de l’orbite et enfin l’inclinaison de l’orbite. La vitesse initiale et la position
initiale de chaque planète sont calculées à partir de ces données.

20