# CAPACITE 1 : MARCHE AU HASARD A TROIS DIMENSIONS ################################################################################ # N particules sont placées à l'origine O à t = 0. # La diffusion est modélisée par une marche au hasard avec un libre parcours moyen l et une durée tau entre chaque pas. # A chaque nouveau pas, la particule k part dans une direction aléatoire d'angle theta par rapport à l'axe Oz, d'angle phi autour de l'axe Oz et se déplace de l. # Toutes les valeurs de theta sur [0,pi] et de phi sur [0,2pi] sont équiprobables import numpy as np import matplotlib.pyplot as plt from mpl_toolkits import mplot3d import random import time N = 500 # Nombre de particules l = 0.75 # Libre parcours moyen tau = 1 # Durée entre deux collisions (durée d'un pas) P = 1000 # Nombre de pas X = np.zeros(N) # Tableau des N abscisses initiales nulles Y = np.zeros(N) # Tableau des N ordonnées initiales nulles Z = np.zeros(N) # Tableau des N cotes initiales nulles r = np.zeros(N) # Tableau des N distances au centre r= sqrt(X²+Y²+Z²) initiales nulles # Après chaque pas, on calcule la distance quadratique moyenne au centre O : rqm = sqrt() (< > moyenne sur les N particules) # rqm va mesurer la distance caractéristique L de la diffusion des particules. # A l'instant t, la théorie de la marche au hasard à trois dimensions montre que rqm = sqrt(6Dt) où D est le coefficent de # diffusion donné par D = l²/6tau . # Donc : L² = rqm² = 6Dt = l²t/tau # On se propose de vérifier cette loi et d'en déduire D. MoyX = np.zeros(P+1) # Initialisation du tableau des P+1 moyennes des N abscisses (1 valeur initiale et P valeurs après chaque pas) MoyY = np.zeros(P+1) # Initialisation du tableau des P+1 moyennes des N ordonnées MoyZ = np.zeros(P+1) # Initialisation du tableau des P+1 moyennes des N cotes rqm = np.zeros(P+1) # Initialisation du tableau des P+1 distances quadratiques moyennes après chaque pas ################################################################################ # Animation de la diffusion des N particules dans l'espace Oxyz start = time.time() # Lancement d'un chronomètre pour mesurer le temps de calcul for i in range(1,P+1): # Pour chaque nouveau pas (P valeurs de i allant de 1 à P) for k in range(N): # Pour chaque particule theta = random.uniform(0,np.pi) # Tirage aléatoire et équiprobable de theta entre 0 et pi phi = random.uniform(0,2*np.pi) # Tirage aléatoire et équiprobable de phi entre 0 et 2pi # Le déplacement de la particule k pendant le pas i est alors et l.sin(theta).cos(phi) selon Ox, l.sin(theta).sin(phi) selon Oy et l.cos(theta) selon Oz X[k] = X[k] + l*np.sin(theta)*np.cos(phi) # Nouvelle abscisse de la particule Y[k] = Y[k] + l*np.sin(theta)*np.sin(phi) # Nouvelle ordonnée de la particule Z[k] = Z[k] + l*np.cos(theta) # Nouvelle cote de la particule r[k] = np.sqrt(X[k]**2+Y[k]**2+Z[k]**2) # Distance au centre r = sqrt(X²+Y²+Z²) MoyX[i] = np.mean(X) # Moyenne des abscisses des N particules après le pas i MoyY[i] = np.mean(Y) # Moyenne des ordonnées des N particules après le pas i MoyZ[i] = np.mean(Z) # Moyenne des cotes des N particules après le pas i rqm[i] = np.sqrt(np.mean(r**2)) # Distance quadratique moyenne au centre O des N particules après le pas i plt.clf() # Effacement de l'image précédente plt.grid() ax = plt.axes(projection='3d') ax.set_xlim3d(-50,50) ax.set_ylim3d(-50,50) ax.set_zlim3d(-50,50) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') ax.plot3D(X,Y,Z,'o',markersize=1) plt.pause(0.01) # Durée entre deux images plt.show() end = time.time() # Fin du chronomètre Tcalcul = end - start # Temps de calcul print('Le temps de calcul est :',Tcalcul) ################################################################################ # Histogrammes des abcisses et ordonnées des N particules après les P pas plt.subplot(1,3,1) plt.title("Histogramme de l'abscisse x") plt.xlabel('x') plt.hist(X,bins='rice') # Histogramme des abscisses après les P pas plt.subplot(1,3,2) plt.title("Histogramme de l'ordonnée y") plt.xlabel('y') plt.hist(Y,bins='rice') # Histogramme des ordonnées après les P pas plt.subplot(1,3,3) plt.title("Histogramme de la cote z") plt.xlabel('z') plt.hist(Z,bins='rice') # Histogramme des cotes après les P pas plt.show() ################################################################################ # Moyennnes des abscisses et des ordonnées, distance quadratique moyenne au centre en fonction du temps, t = np.arange(0,(P+1)*tau,tau) # Tableau des temps [0,tau,2tau,3tau,....,Ptau] plt.subplot(2,2,1) plt.title("Moyenne des abscisses au cours du temps") plt.xlabel('temps') plt.ylabel('x moyen') plt.plot(t,MoyX) plt.grid() plt.subplot(2,2,2) plt.title("Moyenne des ordonnées au cours du temps") plt.xlabel('temps') plt.ylabel('y moyen') plt.plot(t,MoyY) plt.grid() plt.subplot(2,2,3) plt.title("Moyenne des cotes au cours du temps") plt.xlabel('temps') plt.ylabel('z moyen') plt.plot(t,MoyZ) plt.grid() plt.subplot(2,2,4) plt.title("Distance quadratique moyenne au cours du temps") plt.xlabel('temps') plt.ylabel('distance quadratique moyenne') plt.plot(t,rqm) plt.grid() plt.show() ################################################################################ # Vérification de la loi L² = 6Dt pour un phénomène de diffusion à trois dimensions L = rqm # Définition de la distance caractéristique L de diffusion, donc d'étalement du nuage de points plt.plot(t,L**2,'+r',markersize=1) p = np.polyfit(t,L**2,1) # Regression linéaire de L² = f(t) plt.plot(t,p[0]*t+p[1],'k--',linewidth=1,label='Modèle affine') # Tracé de la droite modèle L² = at + b plt.title('L² en fonction du temps') plt.xlabel('temps') plt.ylabel('L²') plt.legend() plt.grid() print('Le coefficient de diffusion est :',p[0]/6) D = l**2/(6*tau) # Valeur théorique du coefficient de diffusion print('La valeur théorique du coefficient de diffusion est :',D) plt.show() ################################################################################