Programme pendule.py

#!/usr/bin/env pythonw2.3

# lignes utilisees pour matplotlib
import matplotlib
matplotlib.use("WXAgg")
import matplotlib.pylab as pylab
import numpy as N
import math


#############################################################################
# debut du programme
#############################################################################

def pendule_lineaire (tt, vt, at, gR, dt) :
	"""calcule t (t+1), v(t+1), a(t+1)
	t'' + gR t = 0"""
	at = - gR * tt
	vt += at * dt
	tt += vt * dt 
	return tt,vt,at

def pendule_non_lineaire (tt, vt, at, gR, dt) :
	"""calcule t (t+1), v(t+1), a(t+1)
	t'' + gR t = 0"""
	at = - gR * math.sin (tt)
	vt += at * dt
	tt += vt * dt 
	return tt,vt,at

def periode (courbe,dt) :
	"""estime la période d'une courbe dont on sait qu'elle est cyclique"""
	lp = []
	p = 0
	for i in range (1, len (courbe)-1) :
		if courbe [i-1] <= courbe [i] and courbe [i] >= courbe [i+1] :
			# maximum
			lp.append (i)
	# moyenne entre les durees entre deux maximas
	s = 0
	for l in range (1,len (lp)) :
		s += lp [l] - lp [l-1]
	return float (s) / (len (lp)-1) * dt
	
        
# initialisation : un angle t0 != 0, pas de vitesse, pas d'accélération
t0			= 1
tt,vt,at 	= t0,0,0
stt,svt,sat = t0,0,0
dt 			= 0.01
gR 			= 0.5
t 			= 0

T 			= []
th 			= []
sth 		= []
cth 		= []

# boucle
print "simultation"
for i in range (0,4000) :
        T.append (t)
        th.append (tt)
        sth.append (stt)
        ctt 		= t0 * math.cos (math.sqrt (gR) * t)
        cth.append (ctt)
        t 			= t + dt
        tt,vt,at 	= pendule_lineaire (tt, vt, at, gR, dt)
        stt,svt,sat = pendule_non_lineaire (stt, svt, sat, gR, dt)
        
print "periode theorique : ", math.pi * 2 / math.sqrt (gR)
print "periode courbe theorique ", periode (cth, dt)
print "periode courbe lineaire ", periode (th, dt)
print "periode courbe non lineaire ", periode (sth, dt)

# pour calculer l'angle a partir duquel la solution lineaire
# n'est plus viable, il faut calculer la difference entre les periodes
# et choisir un seuil au dela duquel la difference est trop grande
        
# affichage des premieres valeurs
print "premieres valeurs "
print "angle initial : ", t0 * 180 / math.pi
for i in range (0,10) :
        print T[i], " : ", th[i], "\t", cth [i], "\t", sth [i]


# matplotlib n'accepte que des donnees stockees dans un type de tableau
# definies dans le package Numeric
aT = N.array (T)
ath = N.array (th)
acth = N.array (cth)
asth = N.array (sth)

# dessin de la courbe
print "dessin de la courbe"
pylab.title("Oscillation d'un pendule")
pylab.xlabel("t")
pylab.ylabel("theta")
pylab.grid (True)
pylab.plot (aT, ath, linewidth=1.0, color = "blue")
pylab.plot (aT, acth, linewidth=1.0, color = "red")
pylab.plot (aT, asth, linewidth=1.0, color = "green")
pylab.show ()