#作者：王日睿
#中国科学技术大学生命科学学院
#2021.11.20
#物理化学实验：实验04 液体饱和蒸气压的测定
# Jupyter lab
import numpy as np
import matplotlib.pyplot as plt
import math
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
# 线性拟合函数
from scipy import optimize
def linear_fit(x, y):
def func(x, k, b):
return k * x + b
k, b = optimize.curve_fit(func, x, y)[0]
return k, b
#求R^2函数
def get_R_square(x, y):
y_mean = np.mean(y)
k, b = linear_fit(x, y)
y_predict = [k * e + b for e in x]
TSS = [(x-y_mean)**2 for x in y]
RSS = [(y_predict[i] - y[i])**2 for i in range(len(y))]
return 1 - sum(RSS) / sum(TSS)
# 原始数据：
# 工作曲线：
U_left = [-173.0, -151.0, -137.6, -119.0, -97.0, -78.3, -58.5, -39.7, -19.8, -0.6]
U_right = [181.0, 163.0, 143.5, 125.0, 102.2, 82.9, 62.2, 42.8, 22.6, 1.2]
U_value = []
for i in range(len(U_left)):
U_value.append((U_left[i] - U_right[i]))
#print(U_value)
#[-354.0, -314.0, -281.1, -244.0, -199.2, -161.2, -120.7, -82.5, -42.4, -1.8]
Digit_value = [-350.5, -316.6, -278.7, -242.0, -197.0, -159.5, -119.6, -82.0, -41.5, 0]
work_k, work_b = linear_fit(Digit_value, U_value)
#print("工作曲线为y=%4.3fx+(%4.3f)"%(work_k, work_b))
#工作曲线为y=1.000x+(-1.268)
work_R_square = get_R_square(Digit_value, U_value)
print(work_R_square)
#0.9993061474052902
工作曲线为y=1.000x-1.268
拟合的R^2=0.9993061474052902
X = np.arange(-360, 10, 0.01)
Y = work_k * X + work_b
plt.figure(1)
plt.figure(figsize=(8,6), dpi=400)
plt.plot(X, Y, "b-")
plt.scatter(Digit_value, U_value, color="r")
plt.title("工作曲线")
plt.xlabel("数字负压仪读数/mmHg")
plt.ylabel("U型水银压力计读数/mmHg")
plt.savefig("work_curve")
plt.show()
# 三次测量的大气压平均值
P_air = 101.51
# 每一组的原始数据压力（digit_value）
digit_1 = [0, -33.4, -73.1, -115.5, -152.7, -195.5, -230.2, -269.7, -309.8, -350.9, -390.7]
digit_2 = [-390.6, -354.3, -313.3, -274.1, -236.6, -197.3, -159.8, -121.3, -80.8, -42.5, 0]
digit_3 = [0, -42.4, -74.8, -115.3, -155.1, -196.2, -235.1, -272.0, -314.3, -352.2, -391.0]
# 每一组的沸点（BT）
BT_1 = [80.639, 79.247, 77.288, 75.232, 73.301, 70.992, 68.993, 66.621, 64.015, 61.451, 58.190]
BT_2 = [58.177, 60.977, 63.811, 66.350, 68.630, 70.901, 72.954, 74.962, 76.968, 78.765, 80.649]
BT_3 = [80.647, 78.749, 77.226, 75.175, 73.151, 70.904, 68.698, 66.438, 63.715, 61.077, 58.121]
# 转为开尔文温度
BT_1 = [t+273.15 for t in BT_1]
BT_2 = [t+273.15 for t in BT_2]
BT_3 = [t+273.15 for t in BT_3]
# 每一组U型管压力（工作曲线得到 U=D-1.268）
# 试一试循环命令变量吧！
# 新技巧get!
for i in range(1, 4):
exec('U_' + str(i) + '=' + '[' + 'p-1.268 for p in digit_' + str(i)+']')
# 液体的饱和蒸气压（SVP）为：大气压-U型管的压差
for i in range(1, 4):
exec('SVP_' + str(i) + '=' + '[' + 'P_air + 101.3 / 760 * u for u in U_' + str(i)+']')
# ln(SVP)
for i in range(1, 4):
exec('ln_SVP_' + str(i) + '=' + '[' + 'math.log(p) for p in SVP_' + str(i)+']')
# 1/T
for i in range(1, 4):
exec('Inv_T_' + str(i) + '=' + '[' + '1/t for t in BT_' + str(i)+']')
# 开始拟合 k1, b1 = linear_fit(Inv_T_1, ln_SVP_1) # (-3776.6207672189535, 15.293833511600848) k2, b2 = linear_fit(Inv_T_2, ln_SVP_2) # (-3761.366318227884, 15.249891300848459) k3, b3 = linear_fit(Inv_T_3, ln_SVP_3) # (-3759.494812223895, 15.24566374088445)
# pearson系数 from scipy.stats import pearsonr Pearson_R = [pearsonr(dot_x[i], dot_y[i]) for i in range(3)] Pearson_R = [list(a) for a in Pearson_R] Pearson_R[0][0]
X = np.arange(0.002826, 0.003030, 1e-6)
dot_x = [Inv_T_1] + [Inv_T_2] + [Inv_T_3] # 早知道搞个二维列表，不用循环命名这么烦
dot_y = [ln_SVP_1] + [ln_SVP_2] + [ln_SVP_3]
k_list = [k1, k2, k3]
b_list = [b1, b2, b3]
R_list = [get_R_square(dot_x[i], dot_y[i]) for i in range(3)]
for i in range(3):
plt.figure(i+1)
Y = k_list[i] * X + b_list[i]
plt.figure(figsize=(6,4), dpi=400)
plt.plot(X, Y, "b-", label=r"$lnP=%5.4ffrac{1}{T}+%4.3f$"%(k_list[i], b_list[i]))
plt.scatter(dot_x[i], dot_y[i], color="r")
plt.xlabel(r"$frac{1}{T}$")
plt.ylabel("lnP")
plt.text(0.003025, 4.4, s=r"$R^2=%3.6f$"%R_list[i], fontsize=15,
verticalalignment="top", horizontalalignment="right")
plt.text(0.003025, 4.3, s=r"Pearson of R=%3.6f"%Pearson_R[i][0], fontsize=12,
verticalalignment="top", horizontalalignment="right")
plt.title(r"第%d组$lnP-frac{1}{T}$拟合曲线"%(i+1))
plt.grid(True)
plt.legend(loc="upper right")
plt.tight_layout()
plt.savefig("output_%d"%(i+1))
plt.show()
# 一些计算： # lnP = A - (delta H)/R * (1/T) R = 8.314 delta_H = [- R * K for K in k_list] delta_H # [31398.825058658378, 31271.99956974663, 31256.439868829464]
ave_H = sum(delta_H)/3 ave_H # 31309.088165744823
# 求沸点：
import sympy as sym
P = 101.325
boiling_t = []
T = sym.symbols("T")
for i in range(3):
t = sym.solve(math.log(P) - (k_list[i] / (T+273.15) + b_list[i]), T)
boiling_t = boiling_t + t
boiling_t
# [80.6152238548020, 80.6425742232915, 80.6072100573596]
ave_bt = sum(boiling_t)/3 ave_bt # 80.6216693784843
# 资料：环己烷 # 沸点 : 80.72 °C #相对误差为： delta = (80.72 - 80.62) / 80.72 * 100 delta # 0.12% 的误差
