[Calcul scientifique / technique par Python] Ajustement par fonction non linéaire, équation d'état, scipy

En utilisant la méthode curvefit de scipy.optimize, les données XY sont équipées d'une fonction non linéaire pour obtenir le paramètre d'ajustement optimal.

Objectif: préparer les données de volume (V) et de pression (P) du matériau sous forme de données XY. Il exprime la relation entre le volume et la pression d'une [équation d'état] solide (https://ja.wikipedia.org/wiki/%E7%8A%B6%E6%85%8B%E6%96%B9%E7%A8 Largement utilisé comme l'un des% 8B% E5% BC% 8F_ (% E7% 86% B1% E5% 8A% 9B% E5% AD% A6)) équation de l'état de Birch-Managan de 3e ordre Trouvez les paramètres d'ajustement en ajustant à (: //en.wikipedia.org/wiki/Birch%E2%80%93Murnaghan_equation_of_state). Cette équation a la forme suivante:


P(V)=\frac{3}{2}K_0[(\frac{V}{V_0})^{-7/3}-(\frac{V}{V_0})^{-5/3}][1+\frac{3}{4}(K_d-4)((\frac{V}{V_0})^{-2/3}-1)]

$ K_0, V_0, K_d $ sont les paramètres d'ajustement. L'axe X est V et l'axe y est P.

** Selon le problème, le paramètre d'ajustement optimal peut être obtenu pour n'importe quelle forme fonctionnelle en modifiant la spécification de la forme de fonction d'ajustement dans le code. Bien entendu, un ajustement linéaire est également possible. C'est un code polyvalent sauf dans des situations particulières. ** **

import numpy as np
import scipy.optimize
import matplotlib.pylab as plt

"""
Ajustement de courbe avec fonction non linéaire
Exemple:Ajuster la relation pression / volume avec une équation cubique de Birch-Managan
"""



x_list=[] # x_définir la liste(Créer une liste vide)
y_list=[] # y_définir la liste(Créer une liste vide)

f=open('EoS.dat','rt') #R le fichier contenant les données que vous souhaitez tracer(Lis) t(texte)Lire en mode

##Lisez les données, x_list et y_Stocker la valeur dans la liste
for line in f:
    data = line[:-1].split(' ')
    x_list.append(float(data[0]))
    y_list.append(float(data[1]))
##
plt.plot(x_list, y_list, 'o',label='Raw data')

#Spécification de la fonction d'adaptation
def fitfunc(V, V0, K0, K0d ):
    return (3.0/2.0)*K0*((V/V0)**(-7.0/3.0)-(V/V0)**(-5.0/3.0))*(1.0+(3.0/4.0)*(K0d-4.0)*((V/V0)**(-2.0/3.0)-1.0))


para_ini =[200.0, 100.0, 4.0] #Conditions initiales des paramètres d'ajustement


para_opt, cov = scipy.optimize.curve_fit(fitfunc, x_list, y_list, para_ini)#Optimiser les paramètres d'ajustement

#Affichage des paramètres d'ajustement obtenus
print ("Fitted V0 =", str(para_opt[0]), "(Å^3)")
print ("Fitted K0 =", str(para_opt[1]), "(GPa)")
print ("Fitted K0' =", str(para_opt[2]))

##
# for plot
V0=para_opt[0]
K0=para_opt[1]
K0d=para_opt[2]

V1 = np.arange(min(x_list) - 1, max(x_list) + 1, 0.1)
y_fit = (3.0/2.0)*K0*((V1/V0)**(-7/3)-(V1/V0)**(-5.0/3.0))*(1.0+(3.0/4.0)*(K0d-4.0)*((V1/V0)**(-2.0/3.0)-1.0))

plt.plot(V1, y_fit, color='RED',linewidth=2.0, label='Fitted curve')

plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.grid(True)
plt.title("Curve fitting for eqaution of states")
plt.legend(loc='upper right')
plt.xlabel('Volume (Å^3)', fontsize=18)
plt.ylabel('Pressure (GPa)', fontsize=18)

plt.show()

résultat

Fitted V0 = 163.627068001 (Å^3) Fitted K0 = 225.52189108 (GPa) Fitted K0' = 4.06160382625

[Annexe] Contenu du fichier EoS.dat utilisé

164.26 0.0 
156.05 11.9 
147.83 27.9 
146.19 31.7 
144.55 35.7 
142.91 39.9 
141.26 44.4 
139.62 49.2 
137.98 54.2 
136.34 59.5 
134.69 65.1 
133.05 71.0 
131.41 77.3 
129.77 83.9 
128.12 90.9 
126.48 98.3 
124.84 106.2
123.20 114.4
121.55 123.2
119.91 132.5
118.27 142.3
116.62 152.6
114.98 163.6
113.34 175.2
111.70 187.5
110.05 200.6