Générateur de nombres aléatoires

In [6]:
# première version
def aleatoire(un):
    return (a*un+c)%m

def affiche(graine,n=10):
    rand = graine
    for i in range(n):
        rand = aleatoire(rand)
        print(rand,end=" ")
    print()
    
# wikipedia non-optimal
a, c, m = 25, 16, 256
affiche(10), affiche(50), affiche(125)
print()
# quickran choice
a, c, m = 8121, 28411, 134456
affiche(10), affiche(50), affiche(125)
print()
# standard minimum
a, c, m = 16807, 0, 2**31-1
affiche(10), affiche(50), affiche(125)

10 10 10 10 10 10 10 10 10 10 
242 178 114 50 242 178 114 50 242 178 
69 205 21 29 229 109 181 189 133 13 

109621 27376 93139 95230 329 11100 85991 131314 58869 112480 
31093 26296 62099 124390 31673 30516 46439 10450 51125 14408 
102344 91499 87934 44609 73636 100535 56314 69549 120640 99435 

168070 677268843 1194115201 1259501992 703671065 407145426 1010275440 1693606898 1702877348 745024267 
840350 1238860568 1675608711 2002542666 1370871678 2035727130 756409906 2025583549 2071935799 1577637688 
2100875 949667773 967796307 711389371 1279695548 794350531 1891024765 1842733402 1958614027 1796610573
Out[6]:
(None, None, None)
In [7]:
# deuxième version avec le mot-clé global
graine = 10213
a, c, m = 16807, 0, 2**31-1
u = graine
def aleatoire():
    global u  # nécessaire pour utiliser le u défini au-dessus et non une variable 'u' locale
    last, u = u, (a*u+c)%m
    return last

print ([ aleatoire() for i in range(10) ])

[10213, 171649891, 849180116, 2141375297, 416176606, 325978763, 494286244, 1002156312, 526892363, 1404868360]
In [8]:
# distribution uniforme

def Uniform(xmin,xmax):
    return (xmax-xmin)*aleatoire()/float(m-1)+xmin

print ([ Uniform(-1,1) for i in range(10) ])

[-0.9629075126377005, 0.41343480852752434, 0.5988158607844412, 0.29815969178281687, -0.8300703659933716, -0.9926425805246853, 0.6561490648017723, -0.10268083783116266, 0.2431515494763401, 0.6480823202506474]
In [10]:
reels = [ Uniform(-1,1) for i in range(1000) ]
# ou 
from random import uniform
reels = [ uniform(-1,1) for i in range(1000) ]
print(reels[:10])

[-0.20115895151860075, -0.34505410414328197, -0.241843647116198, -0.3739170177084783, -0.027002185534297007, 0.8496005459543214, -0.6132072729380802, 0.9707878174903044, -0.7392885381144507, -0.7973081935837594]
In [13]:
moy = sum(reels)/len(reels)
print ("Moyenne : {:.10f}".format(moy))

Moyenne : 0.0038859770
In [15]:
from math import sqrt
stdev = sqrt(sum([ (x-moy)**2 for x in reels ])/len(reels))
print ("Ecart-type : {:.10f}".format(stdev))

Ecart-type : 0.5914029445
In [17]:
import numpy as np
print(np.mean(reels), np.std(reels))

0.00388597698314 0.59140294449

Simulation du tir ballistique

In [18]:
# 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 [19]:
from pylab import *
%matplotlib inline

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 [20]:
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 [21]:
# 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 [22]:
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 [ ]: