从外推拟合曲线计算指数衰减速率常数及粒子寿命

我尝试从外推拟合曲线计算粒子寿命（寿命为衰减常数的倒数），但得到的衰减常数值不符合预期，相关绘图代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import make_smoothing_spline
from matplotlib.ticker import (MultipleLocator, FormatStrFormatter,AutoMinorLocator)
import math
Time=[0.1,0.5,1,5,10,50,100,250,500,1000,2500,5000]
Intensity=[2722.164194,2877.627742,2663.520645,2708.928125,2545.461613,2421.236129,1885.837742,1710.483871,1275.428387,776.0895806,192.4806452,26.35279]
error=[365.0612668,365.9075468,401.1165559,405.5068204,320.8171197,317.7694792,266.8600888,157.3938377,65.85350869,33.65902988,24.90421016,2.410679068]
Time2=[50,100,250,500,1000,2500,5000]
Intensity2=[2421.236129,1885.837742,1710.483871,1275.428387,776.0895806,192.4806452,26.35279]
error2=[317.7694792,266.8600888,157.3938377,65.85350869,33.65902988,24.90421016,2.410679068]
fig, ax = plt.subplots()
ax.scatter(Time, Intensity, label="Measured", c='r',s=10)

x2 = np.logspace(start=-1, stop=4, num=100)
lam=3
smooth_func = make_smoothing_spline(np.log(Time2), Intensity2, lam=lam)
y2 = smooth_func(np.log(x2))
ax.semilogx(x2, y2,label="Curve Fit Extrapolated")

ax.set_ylim([0, 5500])
ax.set_xlabel("Trapping time (ms)")
ax.legend(loc="best")
ax.set_ylabel("Average Intensity (counts/s)")
ax.yaxis.set_minor_locator(MultipleLocator(100))
ax.yaxis.tick_left()
ax.yaxis.set_ticks_position('both')
ax.xaxis.tick_bottom()
ax.xaxis.set_ticks_position('both')
ax.set_title("Average Intentsity Extrapolated vs Trapping Time after 30 Pulses",c='g')
plt.errorbar(Time, Intensity, yerr=error ,fmt="o",c='r')
ax.yaxis.set_minor_locator(MultipleLocator(100))
print(delta)
plt.grid()
plt.show()

问题根源 #

当前使用的make_smoothing_spline是平滑样条插值工具，仅能生成一条贴合数据的平滑曲线，不针对指数衰减的物理模型做参数拟合。粒子强度衰减遵循指数模型：I(t) = I₀ * exp(-t/τ)（τ为寿命，1/τ为衰减常数），用平滑样条无法直接提取准确的衰减参数，这是结果不符合预期的核心原因。

修正方案 #

改用非线性最小二乘法直接拟合指数衰减模型，同时引入测量误差权重，确保拟合结果符合物理规律：

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator)

# 原始数据
Time = [0.1, 0.5, 1, 5, 10, 50, 100, 250, 500, 1000, 2500, 5000]
Intensity = [2722.164194, 2877.627742, 2663.520645, 2708.928125, 2545.461613, 2421.236129, 1885.837742, 1710.483871, 1275.428387, 776.0895806, 192.4806452, 26.35279]
error = [365.0612668, 365.9075468, 401.1165559, 405.5068204, 320.8171197, 317.7694792, 266.8600888, 157.3938377, 65.85350869, 33.65902988, 24.90421016, 2.410679068]
# 用于拟合的长时数据（可替换为全部数据）
Time2 = [50, 100, 250, 500, 1000, 2500, 5000]
Intensity2 = [2421.236129, 1885.837742, 1710.483871, 1275.428387, 776.0895806, 192.4806452, 26.35279]
error2 = [317.7694792, 266.8600888, 157.3938377, 65.85350869, 33.65902988, 24.90421016, 2.410679068]

# 定义指数衰减模型
def exp_decay(t, I0, tau):
    return I0 * np.exp(-t / tau)

# 加权拟合指数模型，传入误差作为权重
popt, pcov = curve_fit(exp_decay, Time2, Intensity2, sigma=error2, p0=[2500, 1000])
I0_fit, tau_fit = popt
decay_constant = 1 / tau_fit

# 输出关键参数
print(f"拟合初始强度I0: {I0_fit:.2f} counts/s")
print(f"粒子寿命τ: {tau_fit:.2f} ms")
print(f"衰减常数: {decay_constant:.6f} ms⁻¹")

# 绘图
fig, ax = plt.subplots()
ax.scatter(Time, Intensity, label="Measured", c='r', s=10)
plt.errorbar(Time, Intensity, yerr=error, fmt="o", c='r')

# 生成拟合曲线数据
x_fit = np.logspace(start=-1, stop=4, num=100)
y_fit = exp_decay(x_fit, I0_fit, tau_fit)
ax.semilogx(x_fit, y_fit, label="Exponential Fit", c='b')

# 绘图样式设置
ax.set_ylim([0, 5500])
ax.set_xlabel("Trapping time (ms)")
ax.legend(loc="best")
ax.set_ylabel("Average Intensity (counts/s)")
ax.yaxis.set_minor_locator(MultipleLocator(100))
ax.yaxis.tick_left()
ax.yaxis.set_ticks_position('both')
ax.xaxis.tick_bottom()
ax.xaxis.set_ticks_position('both')
ax.set_title("Average Intensity vs Trapping Time after 30 Pulses", c='g')
plt.grid()
plt.show()

核心说明 #

1. 模型匹配：用exp_decay函数直接对应粒子衰减的物理规律，确保拟合参数有明确物理意义。
2. 加权拟合：通过sigma=error2传入测量误差，让拟合过程优先信任误差小的数据点，结果更可靠。
3. 初始值设置：p0=[2500, 1000]给出参数初始猜测，帮助拟合算法快速收敛到合理结果。
4. 参数提取：拟合得到的tau_fit就是粒子寿命，衰减常数为其倒数，直接满足你的计算需求。

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