曲线拟合 matplotlib 中函数的函数
Curve fitting a function of functions in matplotlib
我有 3 个函数：E（能量），P（压力）和 H（焓）。
给定数据 data.dat，包含变量 V（体积）和 E（能量）：
# Volume: V_C_I Energy: E_C_I
111.593876 -1.883070511360E+03
113.087568 -1.883074916825E+03
114.632273 -1.883078906679E+03
116.184030 -1.883082373429E+03
117.743646 -1.883085344696E+03
119.326853 -1.883087860954E+03
120.927806 -1.883089938181E+03
122.538335 -1.883091557526E+03
124.158641 -1.883092750745E+03
125.789192 -1.883093540824E+03
125.790261 -1.883093800747E+03
127.176327 -1.883094160364E+03
128.358654 -1.883094017730E+03
128.542807 -1.883094255789E+03
129.977279 -1.883094094751E+03
131.390610 -1.883093689121E+03
132.812287 -1.883093053342E+03
134.242765 -1.883092185844E+03
135.682211 -1.883091101112E+03
137.130792 -1.883089807766E+03
138.588565 -1.883088314435E+03
以下脚本执行：
1) E 的 curve_fit 与 V、
2) 使用def(P)函数计算P（压力），
3) 使用 def(H) 函数计算 H（焓）。 (H = E + PV).
4) 使用拟合曲线执行 3 个图：E 与 P、P 与 V 和 H 与 P .
绘制拟合曲线时，我使用了以下内容：
例如，对于 E 与 V 曲线：
plt.plot(V_C_I_lin, E(V_C_I_lin, *popt_C_I)
其中 V_C_I_lin 是数据点之间的体积的线性空间。
这看起来不错：
同样，对于 P 与 V 曲线，类似的方案：
plt.plot(V_C_I_lin, P(V_C_I_lin, *popt_C_I), color='grey', label='E fit Data' )
产生期望的结果：
脚本将每个数据点的压力保存在文件 ./E_V_P_H__C_I.dat 中。
然而，对于 H 与 P 曲线，类似的方案：
plt.plot(xp_C_I, H(V_C_I_lin, *popt_C_I), color='grey', label='H fit Data' )
其中 xp_C_I 是压力的线性空间，不会产生正确的拟合：
为什么会出现第三种情况？
代码：
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
import sys
import os
# Intial candidates for fit
E0_init = -941.510817926696
V0_init = 63.54960592453
B0_init = 76.3746233515232
B0_prime_init = 4.05340727164527
def E(V, E0, V0, B0, B0_prime):
return E0+ (2.293710449E+17)*(1E-21)*( (9.0/16.0)*(V0*B0) * ( (((V0/V)**(2.0/3.0)-1.0)**3.0)*B0_prime + ((V0/V)**(2.0/3.0)-1)**2 * (6.0-4.0*(V0/V)**(2.0/3.0)) ))
def P(V, E0, V0, B0, B0_prime):
f0=(3.0/2.0)*B0
f1=((V0/V)**(7.0/3.0))-((V0/V)**(5.0/3.0))
f2=((V0/V)**(2.0/3.0))-1
pressure= f0*f1*(1+(3.0/4.0)*(B0_prime-4)*f2)
return pressure
def H(V, E0, V0, B0, B0_prime):
return E(V, E0, V0, B0, B0_prime) + P(V, E0, V0, B0, B0_prime) * V
# Data (Red triangles):
V_not_p_f_unit_C_I, E_not_p_f_unit_C_I = np.loadtxt('data.dat', skiprows = 1).T
# A minor conversion of the data:
nFU_C_I = 2.0
nFU_C_II = 4.0
E_C_I = E_not_p_f_unit_C_I/nFU_C_I
V_C_I = V_not_p_f_unit_C_I/nFU_C_I
######
# Fitting and obtaining the parameters:
init_vals = [E0_init, V0_init, B0_init, B0_prime_init]
popt_C_I, pcov_C_I = curve_fit(E, V_C_I, E_C_I, p0=init_vals)
# Calculation of P:
pressures_per_F_unit_C_I = P(V_C_I, *popt_C_I)
# Calculation of H: H = E + PV
H_C_I = E_C_I + pressures_per_F_unit_C_I * V_C_I
# We save E, P and H into a file:
output_array_3 = np.vstack((E_C_I, V_C_I, pressures_per_F_unit_C_I, H_C_I)).T
np.savetxt('E_V_P_H__C_I.dat', output_array_3, header="Energy / FU (a.u.) \t Volume / FU (A^3) \t Pressure / F.U. (GPa) \t Enthalpy (a.u.)", fmt="%0.13f")
EnergyCI, VolumeCI, PressureCI, EnthalpyCI = np.loadtxt('./E_V_P_H__C_I.dat', skiprows = 1).T
# Plotting E vs V:
fig = plt.figure()
V_C_I_lin = np.linspace(VolumeCI[0], VolumeCI[-1], 100)
# Linspace for plotting the fitting curve:
V_C_I_lin = np.linspace(VolumeCI[0], VolumeCI[-1], 100)
p2, = plt.plot(V_C_I_lin, E(V_C_I_lin, *popt_C_I), color='grey', label='E fit Data' )
# Plotting the scattered points:
p1 = plt.scatter(VolumeCI, EnergyCI, color='red', marker="^", label='Data', s=100)
fontP = FontProperties()
fontP.set_size('small')
plt.legend((p1, p2 ), ("Data", 'E fit Data' ), prop=fontP)
plt.xlabel('V / Formula unit (Angstrom$^{3}$)')
plt.ylabel('E / Formula unit (a.u.)')
plt.ticklabel_format(useOffset=False)
plt.savefig('E_vs_V.pdf', bbox_inches='tight')
# Plotting P vs V:
fig = plt.figure()
# Linspace for plotting the fitting curve:
xp_C_I = np.linspace(PressureCI[-1], PressureCI[0], 100)
# Plotting the fitting curves:
p2, = plt.plot(V_C_I_lin, P(V_C_I_lin, *popt_C_I), color='grey', label='E fit Data' )
# Plotting the scattered points:
p1 = plt.scatter(VolumeCI, PressureCI, color='red', marker="^", label='Data', s=100)
fontP = FontProperties()
fontP.set_size('small')
plt.legend((p1, p2), ("Data", "E fit Data"), prop=fontP)
plt.xlabel('V / Formula unit (Angstrom$^{3}$)')
plt.ylabel('P (GPa)')
plt.ticklabel_format(useOffset=False)
plt.savefig('P_vs_V.pdf', bbox_inches='tight')
# Plotting H vs P:
fig = plt.figure()
xp_C_I = np.linspace(PressureCI[0], PressureCI[-1], 100)
V_C_I_lin = np.linspace(VolumeCI[0], VolumeCI[-1], 100)
# Linspace for plotting the fitting curve:
V_C_I_lin = np.linspace(VolumeCI[0], VolumeCI[-1], 100)
p2, = plt.plot(xp_C_I, H(V_C_I_lin, *popt_C_I), color='grey', label='H fit Data' )
# Plotting the scattered points:
p1 = plt.scatter(pressures_per_F_unit_C_I, H_C_I, color='red', marker="^", label='Data', s=100)
fontP = FontProperties()
fontP.set_size('small')
plt.legend((p1, p2 ), ("Data", 'H fit Data' ), prop=fontP)
plt.xlabel('P / Formula unit (GPa)')
plt.ylabel('H / Formula unit (a.u.)')
plt.ticklabel_format(useOffset=False)
plt.savefig('H_vs_V.pdf', bbox_inches='tight')
plt.show()
您不能只绘制 H(V) 与一些完全不相关的压力 xp_C_I。
plt.plot(xp_C_I, H(V_C_I_lin, *popt_C_I), )
相反，您需要针对 P(V) 绘制 H(V)，这样压力和焓使用完全相同的 V 值：
plt.plot(P(V_C_I_lin, *popt_C_I), H(V_C_I_lin, *popt_C_I), )
我有 3 个函数：E（能量），P（压力）和 H（焓）。
给定数据 data.dat，包含变量 V（体积）和 E（能量）：
# Volume: V_C_I Energy: E_C_I
111.593876 -1.883070511360E+03
113.087568 -1.883074916825E+03
114.632273 -1.883078906679E+03
116.184030 -1.883082373429E+03
117.743646 -1.883085344696E+03
119.326853 -1.883087860954E+03
120.927806 -1.883089938181E+03
122.538335 -1.883091557526E+03
124.158641 -1.883092750745E+03
125.789192 -1.883093540824E+03
125.790261 -1.883093800747E+03
127.176327 -1.883094160364E+03
128.358654 -1.883094017730E+03
128.542807 -1.883094255789E+03
129.977279 -1.883094094751E+03
131.390610 -1.883093689121E+03
132.812287 -1.883093053342E+03
134.242765 -1.883092185844E+03
135.682211 -1.883091101112E+03
137.130792 -1.883089807766E+03
138.588565 -1.883088314435E+03
以下脚本执行：
1) E 的 curve_fit 与 V、
2) 使用def(P)函数计算P（压力），
3) 使用 def(H) 函数计算 H（焓）。 (H = E + PV).
4) 使用拟合曲线执行 3 个图：E 与 P、P 与 V 和 H 与 P .
绘制拟合曲线时，我使用了以下内容：
例如，对于 E 与 V 曲线：
plt.plot(V_C_I_lin, E(V_C_I_lin, *popt_C_I)
其中 V_C_I_lin 是数据点之间的体积的线性空间。
这看起来不错：
同样，对于 P 与 V 曲线，类似的方案：
plt.plot(V_C_I_lin, P(V_C_I_lin, *popt_C_I), color='grey', label='E fit Data' )
产生期望的结果：
脚本将每个数据点的压力保存在文件 ./E_V_P_H__C_I.dat 中。
然而，对于 H 与 P 曲线，类似的方案：
plt.plot(xp_C_I, H(V_C_I_lin, *popt_C_I), color='grey', label='H fit Data' )
其中 xp_C_I 是压力的线性空间，不会产生正确的拟合：
为什么会出现第三种情况？
代码：
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
import sys
import os
# Intial candidates for fit
E0_init = -941.510817926696
V0_init = 63.54960592453
B0_init = 76.3746233515232
B0_prime_init = 4.05340727164527
def E(V, E0, V0, B0, B0_prime):
return E0+ (2.293710449E+17)*(1E-21)*( (9.0/16.0)*(V0*B0) * ( (((V0/V)**(2.0/3.0)-1.0)**3.0)*B0_prime + ((V0/V)**(2.0/3.0)-1)**2 * (6.0-4.0*(V0/V)**(2.0/3.0)) ))
def P(V, E0, V0, B0, B0_prime):
f0=(3.0/2.0)*B0
f1=((V0/V)**(7.0/3.0))-((V0/V)**(5.0/3.0))
f2=((V0/V)**(2.0/3.0))-1
pressure= f0*f1*(1+(3.0/4.0)*(B0_prime-4)*f2)
return pressure
def H(V, E0, V0, B0, B0_prime):
return E(V, E0, V0, B0, B0_prime) + P(V, E0, V0, B0, B0_prime) * V
# Data (Red triangles):
V_not_p_f_unit_C_I, E_not_p_f_unit_C_I = np.loadtxt('data.dat', skiprows = 1).T
# A minor conversion of the data:
nFU_C_I = 2.0
nFU_C_II = 4.0
E_C_I = E_not_p_f_unit_C_I/nFU_C_I
V_C_I = V_not_p_f_unit_C_I/nFU_C_I
######
# Fitting and obtaining the parameters:
init_vals = [E0_init, V0_init, B0_init, B0_prime_init]
popt_C_I, pcov_C_I = curve_fit(E, V_C_I, E_C_I, p0=init_vals)
# Calculation of P:
pressures_per_F_unit_C_I = P(V_C_I, *popt_C_I)
# Calculation of H: H = E + PV
H_C_I = E_C_I + pressures_per_F_unit_C_I * V_C_I
# We save E, P and H into a file:
output_array_3 = np.vstack((E_C_I, V_C_I, pressures_per_F_unit_C_I, H_C_I)).T
np.savetxt('E_V_P_H__C_I.dat', output_array_3, header="Energy / FU (a.u.) \t Volume / FU (A^3) \t Pressure / F.U. (GPa) \t Enthalpy (a.u.)", fmt="%0.13f")
EnergyCI, VolumeCI, PressureCI, EnthalpyCI = np.loadtxt('./E_V_P_H__C_I.dat', skiprows = 1).T
# Plotting E vs V:
fig = plt.figure()
V_C_I_lin = np.linspace(VolumeCI[0], VolumeCI[-1], 100)
# Linspace for plotting the fitting curve:
V_C_I_lin = np.linspace(VolumeCI[0], VolumeCI[-1], 100)
p2, = plt.plot(V_C_I_lin, E(V_C_I_lin, *popt_C_I), color='grey', label='E fit Data' )
# Plotting the scattered points:
p1 = plt.scatter(VolumeCI, EnergyCI, color='red', marker="^", label='Data', s=100)
fontP = FontProperties()
fontP.set_size('small')
plt.legend((p1, p2 ), ("Data", 'E fit Data' ), prop=fontP)
plt.xlabel('V / Formula unit (Angstrom$^{3}$)')
plt.ylabel('E / Formula unit (a.u.)')
plt.ticklabel_format(useOffset=False)
plt.savefig('E_vs_V.pdf', bbox_inches='tight')
# Plotting P vs V:
fig = plt.figure()
# Linspace for plotting the fitting curve:
xp_C_I = np.linspace(PressureCI[-1], PressureCI[0], 100)
# Plotting the fitting curves:
p2, = plt.plot(V_C_I_lin, P(V_C_I_lin, *popt_C_I), color='grey', label='E fit Data' )
# Plotting the scattered points:
p1 = plt.scatter(VolumeCI, PressureCI, color='red', marker="^", label='Data', s=100)
fontP = FontProperties()
fontP.set_size('small')
plt.legend((p1, p2), ("Data", "E fit Data"), prop=fontP)
plt.xlabel('V / Formula unit (Angstrom$^{3}$)')
plt.ylabel('P (GPa)')
plt.ticklabel_format(useOffset=False)
plt.savefig('P_vs_V.pdf', bbox_inches='tight')
# Plotting H vs P:
fig = plt.figure()
xp_C_I = np.linspace(PressureCI[0], PressureCI[-1], 100)
V_C_I_lin = np.linspace(VolumeCI[0], VolumeCI[-1], 100)
# Linspace for plotting the fitting curve:
V_C_I_lin = np.linspace(VolumeCI[0], VolumeCI[-1], 100)
p2, = plt.plot(xp_C_I, H(V_C_I_lin, *popt_C_I), color='grey', label='H fit Data' )
# Plotting the scattered points:
p1 = plt.scatter(pressures_per_F_unit_C_I, H_C_I, color='red', marker="^", label='Data', s=100)
fontP = FontProperties()
fontP.set_size('small')
plt.legend((p1, p2 ), ("Data", 'H fit Data' ), prop=fontP)
plt.xlabel('P / Formula unit (GPa)')
plt.ylabel('H / Formula unit (a.u.)')
plt.ticklabel_format(useOffset=False)
plt.savefig('H_vs_V.pdf', bbox_inches='tight')
plt.show()
您不能只绘制 H(V) 与一些完全不相关的压力 xp_C_I。
plt.plot(xp_C_I, H(V_C_I_lin, *popt_C_I), )
相反，您需要针对 P(V) 绘制 H(V)，这样压力和焓使用完全相同的 V 值：
plt.plot(P(V_C_I_lin, *popt_C_I), H(V_C_I_lin, *popt_C_I), )