激光干涉振动扬声器系统模型曲线拟合失败求助

1. 初始参数数量级严重错误 #

你给出的初始参数p0 = (1.73, 0.94, 0.10*10**-3, 1.5*10**3)中，振膜质量m=0.94kg完全不符合物理常识——扬声器振膜质量通常在**0.1_10克（1e-41e-2 kg）**量级。这个错误会导致：

• 计算项ω²*(m/k)的数值异常巨大（比如ω=2π100≈628 rad/s时，ω²m/k≈628²*0.94/1500≈247）
• 理论模型的分母会变得极大，输出值远小于你的实验数据ampl
• 拟合算法无法找到合理的参数收敛方向，最终只能输出近似常数的直线

建议先修正初始参数的数量级，比如把m改成0.00094（0.94克），b调整到0.01（更符合实际的阻尼系数），示例初始参数：

p0 = (1.73, 0.00094, 0.01, 1.5*10**3)

2. 输入数据未转为numpy数组 #

你的omega和ampl都是Python列表，虽然curve_fit兼容列表，但转为numpy数组能避免潜在的类型运算问题，提升拟合稳定性：

# 替换原omega和ampl的定义
omega = np.array([x*2*np.pi for x in freq])
x0 = np.array([cycle*633/4 for cycle in cycles])
i0 = np.array([curr*np.sqrt(2) for curr in current])
ampl = x0 / i0

3. 增加参数边界约束 #

为了避免拟合出物理上不合理的参数（如质量、阻尼、刚度为负），可以给curve_fit添加边界条件：

# 设置所有参数大于0
param, covariance = optimize.curve_fit(
    form, omega, ampl, p0=p0,
    bounds=(0, [np.inf, np.inf, np.inf, np.inf])
)

4. 检查干涉振幅计算的正确性 #

激光干涉中，振膜移动距离x0与条纹变化数cycles的关系应为：

光程差变化 = 2x0 = cyclesλ → x0 = cycles*λ/2

你当前用的是x0 = cycle*633/4，相当于把振幅算成了实际值的一半。虽然这不会直接导致拟合出直线，但会让理论模型与实验数据的数值尺度不匹配，影响拟合效果，建议修正为：

x0 = np.array([cycle*633/2 for cycle in cycles])

修正后的完整拟合代码示例 #

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

freq = [100, 140, 180, 200, 220, 240, 260, 280, 290, 300, 310, 320, 340, 360, 400, 440, 480]
cycles = [18, 20.5, 23, 27, 26.5, 29, 26, 22.5, 28, 27.5, 27.5, 23, 14, 20, 12, 15.5, 13]
current = [1.20, 1.22, 1.24, 1.25, 1.05, 0.91, 0.61, 0.34, 0.32, 0.25, 0.25, 0.26, 0.45, 0.78, 0.78, 1.63, 1.90]

# 修正数据类型与振幅计算
omega = np.array([x*2*np.pi for x in freq])
x0 = np.array([cycle*633/2 for cycle in cycles])  # 修正干涉振幅公式
i0 = np.array([curr*np.sqrt(2) for curr in current])
ampl = x0 / i0

# 理论模型函数
def form(omega, Bl, m, b, k):
    return (Bl/k)/np.sqrt((1 - np.power(omega,2)*(m/k))**2 + np.power(omega,2)*(b/k)**2)

# 修正初始参数数量级
p0 = (1.73, 0.00094, 0.01, 1.5*10**3)

# 添加边界约束的拟合
param, covariance = optimize.curve_fit(form, omega, ampl, p0=p0, bounds=(0, [np.inf]*4))
Bl, m, b, k = param

# 绘图
plt.scatter(omega, ampl, label='实验数据')
plt.plot(omega, form(omega, *param), 'b-', label='拟合曲线')
plt.title("振幅 [x₀/i₀] × 角频率 [ω]")
plt.ylabel("振动振幅与电流振幅比 [x₀/i₀]")
plt.xlabel("角频率 [ω]")
plt.legend()
plt.show()

# 输出拟合参数
print(f"拟合参数：Bl={Bl:.4f}, m={m:.6f}kg, b={b:.4f}N·s/m, k={k:.2f}N/m")

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