使用rootSolve求解含σ²₂和μ₂的非线性方程组遇异常问题

方程组说明 #

需要求解的非线性方程组为：

1. $\ln\left( \frac{\exp(2\mu_1+\sigma_1^{2)(\exp(\sigma_1}2)-1) + \exp(2\mu_2+\sigma_2^{2)(\exp(\sigma_2}2)-1)}{(\exp(\mu_1+\sigma_1^2/2) + \exp(\mu_2+\sigma_2^2/2))2} + 1 \right) = \sigma_Z^2$
2. $\ln\left( \exp(\mu_1+\sigma_1^2/2) + \exp(\mu_2+\sigma_2^2/2) \right) = \mu_Z + \frac{\sigma_Z^2}{2}$

已知参数：

• $\mu_1 = -2.931476$
• $\sigma_1^2 = (1.084857)^2$
• $\mu_Z = -6.717234$
• $\sigma_Z^2 = 0.104320$

用户提供的R代码 #

library(rootSolve)
mu_1 <- -2.931476
sigmasq_1 <- (1.084857)^2
mu_Z <- -6.717234
sigmasq_Z <- 0.104320

model <- function(x) c(F1 = log(((exp(2*mu_1+sigmasq_1)*(exp(sigmasq_1)-1))+(exp(2*x[1]+x[2])*(exp(x[2])-1)))/
                                 ((exp(mu_1+sigmasq_1/2)+exp(x[1]+x[2]/2))^2)+1)-sigmasq_Z, 
                       F2 = log(exp(mu_1+sigmasq_1/2)+exp(x[1]+x[2]/2))-sigmasq_Z/2-mu_Z)
(ss <- multiroot(f = model, start = c(1.5, 0)))

运行后警告信息 #

Warning messages:
1: In stode(y, times, func, parms = parms, ...) :
error during factorisation of matrix (dgefa); singular matrix
2: In stode(y, times, func, parms = parms, ...) : steady-state not reached

问题排查与解决方法 #

1. 核心问题分析 #

原代码出现奇异矩阵警告，原因是二维求解时雅可比矩阵线性相关，且初始值选择不当（$\sigma_2^2=0$会导致部分项为0，破坏矩阵可逆性）。同时直接计算高次指数容易出现数值溢出/下溢问题。

2. 降维简化求解 #

观察第二个方程，可先解出中间变量$S = \exp(\mu_1+\sigma_1^2/2) + \exp(\mu_2+\sigma_2^2/2)$，由第二个方程直接得$S = \exp(\mu_Z + \sigma_Z^{2/2)$。据此可将$\mu_2$用$\sigma_2}2$表示，把二维问题降为一维，避免雅可比矩阵奇异问题。

3. 修改后的实现代码 #

library(rootSolve)

# 已知参数
mu_1 <- -2.931476
sigmasq_1 <- (1.084857)^2
mu_Z <- -6.717234
sigmasq_Z <- 0.104320

# 预计算已知项，减少重复计算并规避数值问题
term1_exp <- exp(mu_1 + sigmasq_1/2)
term1_var <- exp(2*mu_1 + sigmasq_1)*(exp(sigmasq_1)-1)
S <- exp(mu_Z + sigmasq_Z/2)
term2_exp <- S - term1_exp

# 检查term2_exp合法性，确保对数有意义
if(term2_exp <= 0) stop("term2_exp非正，无解或输入参数有误")

# 构建一维方程，未知量为σ2²
model_1d <- function(sigmasq_2) {
  mu_2 <- log(term2_exp) - sigmasq_2/2
  term2_var <- exp(2*mu_2 + sigmasq_2)*(exp(sigmasq_2)-1)
  numerator <- term1_var + term2_var
  left_side <- log(numerator / S^2 + 1)
  left_side - sigmasq_Z
}

# 求解一维方程，选择合理初始区间
result_sigmasq2 <- uniroot(model_1d, interval = c(0.01, 5))$root
result_mu2 <- log(term2_exp) - result_sigmasq2/2

# 输出结果
cat("求解结果：\n")
cat("μ2 =", result_mu2, "\n")
cat("σ2² =", result_sigmasq2, "\n")

4. 结果验证 #

将求解结果代入原方程组验证残差：

# 验证F1和F2的残差
mu2 <- result_mu2
sigmasq2 <- result_sigmasq2
F1_check <- log(((exp(2*mu_1+sigmasq_1)*(exp(sigmasq_1)-1))+(exp(2*mu2+sigmasq2)*(exp(sigmasq2)-1)))/
                  ((exp(mu_1+sigmasq_1/2)+exp(mu2+sigmasq2/2))^2)+1) - sigmasq_Z
F2_check <- log(exp(mu_1+sigmasq_1/2)+exp(mu2+sigmasq2/2)) - sigmasq_Z/2 - mu_Z

cat("F1残差：", F1_check, "\n")
cat("F2残差：", F2_check, "\n")

残差应接近0，说明求解正确。

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