为何curve_fit函数无法正确拟合单摆衰减数据？

我尝试通过拟合单摆实验数据来获取衰减时间常数，但使用scipy.optimize.curve_fit拟合全量数据时结果完全不对，截断前1500组数据后拟合效果却很理想。

使用的数据与代码 #

实验数据包含时间t、x方向位移x、y方向位移y序列，代码如下：

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

filepath1 = '1.txt'

t1, x1, y1 = np.genfromtxt(filepath1, dtype=float, delimiter=' ', skip_header=2, unpack=True)

rad1 = np.arctan(x1/y1)

# 取消注释后使用截断数据拟合
# t1, rad1 = t1[:1500], rad1[:1500]

error = 0.05*np.ones(len(t1))

def angle(t, theta_0, tau, T, phi):
    return theta_0 * np.exp(-t/tau) * np.cos(((2*np.pi)/T) * t + phi)

popt1, pcov1 = curve_fit(angle, t1, rad1, sigma=error, p0=(0.45,50,1.3,0.05), absolute_sigma=True, method='lm')

theta_0, tau, T, phi = popt1

print(tau,T)
plt.plot(t1,rad1)
plt.plot(t1,angle(t1,theta_0, tau, T, phi))

拟合效果对比 #

• 全量数据拟合结果：
  
• 截断数据（前1500组）拟合结果：
  

问题原因与解决思路 #

1. 角度计算错误 #

np.arctan(x1/y1)仅能输出[-π/2, π/2]范围内的角度，单摆运动中当y坐标过零点时，会出现角度跳变，导致数据不连续。改用np.arctan2(x1, y1)可计算[-π, π]全范围的连续角度，这是最可能的核心问题：

rad1 = np.arctan2(x1, y1)

2. 初始值与算法局限性 #

LM算法对初始参数p0敏感，全量数据包含更多衰减周期，容易陷入局部最优。解决方法：

• 先用截断数据拟合得到的参数作为全量拟合的初始值（比如从截断拟合结果中提取tau、T等参数代入p0）
• 改用鲁棒性更强的method='trf'算法，替代默认的LM算法

3. 噪声权重不合理 #

固定值error=0.05不符合实际数据的噪声分布，后期数据噪声通常更大。可通过滚动标准差估算动态噪声：

window_size = 50
# 计算相邻数据的波动作为噪声估计
diff = np.abs(np.diff(rad1))
# 滚动平均平滑噪声
error = np.convolve(diff, np.ones(window_size)/window_size, mode='same')
# 补全第一个元素
error = np.insert(error, 0, error[0])

4. 模型适配性问题 #

全量数据后期可能因摩擦非线性、传感器漂移等原因，不再符合小角度阻尼单摆的指数衰减余弦模型。如果修正角度和参数后仍无法拟合全量数据，可考虑：

• 仅使用前期符合模型的数据进行拟合
• 引入更复杂的阻尼模型（比如二次阻尼项）

内容的提问来源于stack exchange，提问作者A.S