扩展绘图
Extending plotting
我有以下代码用于使用实验数据进行理论分析。我正在尝试拟合曲线并推断 S_u0 值以确定初始温度。
# Determination of laminar burning velocities from experimental data using constant volume pressure rise method.
import numpy as np
import openpyxl as xl
import matplotlib.pylab as plt
from lmfit import Model
plt.rc('font',family='Book Antiqua',size=12)
p_i=float(2.0) # Initial pressure
c=float(2.0265) # Constant
V_b=float(1.94e-3) # Volume of combustion bomb
k_u=float(1.3) # Specific heat ratio of unburnt mass
k_b=float(1.27) # Specific heat ratio of burnt mass
# Extracting experimental data from excel.
wb = xl.load_workbook("G:/Laminar Paper/Data/T0444CH1.xlsx")
ws = wb.active
# Calculation of initial voltage from experimental values.
num_cells=list(ws['B2':'B1500'])
V_initial=[]
for row in num_cells:
for cell in row:
V_initial.append(cell.value)
V_av=np.absolute(sum(V_initial)/len(V_initial)) # Average voltage
# Determining the final voltage from the experimental values.
num_cells=list(ws['B3':'B13252'])
V_final=[] # Experimental testing voltages
for row in num_cells:
for cell in row:
V_final.append(cell.value)
p=(V_final/np.floor(c))+(p_i-(V_av/np.floor(c))) # Combustion pressure
# The time elapse.
num_cells=list(ws['A3':'A13252'])
time=[]
for row in num_cells:
for cell in row:
time.append(cell.value)
# Curve fitting.
Func_1=np.polyfit(time, p, 15)
Func_2=np.poly1d(Func_1)
time_new=np.linspace(time[0], time[-1], 50) # New time values
p_new=Func_2(time_new) # New pressure values/filtered pressure
dp_dt=np.gradient(p_new,time_new) # Derivative of pressure wrt time
# Linear approximation method.
x_linear=(p_new-p_i)/(max(p_new)-p_i) # Burnt mass fraction
R_b=(pow(((V_b*3)/(4*np.pi)), 1.0/3))
S_u1=((1-(1-x_linear)*(p_i/p_new)**1/k_u)**-2/3)
S_u2=S_u1*((p_i/p_new)*1/k_u)
S_u3=S_u2*(1/(max(p_new)-p_i))*dp_dt
S_u=(R_b/3)*S_u3 # Linear approximation laminar flame speed
# Two-zone approach.
func_p=(k_b-1)/(k_u-1)+(((k_u-k_b)/(k_u-1))*(p_new/p_i)**((k_u-1)/k_u))
x_2zone=(p_new-p_i*func_p)/(max(p)-p_i*func_p) # Burnt mass fraction
S_uz1=((1-(1-x_2zone)*(p_i/p_new)**1/k_u)**-2/3)
S_uz2=S_uz1*((p_i/p_new)*1/k_u)
S_uz3=S_uz2*(1/(max(p_new)-p_i))*dp_dt
S_uz=(R_b/3)*S_uz3 # Two-zone laminar flame speed
# Model function.
def mod_m(n,su0=1,alpha=1):
return su0*(n)**alpha
# S_u and pressure ratio array data for fitting and calculating S_u0 and T_0.
n=np.array([1.78759326,1.91173221,2.05673713,2.22911249,2.4360555,2.68433232,2.97885016])
m=np.array([0.20791161,0.20332239,0.20206572,0.20260075,0.20292399,0.20115307,0.19601662])
# Fitting model.
model=Model(mod_m)
# Making a set of parameters (for 'su0' and 'alpha'):
params=model.make_params(s_u0=10)
# Setting min/max bounds on parameters:
params['alpha'].min=0.0
params['su0'].min=0.0
params['su0'].max=1e2
# Run the fit with Model.fit(Data_Array, Parameters, independent vars).
result=model.fit(m,params,n=n)
# Plotting.
plt.plot(n,result.best_fit,'b*-', label='Fit')
plt.plot(p_new/p_i,S_u,'k--',linewidth='2',label='Linear')
plt.plot(p_new/p_i,S_uz,'r-',label='Two-zone')
plt.xlabel('Pressure ratio $p/p_i$')
plt.ylabel('Laminar flame speed $(S_u)$ [m/s]')
plt.legend(loc='lower right')
plt.savefig('Laminar.png')
plt.show()
# Printing report with results and fitting statistics.
print(result.fit_report())
绘图后我得到了图 1 但我想得到像图 2 的图。请问我怎样才能得到像图 2 的图。提前致谢。
图 1
图2
如果我理解这个问题，您希望绘制超出拟合范围的模型预测值。
n 和 m 数组的范围有限（至少与其他数据相比）
result = model.fit(m, params, n=n)
plt.plot(n, result.best_fit, 'b*-', label='Fit')
将仅在您用于拟合的有限范围内绘制。
您可以使用一组参数（您可能需要最适合的参数）和自变量的任何值来评估模型 n
result.eval():
new_m = result.eval(result.params, n=new_n)
看来你想要的可能是这样的
new_n = np.linspace(0, 5, 51)[1:] # to avoid 0, where your model fails
plt.plot(new_n, result.eval(result.params, n=new_n), label='Predicted')
我有以下代码用于使用实验数据进行理论分析。我正在尝试拟合曲线并推断 S_u0 值以确定初始温度。
# Determination of laminar burning velocities from experimental data using constant volume pressure rise method.
import numpy as np
import openpyxl as xl
import matplotlib.pylab as plt
from lmfit import Model
plt.rc('font',family='Book Antiqua',size=12)
p_i=float(2.0) # Initial pressure
c=float(2.0265) # Constant
V_b=float(1.94e-3) # Volume of combustion bomb
k_u=float(1.3) # Specific heat ratio of unburnt mass
k_b=float(1.27) # Specific heat ratio of burnt mass
# Extracting experimental data from excel.
wb = xl.load_workbook("G:/Laminar Paper/Data/T0444CH1.xlsx")
ws = wb.active
# Calculation of initial voltage from experimental values.
num_cells=list(ws['B2':'B1500'])
V_initial=[]
for row in num_cells:
for cell in row:
V_initial.append(cell.value)
V_av=np.absolute(sum(V_initial)/len(V_initial)) # Average voltage
# Determining the final voltage from the experimental values.
num_cells=list(ws['B3':'B13252'])
V_final=[] # Experimental testing voltages
for row in num_cells:
for cell in row:
V_final.append(cell.value)
p=(V_final/np.floor(c))+(p_i-(V_av/np.floor(c))) # Combustion pressure
# The time elapse.
num_cells=list(ws['A3':'A13252'])
time=[]
for row in num_cells:
for cell in row:
time.append(cell.value)
# Curve fitting.
Func_1=np.polyfit(time, p, 15)
Func_2=np.poly1d(Func_1)
time_new=np.linspace(time[0], time[-1], 50) # New time values
p_new=Func_2(time_new) # New pressure values/filtered pressure
dp_dt=np.gradient(p_new,time_new) # Derivative of pressure wrt time
# Linear approximation method.
x_linear=(p_new-p_i)/(max(p_new)-p_i) # Burnt mass fraction
R_b=(pow(((V_b*3)/(4*np.pi)), 1.0/3))
S_u1=((1-(1-x_linear)*(p_i/p_new)**1/k_u)**-2/3)
S_u2=S_u1*((p_i/p_new)*1/k_u)
S_u3=S_u2*(1/(max(p_new)-p_i))*dp_dt
S_u=(R_b/3)*S_u3 # Linear approximation laminar flame speed
# Two-zone approach.
func_p=(k_b-1)/(k_u-1)+(((k_u-k_b)/(k_u-1))*(p_new/p_i)**((k_u-1)/k_u))
x_2zone=(p_new-p_i*func_p)/(max(p)-p_i*func_p) # Burnt mass fraction
S_uz1=((1-(1-x_2zone)*(p_i/p_new)**1/k_u)**-2/3)
S_uz2=S_uz1*((p_i/p_new)*1/k_u)
S_uz3=S_uz2*(1/(max(p_new)-p_i))*dp_dt
S_uz=(R_b/3)*S_uz3 # Two-zone laminar flame speed
# Model function.
def mod_m(n,su0=1,alpha=1):
return su0*(n)**alpha
# S_u and pressure ratio array data for fitting and calculating S_u0 and T_0.
n=np.array([1.78759326,1.91173221,2.05673713,2.22911249,2.4360555,2.68433232,2.97885016])
m=np.array([0.20791161,0.20332239,0.20206572,0.20260075,0.20292399,0.20115307,0.19601662])
# Fitting model.
model=Model(mod_m)
# Making a set of parameters (for 'su0' and 'alpha'):
params=model.make_params(s_u0=10)
# Setting min/max bounds on parameters:
params['alpha'].min=0.0
params['su0'].min=0.0
params['su0'].max=1e2
# Run the fit with Model.fit(Data_Array, Parameters, independent vars).
result=model.fit(m,params,n=n)
# Plotting.
plt.plot(n,result.best_fit,'b*-', label='Fit')
plt.plot(p_new/p_i,S_u,'k--',linewidth='2',label='Linear')
plt.plot(p_new/p_i,S_uz,'r-',label='Two-zone')
plt.xlabel('Pressure ratio $p/p_i$')
plt.ylabel('Laminar flame speed $(S_u)$ [m/s]')
plt.legend(loc='lower right')
plt.savefig('Laminar.png')
plt.show()
# Printing report with results and fitting statistics.
print(result.fit_report())
绘图后我得到了图 1 但我想得到像图 2 的图。请问我怎样才能得到像图 2 的图。提前致谢。
图2
如果我理解这个问题，您希望绘制超出拟合范围的模型预测值。
n 和 m 数组的范围有限（至少与其他数据相比）
result = model.fit(m, params, n=n)
plt.plot(n, result.best_fit, 'b*-', label='Fit')
将仅在您用于拟合的有限范围内绘制。
您可以使用一组参数（您可能需要最适合的参数）和自变量的任何值来评估模型 n
result.eval():
new_m = result.eval(result.params, n=new_n)
看来你想要的可能是这样的
new_n = np.linspace(0, 5, 51)[1:] # to avoid 0, where your model fails
plt.plot(new_n, result.eval(result.params, n=new_n), label='Predicted')