# -*- coding: utf-8 -*- """ Created on Wed Oct 19 11:14:53 2016 @author: adomps """ import numpy as np import matplotlib.pyplot as plt import random as rd import time # Attention # rd.randint(a, b) : a et b inclus # np.random.randint(a, n) : de a jusqu'à b EXCLU ! # Les particules marchent depuis l'origine sur un réseau carré. # mouvement sur réseau avec 4 directions possibles NSEW = [np.array([0, 1]), np.array([-1, 0]), np.array([1, 0]) , np.array([0, -1])] D = 1**2 / 4 / 1. # Déplacement d'une particule. Version très naïve def marche_alea_2d(N_pas, option) : trajectoire = np.zeros([N_pas+1, 2]) if option == 0 : for i in range(N_pas) : trajectoire[i+1] = trajectoire[i] + rd.choice(NSEW) else : # plus rapide : on définit d'abord les déplacements puis son les #somme avec cumsum deplacements = np.zeros([N_pas, 2]) if option == 1 : # déplacements tirés itérativement for i in range(0, N_pas) : deplacements[i] = rd.choice(NSEW) elif option == 2 : # déplacements tirés en groupe, encore plus rapide. dep_x = np.random.choice( (-1, 1, 0, 0), size = N_pas) dep_y = np.random.choice((-1, 1), size = N_pas) deplacements[:, 0] = dep_x masque = deplacements[:, 0] == 0 # on cherche les cases pour lesquelles Delta x = 0 deplacements[masque, 1] = dep_y[masque] # pour ces cas, Delta y = +1 ou -1 trajectoire[1:] = np.cumsum(deplacements, axis = 0) return trajectoire def trace_marche_2d(tj, numfig) : plt.figure(numfig) plt.plot(tj[:,0], tj[:,1]) plt.xlabel('x') plt.ylabel('y') N_pas = 10000 t1 = time.time() tj0 = marche_alea_2d(N_pas, 0) t2 = time.time() print('durée', t2 - t1) t1 = time.time() tj1 = marche_alea_2d(N_pas, 1) t2 = time.time() print('durée', t2 - t1) t1 = time.time() tj2 = marche_alea_2d(N_pas, 2) t2 = time.time() print('durée', t2 - t1) plt.figure(0) plt.clf() trace_marche_2d(tj0, 0) trace_marche_2d(tj1, 0) trace_marche_2d(tj2, 0) # marche d'un ensemble de particules # On stocke à tous les instants les positions de toutes les particules. # On utilise pour cela un tableau positions de taille N_pas x N_part x 2. # positions[i, j, k]. i : instant. j : particule. k = 0/1 pour x ou y def mouvement_ensemble_2d(N_pas, N_part) : ''' Renvoie un tableau de taille N_pas x N_part x 2. position[i, j, k] : i désigne l'instant, j la particule, k = 0/1 pour x ou y. ''' positions = np.zeros((N_pas+1, N_part, 2)) for j in range(N_part) : positions[:, j, :] = marche_alea_2d(N_pas, 2) return positions N_pas = 1000 N_part = 1000 t1 = time.time() Z = mouvement_ensemble_2d(N_pas, N_part) t2 = time.time() print('durée', t2 - t1) ### Représentation des particules à un instant particulier def dessine_ensemble(positions, instant, numfig) : X = positions[instant, : , 0] Y = positions[instant, : , 1] plt.xlabel('x') plt.ylabel('y') plt.figure(numfig) plt.scatter(X, Y, s = 10, marker = 'o', alpha = 1) # s : taille dates_choisies = [0, N_pas // 5, 2 * N_pas // 5, 3*N_pas //5, 4 * N_pas // 5, N_pas] x_max = 3 * np.sqrt(N_pas * D) plt.figure(1) for i, t_choisi in enumerate(dates_choisies) : plt.subplot(3, 2, i+1) plt.xlim(-x_max, x_max) plt.ylim(-x_max, x_max) ax = plt.gca() ax.set_box_aspect(1) dessine_ensemble(Z, t_choisi, 1) #### Croissance du rayon carré moyen r2_particules = np.sum(Z**2, axis = 2) r2_moyen = np.average(r2_particules, axis = 1) plt.figure(3) plt.clf() t = np.array([i for i in range(N_pas+1)]) plt.plot(t, r2_moyen) plt.plot(t, 4 * D * t, color = 'red') def distribution_particules_2d(positions, instant, N_boites, d_max, numfig) : X = positions[instant, : , 0] Y = positions[instant, : , 1] x_min, x_max = - d_max, d_max Delta_x = (x_max - x_min) / N_boites bords = [ x_min + i * Delta_x for i in range(N_boites + 1)] plt.figure(numfig) plt.clf() plt.hist2d(X, Y, bins = (bords, bords), density=False, cmap = 'Greys') plt.colorbar() popu, bords_x, bords_y = np.histogram2d(X, Y, (bords, bords), density=False) #x_centres = (bords_x[1:] + bords_x[0:-1]) / 2 #y_centres = (bords_y[1:] + bords_y[0:-1]) / 2 return popu # x_centres # y_centres N_part = 200 Z = mouvement_ensemble_2d(N_pas, N_part) t_choisi = N_pas dessine_ensemble(Z, t_choisi, 0) d_max = 4 * np.sqrt(D * t_choisi) # ordre de grandeur N_boites = 7 distribution_particules_2d(Z, t_choisi, N_boites, d_max, 1) ## Évolution en plusieurs images N_part = 10000 N_pas = 10000 N_boites = 21 Z = mouvement_ensemble_2d(N_pas, N_part) dates_choisies = [0, N_pas // 5, 2 * N_pas // 5, 3*N_pas //5, 4 * N_pas // 5, N_pas] d_max = 3 * np.sqrt(N_pas * D) plt.figure(3) for i, t_choisi in enumerate(dates_choisies) : #plt.subplot(3, 2, i+1) Je n'arrive pas à faire des subplot avec plt.hist plt.xlim(-d_max, d_max) plt.ylim(-d_max, d_max) ax = plt.gca() ax.set_box_aspect(1) # axes normés distribution_particules_2d(Z, t_choisi, N_boites, d_max, 3+i) ############################################################## #### Répartition radiale ############################################################# def repartition_radiale(positions, instant, N_classes, numfig): r2_particules = np.sum(positions[instant]**2, axis = 1) r_particules = np.sqrt(r2_particules) plt.hist(r_particules, bins = N_classes) popu_radiale, bords_r = np.histogram(r_particules, bins = N_classes, density = False) return popu_radiale, bords_r def trace_popu_radiale_theorique(bords_r, instant) : r_centres = (bords_r[1:] + bords_r[0:-1]) / 2 n_theo = N_part / (4 * np.pi * D * instant) * np.exp(- r_centres**2 / (4 * D * instant)) surfaces_zones = np.pi * bords_r[1:]**2 - np.pi * bords_r[:-1]**2 plt.plot(r_centres, n_theo * surfaces_zones, color = 'red') plt.xlabel('r') plt.text(bords_r[-5], 200, 't = ' + str(instant)) # N_boites = 20 plt.figure(2) plt.clf() popu_r, bords_r = repartition_radiale(Z, t_choisi, N_boites, 2) trace_popu_radiale_theorique(bords_r, t_choisi) # On peut inversement tracer n en divisant les populations par l'aire de chaque zone def trace_densite_radiale_theo(bords_r, instant) : r_centres = (bords_r[1:] + bords_r[0:-1]) / 2 n_theo = N_part / (4 * np.pi * D * instant) * np.exp(- r_centres**2 / (4 * D * instant)) plt.plot(r_centres, n_theo ) return def trace_densite_radiale_expe(bords_r, popu_radiale) : surfaces_zones = np.pi * bords_r[1:]**2 - np.pi * bords_r[:-1]**2 densite = popu_radiale / surfaces_zones r_centres = (bords_r[1:] + bords_r[0:-1]) / 2 plt.plot(r_centres, densite, 'o') return plt.figure(9) plt.clf() trace_densite_radiale_theo(bords_r, t_choisi) trace_densite_radiale_expe(bords_r, popu_r) plt.xlabel('r') plt.ylabel('n')