Ce document fait suite à Modèle des sphères dures à deux dimensions, où est présenté l'algorithme de Metropolis appliqué au modèle des sphères dures.
On s'intéresse ici à l'obtention de la pression en fonction de la densité, en suivant la méthode utilisée dans [1].
Il s'agit de calculer la fonction de distribution radiale g(r) pour des valeurs proches de r=d (diamètre des sphères), puis d'obtenir g(d) par extrapolation.
import sys sys.path.append("../spheres2d") from spheres2dMetropolis import * from matplotlib.pyplot import *
On calcule la fonction de distribution jusqu'à r=2d :
N = 7 densite = 0.4 systeme = Systeme(N,densite) systeme.initialiser() systeme.boucle_metropolis(100000) figure(figsize=(8,6)) (r,g) = systeme.calculer_distribution_radiale(3,20,10000,100) plot(r,g,'o') xlabel('r') ylabel('g') title('densite = %f'%densite) axis([0,3,0,g.max()])Figure pleine page
Pour obtenir la valeur extrapolée de g(d), on utilise une interpolation polynomiale avec l'algorithme de Neville. La fonction suivante fait ce calcul, en prenant les P premiers points :
def extrapoler(r,g,P): def recurrence(i,m,x): if m==0: return g[i] else: return ((x-r[i+m])*recurrence(i,m-1,x)+(r[i]-x)*recurrence(i+1,m-1,x))/(r[i]-r[i+m]) return recurrence(0,P-1,1.0)
print(extrapoler(r,g,5)) --> 2.2831652360918731
L'équation d'état est donnée par la relation suivante :
On trace donc ce rapport en fonction de la densité (aire occupée par les sphères divisée par l'aire totale).
import math densite = [0.1,0.2,0.3,0.4,0.5] pression = [] for k in range(len(densite)): systeme = Systeme(N,densite[k]) systeme.initialiser() systeme.boucle_metropolis(100000) (r,g) = systeme.calculer_distribution_radiale(2,20,10000,100) pression.append(1.0+math.pi/2*systeme.N/(systeme.L**2)*extrapoler(r,g,5)) figure(figsize=(8,6)) plot(densite,pression,'o') xlabel('densite') ylabel('PA/NkT') axis([0,0.5,0.0,max(pression)])Figure pleine page