如何用Python拟合Arrhenius方程求解A、B参数？

需要拟合阿伦尼乌斯方程 (k_{gas}=Ae^{(-B/(RT))}) 求解参数A和B（T为开尔文温度，A单位为mL/min），最终外推0℃至35℃的速率函数分析酵母活性。但多次尝试拟合均得到偏差极大的直线结果，相关原始数据与尝试代码如下：

原始数据绘图代码 #

import matplotlib.pyplot as plt
import numpy as np
import scipy as scipy

temperature = np.array([20,25,30,33,35,37,40,45])                  # 单位：℃
prod_rate = np.array([6.9,9.9,21.2,26.1,32.9,30.4,18.8,18.5])/30   # 单位：mL/min

plt.plot(temperature,prod_rate,'ok')
plt.ylabel('Production (mL/min)')
plt.xlabel('Temperature (C)')

首次拟合尝试代码 #

import matplotlib.pyplot as plt
import numpy as np
import scipy as scipy

def arrhenius(T, A, B):
    R = 8.314  # Gas constant, J/(mol·K)
    return A * np.exp(-B / (R * T))

temperature = np.array([20,25,30,33,35,37,40,45])                  # 单位：℃
prod_rate = np.array([6.9,9.9,21.2,26.1,32.9,30.4,18.8,18.5])/30   # 单位：mL/min

plt.plot(temperature, prod_rate, 'ok')
plt.ylabel('Production (mL/min)')
plt.xlabel('Temperature (C)')

popt, pcov = curve_fit(arrhenius, temperature, prod_rate)
print(popt)

A_optimal, B_optimal = popt
fitted_curve = arrhenius(temperature, A_optimal, B_optimal)

plt.plot(temperature, prod_rate, 'ok', label='Data')
plt.plot(temperature, fitted_curve, '-r', label='Fitted Curve')
plt.ylabel('Production (mL/min)')
plt.xlabel('Temperature (°C)')
plt.legend()
plt.show()

print("Optimal value for A:", A_optimal, "mL/min")
print("Optimal value for B:", B_optimal)

转开尔文后的拟合代码 #

import matplotlib.pyplot as plt
import numpy as np
import scipy as scipy
from scipy.optimize import curve_fit

def arrhenius(T, A, B):
    R = 8.314  # Gas constant, J/(mol·K)
    return A * np.exp(-B / (R * T))

temperature = np.array([20,25,30,33,35,37,40,45])  
temperature_K = temperature + 273.15                # 转换为开尔文
prod_rate = np.array([6.9,9.9,21.2,26.1,32.9,30.4,18.8,18.5])/30   # 单位：mL/min

plt.plot(temperature_K, prod_rate, 'ok')
plt.ylabel('Production (mL/min)')
plt.xlabel('Temperature (C)')

popt, pcov = curve_fit(arrhenius, temperature_K, prod_rate)
print(popt)

A_optimal, B_optimal = popt
fitted_curve = arrhenius(temperature_K, A_optimal, B_optimal)

plt.plot(temperature_K, prod_rate, 'ok', label='Data')
plt.plot(temperature_K, fitted_curve, '-r', label='Fitted Curve')
plt.ylabel('Production (mL/min)')
plt.xlabel('Temperature (°C)')
plt.legend()
plt.show()

print("Optimal value for A:", A_optimal, "mL/min")
print("Optimal value for B:", B_optimal)

问题原因分析 #

1. 模型适用范围错误：标准阿伦尼乌斯方程描述温度升高反应速率单调上升的过程，但酵母活性是钟型曲线——温度超过35℃后酶变性导致活性下降，这部分数据完全不符合阿伦尼乌斯模型，用全部数据拟合会干扰优化器找到正确的指数趋势。
2. 缺少初始参数猜测：curve_fit默认初始参数为[1,1]，与真实解差距极大，容易陷入局部最优，最终拟合出无意义的直线。

修正方案与代码 #

核心修正点 #

• 仅使用35℃及以下的有效数据（符合阿伦尼乌斯模型的温度区间）
• 为curve_fit提供合理的初始参数猜测，引导优化器找到正确解
• 生成0℃至35℃的外推曲线，满足分析需求

修正后代码 #

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

def arrhenius(T, A, B):
    R = 8.314  # 气体常数，J/(mol·K)
    return A * np.exp(-B / (R * T))

# 筛选35℃及以下的有效数据
temp_c_valid = np.array([20,25,30,33,35])
temp_k_valid = temp_c_valid + 273.15
prod_rate_valid = np.array([6.9,9.9,21.2,26.1,32.9])/30  # mL/min

# 设置初始参数猜测：A略大于最大产率，B取酶促反应典型活化能范围
p0 = [2, 10000]
popt, pcov = curve_fit(arrhenius, temp_k_valid, prod_rate_valid, p0=p0)

A_opt, B_opt = popt
print(f"最优A值: {A_opt:.2f} mL/min")
print(f"最优B值: {B_opt:.2f} J/mol")

# 生成0℃到35℃的外推区间
extrapolate_c = np.linspace(0, 35, 100)
extrapolate_k = extrapolate_c + 273.15
extrapolate_rate = arrhenius(extrapolate_k, A_opt, B_opt)

# 绘图展示
plt.figure(figsize=(8,5))
# 全部实验数据点
all_temp_c = np.array([20,25,30,33,35,37,40,45])
all_rate = np.array([6.9,9.9,21.2,26.1,32.9,30.4,18.8,18.5])/30
plt.scatter(all_temp_c, all_rate, color='k', label='全部实验数据')
# 用于拟合的有效数据点
plt.scatter(temp_c_valid, prod_rate_valid, color='blue', s=100, label='拟合用有效数据')
# 拟合曲线与外推结果
plt.plot(extrapolate_c, extrapolate_rate, 'r-', label='阿伦尼乌斯拟合+外推')

plt.xlabel('温度 (℃)')
plt.ylabel('产气速率 (mL/min)')
plt.legend()
plt.grid(True)
plt.show()

结果说明 #

• 拟合得到的A、B符合酶促反应的典型参数范围，曲线能准确贴合35℃以下的活性上升趋势
• 外推的0℃至35℃速率函数可直接用于分析低温到最适温度区间的酵母活性
• 37℃及以上数据因酶变性偏离模型，禁止纳入阿伦尼乌斯拟合

内容的提问来源于stack exchange，提问作者bob micheal