[Scientific / technical calculation by Python] Fitting by nonlinear function, equation of state, scipy

Using the curvefit method of scipy.optimize, the XY data is fitted with a non-linear function to obtain the optimum fitting parameter.

Purpose: Prepare the volume (V) and pressure (P) data of the substance as XY data. It expresses the relationship between the volume and pressure of a solid [Equation of state](https://ja.wikipedia.org/wiki/%E7%8A%B6%E6%85%8B%E6%96%B9%E7%A8 Widely used as one of the% 8B% E5% BC% 8F_ (% E7% 86% B1% E5% 8A% 9B% E5% AD% A6)) 3rd order Birch-Marnagan equation of state Find the fitting parameters by fitting to (: //en.wikipedia.org/wiki/Birch%E2%80%93Murnaghan_equation_of_state). This equation has the following 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 $ are the fitting parameters. The X-axis is V and the y-axis is P.

** Depending on the problem, the optimal fit parameter can be obtained for any functional form by changing the specification of the fit function form in the code. Of course, a linear fit is also possible. It is a versatile code except in special situations. ** **

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

"""
Curve fitting with non-linear function
Example:Fitting the pressure vs. volume relationship with a cubic Birch-Managan equation
"""



x_list=[] # x_define list(Create an empty list)
y_list=[] # y_define list(Create an empty list)

f=open('EoS.dat','rt') #R the file that contains the data you want to plot(Read) t(text)Read in mode

##Read the data, x_list and y_Store the value in list
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')

#Specifying fitting function
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] #Initial conditions for fitting parameters


para_opt, cov = scipy.optimize.curve_fit(fitfunc, x_list, y_list, para_ini)#Optimizing fitting parameters

#Display of the obtained fitting parameters
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()

result

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

[Appendix] Contents of the used EoS.dat file

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