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
75
76
77
78
79
80
81
82
83
84
85
86
87
88
| # -*- coding: utf-8 -*-
from numpy import*
from scipy import*
import matplotlib.pylab as plt
#Condition initiales
l1,l2,m1,m2=1,1,1,1 #Longueur et masses des pendules
t0,tf= 0,2 #Temps de l'experience
th1_0,th2_0,om1_0,om2_0=0,3,0,0 #Angles et vitesses initiales
g=10
#Les deux equations du pendule double
dom1=lambda t,om1,om2,th1,th2 : (-m2*l1*sin(th1-th2)*cos(th1-th2)*om1**2-m2*l2*sin(th1-th2)*om2**2-m1*g*sin(th1-th2)*cos(th2))/l1*(m1+m2*(sin(th1-th2))**2)
dom2=lambda t,om1,om2,th1,th2 : ((m1+m2)*l1*sin(th1-th2)*om1**2+m2*l2*sin(th1-th2)*cos(th1-th2)*om2**2+(m1+m2)*g*sin(th1-th2)*cos(th1))/l2*(m1+m2*(sin(th1-th2))**2)
hmin,h,hmax=0.01,0.1,0.5 #Pas initial
t,th1,th2,om1,om2 = t0,th1_0,th2_0,om1_0,om2_0 #Initialisation des variables
th1_1,th2_1,om1_1,om2_1=th1_0,th2_0,om1_0,om2_0
p=0.001 #Precision demandée
#Tableaux des données
T=array([t])
Th1=array([th1])
Th2=array([th2])
Om1=array([om1])
Om2=array([om2])
#Boucle Runge-Kutta 4
while t<tf:
if t+h>tf: h=tf-t
th1 += h*om1
th2 += h*om2
k1 = h*dom1 (t,om1,om2,th1,th2)
j1 = h*dom2 (t,om1,om2,th1,th2)
k2 = h*dom1 (t+h/2,om1+h*k1/2,om2+h*j1/2,th1,th2)
j2 = h*dom2 (t+h/2,om1+h*k1/2,om2+h*j1/2,th1,th2)
k3 = h*dom1 (t+h/2,om1+h*k2/2,om2+h*j2/2,th1,th2)
j3 = h*dom2 (t+h/2,om1+h*k2/2,om2+h*j2/2,th1,th2)
k4 = h*dom1 (t+h,om1+h*k3,om2+h*j2,th1,th2)
j4 = h*dom2 (t+h,om1+h*k3,om2+h*j2,th1,th2)
om1 += k1/6 + k2/3 + k3/3 + k4/6
om2 += j1/6 + j2/3 + j3/3 + j4/6
h=h/2 #Division du pas par 2
for k in range (2): #Double boucle Runge Kutta avec un pas divisé par 2
th1_1 += h*om1_1 #===> Meilleur precision
th2_1 += h*om2_1
k1_1 = h*dom1 (t,om1_1,om2_1,th1_1,th2_1)
j1_1 = h*dom2 (t,om1_1,om2_1,th1_1,th2_1)
k2_1 = h*dom1 (t+h/2,om1_1+h*k1_1/2,om2_1+h*j1_1/2,th1_1,th2_1)
j2_1 = h*dom2 (t+h/2,om1_1+h*k1_1/2,om2_1+h*j1_1/2,th1_1,th2_1)
k3_1 = h*dom1 (t+h/2,om1_1+h*k2_1/2,om2_1+h*j2_1/2,th1_1,th2_1)
j3_1 = h*dom2 (t+h/2,om1_1+h*k2_1/2,om2_1+h*j2_1/2,th1_1,th2_1)
k4_1 = h*dom1 (t+h,om1_1+h*k3_1,om2_1+h*j2_1,th1_1,th2_1)
j4_1 = h*dom2 (t+h,om1_1+h*k3_1,om2_1+h*j2_1,th1_1,th2_1)
om1_1 += k1_1/6 + k2_1/3 + k3_1/3 + k4_1/6
om2_1 += j1_1/6 + j2_1/3 + j3_1/3 + j4_1/6
h=h*2 #On retablit le pas
err=array([abs(om1-om1_1),abs(om2-om2_1)]) #Calcule de l'erreur
err=max(err)
if err<=p: #Si l'erreur est acceptable
T=append(T,t) #Stockage des données
Th1=append(Th1,th1)
Th2=append(Th2,th2)
Om1=append(Om1,om1)
Om2=append(Om2,om2)
#Augmentation du pas
t+=h #et on passe a l'etape suivante
else: #Sinon diminution du pas
#Et on recommence la meme etape
h = h *min( max( 0.9 * (p/err) **0.2,0.1),4.0)
if h > hmax:
h = hmax
if h < hmin:
h=hmin |
