# -*- 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 # 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])] # Valeur attendue du coefficient de diffusion 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([...]) if option == 0 : for i in range(N_pas) : .... else : # plus rapide : on définit d'abord les déplacements puis on les somme deplacements = np.zeros([...]) if option == 1 : # déplacements tirés itérativement for i in range(0, N_pas) : deplacements[i] = .... elif option == 2 : # déplacements tirés en groupe, encore plus rapide. dep_x = np.random.choice( [...], size = N_pas) # attention aux proba 1/4, 1/4, 1/2 dep_y = np.random.choice([...], size = N_pas) # on prend seulemen +1 ou -1 deplacements[:, 0] = ... # déplacements selon x # << masque >> est un tableau de booléens. Ses éléments True # identifient les cases pour lesquelles Delta x = 0 masque = deplacements[:, 0] == 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) # On ouvre un novelle figure plt.plot(..., ...) # toutes les abscisses, puis toutes les ordonnées plt.xlabel('x') plt.ylabel('y') ### Applications 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((....)) # Donner les tailles du tableau for j in range(N_part) : positions[:, j, :] = .... # trajectoire de la particule d'indice j return positions ## Application N_pas = 1000 N_part = 100 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 #### Exemples 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 ... ... ... ... ######################################################################### def distribution_particules_2d(positions, instant, N_boites, d_max, numfig) : X = ... # toutes les abscisses à cet instant Y = ... # toutes les ordonnées à cet instant 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(..., ... , 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 ## Application N_part = 2000 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.xlim(-x_max, x_max) plt.ylim(-x_max, x_max) ax = plt.gca() ax.set_box_aspect(1) 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 = .... # carrés des distances de toues les particules à l'origine à cet instant r_particules = np.sqrt(r2_particules) plt.hist(..., bins = N_classes) popu_radiale, bords_r = np.histogram(... , 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 = .... popu_theorique = .... plt.plot(r_centres, popu_theorique, 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) plt.ylabel('population')