R语言拟合普朗克黑体方程报错：初始参数梯度矩阵奇异

报错原因分析 #

Error in nlsModel(...) : singular gradient matrix at initial parameter estimates 本质是参数初始值不合理+量纲不匹配：

• 你的频率数据（2.27~21.33）应该是太赫兹（THz），但代码直接代入公式时未转换为国际单位（Hz），导致(h*frequency)/(k*t)的计算值极小，exp()结果接近1，分母趋近于0，数值稳定性极差，梯度计算失效。
• 初始温度设为2.5K完全不符合数据对应的黑体辐射特征，根据维恩位移定律估算，数据峰值对应的温度约90K左右。

修正后的代码 #

data <- data.frame(x=c(2.27, 2.72, 3.18, 3.63, 4.08, 4.54, 4.99, 5.45, 5.90, 6.35, 6.81, 7.26, 7.71, 8.17, 8.62, 9.08, 9.53, 9.98, 10.44, 10.89, 11.34, 11.80, 12.25, 12.71, 13.16, 13.61, 14.07, 14.52, 14.97, 15.43, 15.88, 16.34, 16.79, 17.24, 17.70, 18.15, 18.61, 19.06, 19.51, 19.97, 20.42, 20.87, 21.33),
                  y=c(200.723, 249.508, 293.024, 327.770, 354.081, 372.079, 381.493, 383.478, 378.901, 368.833, 354.063, 336.278, 316.076, 293.924, 271.432, 248.239, 225.940, 204.327, 183.262, 163.830, 145.750, 128.835, 113.568, 99.451, 87.036, 75.876, 65.766, 57.008, 49.223, 42.267, 36.352, 31.062, 26.580, 22.644, 19.255, 16.391, 13.811, 11.716, 9.921, 8.364, 7.087, 5.801, 4.523))
frequency <- data$x * 1e12  # 将THz转换为Hz，匹配公式单位
brightness <- data$y * 2.71057477e-3

# 定义普朗克黑体辐射函数
B <- function(frequency, t) {
  h <- 6.62607015e-34
  c <- 299792458
  k <- 1.380649e-23
  (2 * h * frequency^3 * c^-2) / (exp((h * frequency) / (k * t)) - 1)
}

library(stats)
# 根据维恩位移定律估算初始温度，设置合理初始值
fit <- nls(brightness ~ B(frequency, t), data = data, start = list(t = 90)) 

# 查看拟合结果
summary(fit)

library(ggplot2)
# 添加拟合值到数据集
data$fitted_brightness <- predict(fit)

# 绘制散点图+拟合曲线
ggplot(data, aes(x = x, y = brightness)) +
  geom_point(color = "steelblue", size = 2) +
  geom_line(aes(y = fitted_brightness), color = "firebrick", linewidth = 1) +
  labs(x = "Frequency (THz)", y = "Brightness", title = "Planck Blackbody Fit to Measured Data") +
  theme_minimal()

关键修改点 #

• 单位统一：将频率转换为Hz，确保公式中物理量的量纲匹配
• 初始值修正：利用维恩位移定律估算合理的初始温度，避免梯度奇异问题
• 绘图优化：添加清晰的轴标签和主题，提升图表可读性

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