Corrigé TD2 ¶
Création d'un histogramme¶
In [1]:
from pylab import *
%matplotlib inline
In [2]:
# une première version sans filtrage
def Histogram(Liste,(xmin, xmax),N=10):
Deltax = (xmax-xmin)/float(N)
H = [0]*N
for x in Liste:
ind = min(N-1,int((x-xmin)//Deltax))
H[ind] += 1
# abcisses du début des boîtes
bins = [xmin + n*Deltax for n in range(N+1)]
# abscisses pour l'histogramme avec dédoublement
X = sorted(bins + bins[1:-1])
# ordonnées pour l'histogramme avec normalisation
norm = Deltax*float(sum(H))
for i in range(len(H)): H[i] /= norm
# dédoublement des indices pour les ordonnées
Y = [ H[i] for i in sorted(range(len(H))*2) ]
plot(X,Y,'b')
return (H,bins)
In [3]:
# génération de liste aléatoires
from random import uniform
M = 100000 # nombre d'échantillons
Nbins = 100
xi_uniforme = [ uniform(-0.75,0.75) for i in range(M) ]
H,bin = Histogram(xi_uniforme,(-1.0,1.0),Nbins)
xlabel('x')
ylabel('H(x)')
title('Histogramme des valeurs de $x$')
show()
In [4]:
# une deuxième version avec filtrage
def Histogram(Liste,(xmin, xmax),N=10):
Deltax = (xmax-xmin)/float(N)
# on filtre les valeurs dans l'intervalle
Liste_filtree = filter(lambda x: xmin <= x <= xmax, Liste)
H = [0]*N
for x in Liste_filtree:
ind = min(N-1,int((x-xmin)//Deltax))
H[ind] += 1
# abcisses du début des boîtes
bins = [xmin + n*Deltax for n in range(N+1)]
# abscisses pour l'histogramme avec dédoublement
X = sorted(bins + bins[1:-1])
# ordonnées pour l'histogramme avec normalisation
# tient compte des données rejetées
norm = Deltax*float(sum(H)) * len(Liste)/float(len(Liste_filtree))
for i in range(len(H)): H[i] /= norm
# dédoublement des indices pour les ordonnées
Y = [ H[i] for i in sorted(range(len(H))*2) ]
plot(X,Y,'b')
return (H,bins)
In [5]:
# génération de liste aléatoires
from random import gauss
M = 10000 # nombre d'échantillons
Nbins = 20
sigma = 0.2
xi_gauss = [ gauss(0.0,sigma) for i in range(M) ]
H,bin = Histogram(xi_gauss,(-1.0,1.0),Nbins)
print H
print bin
# expression analytique
y = [ 1./sqrt(2*pi*sigma**2)*exp(-0.5*(x/sigma)**2) for x in bin ]
plot(bin,y,'r-')
xlabel('x')
ylabel('H(x)')
title('Histogramme des valeurs de $x$')
show()
In [6]:
help(hist)
In [7]:
# avec la fonction hist de pylab, attention à bien donner range=() et normed=1
n2, bin2, patches = hist(xi_gauss, Nbins, range=(-1.0,1.0), normed=1, facecolor='white')
print n2
print bin2
xlabel('x')
ylabel('H(x)')
title('Histogramme des valeurs de $x$')
y = [ 1./sqrt(2*pi*sigma**2)*exp(-0.5*(x/sigma)**2) for x in bin ]
plot(bin,y,'r-')
show()
La pommme d'arrosage¶
In [8]:
# Algorithme de Verlet
def Verlet(r0,v0,tf,F,N=100):
Fx, Fz = F
dt = tf/float(N)
x, z = [0.0]*N, [0.0]*N
x[0], z[0] = r0[0], r0[1]
x[1] = x[0] + v0[0]*dt + 0.5*Fx(x[0],z[0])*dt**2
z[1] = z[0] + v0[1]*dt + 0.5*Fz(x[0],z[0])*dt**2
for i in range(2,N):
x[i] = 2*x[i-1]-x[i-2] + Fx(x[i-1],z[i-1])*dt**2
z[i] = 2*z[i-1]-z[i-2] + Fz(x[i-1],z[i-1])*dt**2
return (x,z)
In [9]:
from math import cos, sin, tan, pi
g, v0, h = 9.81, 0.2, 1.0
def Fx(x,z): return 0.0
def Fz(x,z): return -g
N, tf = 100, 0.5
for theta in [pi/5,pi/4,pi/3]:
(x,z) = Verlet((0.0,h),(v0*sin(theta),v0*cos(theta)),tf,(Fx,Fz),N)
plot(x, z)
def Exact(x): return h + x/tan(theta) -0.5*g*(x/(v0*sin(theta)))**2
plot(x, [ Exact(pos) for pos in x ],'+')
show()
In [10]:
theta = pi/4
(x,z) = Verlet((0.0,h),(v0*sin(theta),v0*cos(theta)),tf,(Fx,Fz),N)
def Exact(x): return h + x/tan(theta) -0.5*g*(x/(v0*sin(theta)))**2
diff = [ abs(Exact(x[i])-z[i])**0.5 for i in range(len(x)) ]
plot(diff)
show()
In [12]:
# une fonction qui détermine la position au sol à partir de (x,z)
def Position(x,z):
for i in range(len(z)):
if z[i] < 0.0:
break
# return x[i] # version sans interpolation assez sommaire
return x[i] - z[i]*(x[i]-x[i-1])/(z[i]-z[i-1]) # avec interpolation linéaire
Npoints, thetamin, thetamax = 1000, pi/5, pi/3
dtheta = (thetamax-thetamin)/float(Npoints-1)
Thetas = [thetamin + i*dtheta for i in range(Npoints) ]
pos = []
for theta in Thetas:
(x,z) = Verlet((0.0,h),(v0*sin(theta),v0*cos(theta)),tf,(Fx,Fz),N)
pos.append(Position(x,z))
In [13]:
n, bin, patches = hist(pos, bins=100, range=(0.05,0.085), normed=1, facecolor='white')
xlabel('x')
ylabel('H(x)')
title('Histogramme des valeurs de $x$')
show()
In [ ]: