[Wissenschaftlich-technische Berechnung von Python] Anpassung durch nichtlineare Funktion, Zustandsgleichung, scipy
Mit der Curvefit-Methode von scipy.optimize werden XY-Daten mit einer nichtlinearen Funktion versehen, um den optimalen Anpassungsparameter zu erhalten.
Zweck: Bereiten Sie die Daten zu Materialvolumen (V) und Druck (P) als XY-Daten vor. Es drückt die Beziehung zwischen dem Volumen und dem Druck einer Festkörpergleichung aus (https://ja.wikipedia.org/wiki/%E7%8A%B6%E6%85%8B%E6%96%B9%E7%A8) Weit verbreitet als eine von% 8B% E5% BC% 8F_ (% E7% 86% B1% E5% 8A% 9B% E5% AD% A6)) Birken-Managan-Zustandsgleichung 3. Ordnung Finden Sie die Anpassungsparameter durch Anpassen an (: //en.wikipedia.org/wiki/Birch%E2%80%93Murnaghan_equation_of_state). Diese Gleichung hat die folgende Form:
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 $ sind die Anpassungsparameter. Die X-Achse ist V und die y-Achse ist P.
** Je nach Problem kann der optimale Anpassungsparameter für jede Funktionsform erhalten werden, indem die Spezifikation der Anpassungsfunktionsform im Code geändert wird. Natürlich ist auch eine lineare Anpassung möglich. Es ist ein vielseitiger Code, außer in besonderen Situationen. ** **.
import numpy as np
import scipy.optimize
import matplotlib.pylab as plt
"""
Kurvenanpassung mit nichtlinearer Funktion
Beispiel:Passen Sie die Beziehung zwischen Druck und Volumen an eine kubische Birch-Managan-Gleichung an
"""
x_list=[] # x_Liste definieren(Erstellen Sie eine leere Liste)
y_list=[] # y_Liste definieren(Erstellen Sie eine leere Liste)
f=open('EoS.dat','rt') #R die Datei, die die Daten enthält, die Sie zeichnen möchten(Lesen) t(Text)Einlesemodus
##Lesen Sie die Daten, x_Liste und y_Speichern Sie den Wert in der 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')
#Anpassen der Anpassungsfunktion
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] #Ausgangsbedingungen für Anpassungsparameter
para_opt, cov = scipy.optimize.curve_fit(fitfunc, x_list, y_list, para_ini)#Anpassungsparameter optimieren
#Anzeige der erhaltenen Anpassungsparameter
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()
Ergebnis
Fitted V0 = 163.627068001 (Å^3) Fitted K0 = 225.52189108 (GPa) Fitted K0' = 4.06160382625
[Anhang] Inhalt der verwendeten EoS.dat-Datei
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
Recommended Posts