使用Python进行指数曲线拟合

Question

我想知道是否有人知道如何使用Python找到阿伦尼乌斯方程（𝑘𝑔𝑎𝑠=𝐴𝑒^（−𝐵𝑅𝑇））中A和B的值。在此方程中，𝑇 的单位为开尔文，A 的单位为 mL/min。我们希望获得 0°C 到 35°C 范围内温度函数的速率描述，假设我们可以推断。

基本上，我们试图找出酵母随温度的活性，但首先我们必须找到 A 和 B。我尝试了几种曲线拟合技术，但都失败了。我的身材是一条直线，而且我的数字真的很差。

这是图像中给出的代码/信息：

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

temperature = np.array([20,25,30,33,35,37,40,45])                  # degrees C
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])                  # degrees C
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                # degrees C
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)

Answer 1

假设你的实验做错了。选择一个子数组来拟合，并以对数尺度进行拟合，因为线性尺度很大。请使用开尔文，而不是摄氏度。这里我假设 A 的单位是任意的，因为生产率的单位尚未修正。

import matplotlib.pyplot as plt
import numpy as np

# Gas constant, J/(mol·K); and 0C in K
from scipy.constants import R, zero_Celsius
from scipy.optimize import curve_fit

temperature_C = np.array((20, 25, 30, 33, 35, 37, 40, 45))
temperature_K = temperature_C + zero_Celsius
prod_rate_mL_min = np.array((6.9,9.9,21.2,26.1,32.9,30.4,18.8,18.5))/30   # mL/min


def arrhenius(T: np.ndarray, A: float, B: float) -> np.ndarray:
    return A * np.exp(-B/R/T)


def arrhenius_log(T: np.ndarray, logA: float, B: float) -> np.ndarray:
    return logA - B/R/T


# The data are garbage; cut off the last three points during fit
valid_end = -3

# Estimate calculation, partially based on experimental data
B0 = 80_000
logA0 = np.log(
    prod_rate_mL_min[valid_end - 1]
) + B0/R/temperature_K[valid_end - 1]
p0 = logA0, B0

(logA, B), _ = curve_fit(
    f=arrhenius_log, p0=p0,
    xdata=temperature_K[:valid_end],
    ydata=np.log(prod_rate_mL_min[:valid_end]),
)

# Keep A in log space for accuracy's sake
a_pow = np.floor(logA / np.log(10))
a_mantissa = np.exp(logA - a_pow*np.log(10))
print(f'A = {a_mantissa:.6f}e{a_pow:.0f}')
print(f'B = {B:.3f}')

fig, ax = plt.subplots()
ax.scatter(
    temperature_C[:valid_end], prod_rate_mL_min[:valid_end],
    label='experiment (OK)',
)
ax.scatter(
    temperature_C[valid_end:], prod_rate_mL_min[valid_end:],
    label='experiment (bad)',
)
T_hires = np.linspace(
    start=temperature_K.min(),
    stop=temperature_K.max(),
    num=201,
)
ax.plot(
    T_hires - zero_Celsius,
    np.exp(arrhenius_log(T_hires, *p0)),
    label='guess',
)
ax.plot(
    T_hires - zero_Celsius,
    arrhenius(
        T=T_hires, A=a_mantissa * 10**a_pow, B=B,
    ),
    label='fit',
)
ax.set_xlabel('T (C)')
ax.set_ylabel('Rate (mL/min)')
ax.legend()
plt.show()

A = 6.119210e13
B = 81109.535

线性对数空间拟合

请注意，在对数空间中，表达式是线性的，因此您可以简化为计算上更简单的线性拟合，不需要估计：

import numpy as np

# Gas constant, J/(mol·K); and 0C in K
from scipy.constants import R, zero_Celsius

temperature_C = np.array((20, 25, 30, 33, 35, 37, 40, 45))
temperature_K = temperature_C + zero_Celsius
prod_rate_mL_min = np.array((6.9,9.9,21.2,26.1,32.9,30.4,18.8,18.5))/30   # mL/min

# The data are garbage; cut off the last three points during fit
valid_end = -3

(logA, B), *_ = np.linalg.lstsq(
    a=np.stack((
        np.ones(temperature_K.size + valid_end),
        -1/R/temperature_K[:valid_end],
    ), axis=1),
    b=np.log(prod_rate_mL_min[:valid_end]),
    rcond=None,
)

# Keep A in log space for accuracy's sake
a_pow = np.floor(logA / np.log(10))
a_mantissa = np.exp(logA - a_pow*np.log(10))
print(f'A = {a_mantissa:.6f}e{a_pow:.0f}')
print(f'B = {B:.3f}')