Outils pour utilisateurs

Outils du site


teaching:progappchim:fit_modele_einstein

Différences

Ci-dessous, les différences entre deux révisions de la page.

Lien vers cette vue comparative

teaching:progappchim:fit_modele_einstein [2014/03/14 15:06]
villersd créée
teaching:progappchim:fit_modele_einstein [2015/04/01 15:47] (Version actuelle)
villersd
Ligne 7: Ligne 7:
 Le programme Python nécessite les librairies scipy (optimisation),​ numpy (manipulation des données) et matplotlib (représentation du fit). Le programme Python nécessite les librairies scipy (optimisation),​ numpy (manipulation des données) et matplotlib (représentation du fit).
  
-<sxh python; title : fit-Cv-diamant-Einstein-03.py>+<sxh python; title : fit-Cv-diamant-Einstein-04.py>
 #​!/​usr/​bin/​env python #​!/​usr/​bin/​env python
 # -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
Ligne 13: Ligne 13:
 Fit des données (température absolue, chaleur spécifique molaire à volume constant Fit des données (température absolue, chaleur spécifique molaire à volume constant
 Diamant Diamant
-Modèle d'​Einstein+Modèle d'​Einstein ​: http://​fr.wikipedia.org/​wiki/​Mod%C3%A8le_d%27Einstein
 """​ """​
 import numpy as np import numpy as np
Ligne 20: Ligne 20:
 import matplotlib.pyplot as plt import matplotlib.pyplot as plt
  
-def cvEinstein(T,TE):+def cvEinstein(Tr):
     # fonction à fitter     # fonction à fitter
-    return 3.*NA*kB*(TE/T)**2*sp.exp(TE/T)/(sp.exp(TE/T)-1)**2+    ​# Tr = température réduite = T/TE 
 +    ​return 3.*NA*kB*(1./Tr)**2*sp.exp(1./Tr)/(sp.exp(1./Tr)-1)**2
  
 def residuals(p,​x,​y):​ def residuals(p,​x,​y):​
     # erreur entre y donné et y calculé     # erreur entre y donné et y calculé
-    return cvEinstein(x,p[0]) - y+    return cvEinstein(x/p[0]) - y
  
  
Ligne 71: Ligne 72:
  
 p=[1000.] # valeur initiale de l'​unique paramètre recherché (TE) p=[1000.] # valeur initiale de l'​unique paramètre recherché (TE)
-plsq = leastsq(residuals,​ p, args=(data[0],​ data[1])) ​ # fit par moinde ​carré+plsq = leastsq(residuals,​ p, args=(data[0],​ data[1])) ​ # fit par moindre ​carré
 # Résultat : # Résultat :
 print "​Température caractéristique suivant le modèle d'​Einstein:​ ",​plsq[0] print "​Température caractéristique suivant le modèle d'​Einstein:​ ",​plsq[0]
Ligne 77: Ligne 78:
 # graphe # graphe
 x=np.linspace(1.,​1.5*plsq[0],​200) x=np.linspace(1.,​1.5*plsq[0],​200)
-plt.plot(x,​cvEinstein(x,plsq[0]),​data[0],​ data[1],'​o'​)+plt.plot(x,​cvEinstein(x/plsq[0]),​data[0],​ data[1],'​o'​)
 plt.title(u'​Capacité calorifique du diamant : optimisation par moindre carré'​) plt.title(u'​Capacité calorifique du diamant : optimisation par moindre carré'​)
 #​plt.legend(['​Fit',​ '​Experimental'​]) #​plt.legend(['​Fit',​ '​Experimental'​])
Ligne 93: Ligne 94:
 {{:​teaching:​progappchim:​einstein-c-diamond-fit-lowt.png?​direct|}} {{:​teaching:​progappchim:​einstein-c-diamond-fit-lowt.png?​direct|}}
  
-Le modèle de Debye, qui utilise un spectre de fréquences plutôt qu'une fréquence unique de vibration, pourra remplacer le modèle d'​Einstein.+Le modèle de Debye, qui utilise un spectre de fréquences plutôt qu'une fréquence unique de vibration, pourra remplacer le modèle d'​Einstein. ​Le calcul nécessite d'​utiliser les fonctions d'​intégration scipy.integrate.quad,​ qui oblige de contourner l'​utilisation des tableau numpy : les valeurs d'un tableau à traiter sont considérées une par une et le tableau numpy est recréé immédiatement. Voici un programme : 
 + 
 +<sxh python; title : fit-Cv-diamant-Debye-02.py>​ 
 +#​!/​usr/​bin/​env python 
 +# -*- coding: utf-8 -*- 
 +"""​  
 +Fit des données (température absolue, chaleur spécifique molaire à volume constant 
 +Diamant 
 +Modèle de Debye : http://​fr.wikipedia.org/​wiki/​Mod%C3%A8le_de_Debye 
 +"""​ 
 +import numpy as np 
 +import scipy as sp 
 +from scipy.optimize import leastsq 
 +from scipy import integrate 
 +from scipy.integrate import simps 
 +import matplotlib.pyplot as plt 
 + 
 +def f(x): 
 +    return x**4.*sp.exp(x)/​(sp.exp(x)-1.)**2. 
 + 
 +def cvDebye(Tr):​ 
 +    # Chaleur spécifique théorique (modèle de Debye) 
 +    # Tr = température réduite = T/TD 
 +    # problème : integrate.quad n'​accepte pas le broadcast pour les limites a, b 
 +    # http://​newscentral.exsees.com/​item/​12e4126a47e3411a138f5a3820eb62d1-71278f2740f3cd9c3dff2da0d78b6bc7 
 +    #solution : 
 +    if '​array'​ in str(type(Tr)):​ 
 +        return 9.*NA*kB*(Tr)**3 *np.array([integrate.quad(f,​0,​1./​each_Tr)[0] for each_Tr in Tr])  # array 
 +    else: 
 +        return 9.*NA*kB*(Tr)**3 * integrate.quad(f,​0,​1./​Tr)[0] ​  # float 
 + 
 +def residuals(p,​x,​y):​ 
 +    # erreur entre y donné et y calculé 
 +    return cvDebye(x/​p[0]) - y 
 + 
 + 
 +kB = 1.3806488E-23 ​ # Constante de Boltzmann 
 +NA = 6.02214129E23 ​ # Nombre d'​Avogadro 
 +# données expérimentales (T, Cv) : 
 +liste_data=[ 
 +[12.9,​0.00053],​ 
 +[16.1,​0.00081],​ 
 +[19.8,​0.00138],​ 
 +[24.1,​0.00257],​ 
 +[30.1,​0.00494],​ 
 +[33.4,​0.0074],​ 
 +[41.3,​0.0133],​ 
 +[47.7,​0.02],​ 
 +[57.2,​0.0365],​ 
 +[67,​0.0595],​ 
 +[76.1,​0.092],​ 
 +[87,​0.147],​ 
 +[100.4,​0.24],​ 
 +[113.1,​0.378],​ 
 +[126.3,​0.56],​ 
 +[143.4,​0.88],​ 
 +[159,​1.19],​ 
 +[176,​1.66],​ 
 +[197,​2.21],​ 
 +[215,​2.61],​ 
 +[264,​4.18],​ 
 +[273,​5.2],​ 
 +[280,​5.43],​ 
 +[306,​6.59],​ 
 +[335,​7.81],​ 
 +[363,​9.02],​ 
 +[412,​11.8],​ 
 +[471,​14.6],​ 
 +[516,​15.6],​ 
 +[874,​22.3],​ 
 +[1079,​22.4],​ 
 +[1238,​22.7]] 
 + 
 +# transformation en tableau numpy 
 +data=np.array(zip(*liste_data)) 
 +# cf. http://​docs.python.org/​2.7/​library/​functions.html#​zip 
 + 
 +p=[1000.] # valeur initiale de l'​unique paramètre recherché (TE) 
 +plsq = leastsq(residuals,​ p, args=(data[0],​ data[1])) ​ # fit par moindre carré 
 +# Résultat : 
 +print "​Température caractéristique suivant le modèle de Debye: ",​plsq[0] 
 + 
 +# graphe 
 +x=np.linspace(1.,​1.5*plsq[0],​200) 
 +plt.plot(x,​cvDebye(x/​plsq[0]),​data[0],​ data[1],'​o'​) 
 +plt.title(u'​Capacité calorifique du diamant : optimisation par moindre carré'​) 
 +#​plt.legend(['​Fit',​ '​Experimental'​]) 
 +plt.show() 
 +</​sxh>​ 
 ===== Références ===== ===== Références =====
   * [[http://​en.wikipedia.org/​wiki/​Einstein_solid]]   * [[http://​en.wikipedia.org/​wiki/​Einstein_solid]]
teaching/progappchim/fit_modele_einstein.txt · Dernière modification: 2015/04/01 15:47 par villersd