求解单未知量非线性方程组：技能技术行业边界定位

问题描述 #

我正在劳动市场任务模型（Task-Based model）中确定高/低技能技术行业的边界，该边界对应两类生产价格相等的任务I，满足：

( P1 - P2 = 0 )
其中：

• ( P1 = \mu_L[i]^{-1} \times (w_L/\delta)^\delta )
• ( P2 = \mu_H[i]^{-1} \times (w_H/\alpha)^\alpha )

这是动态联立问题：工资依赖于I，同时I也依赖于工资。相关参数与工资方程（基于实际数据设定）如下：

alpha <-  0.55 
delta <- 0.83
L_L <- 0.8  # 低技能劳动力供给
L_H <- 0.2  # 高技能劳动力供给
Y <- 40000  # 总产出
I_start <- 0.5
I <- I_start
   
w_H <- alpha * I * (Y / L_H) 
w_L <- delta*(1-I)*(Y / L_L)

I的定义是使低技能生产价格等于高技能生产价格的任务i，其中( \mu_H = i )，( \mu_L = 1-i )，i是0到1之间均匀分布的任务值，需要找到满足条件的i（即I）。

现有尝试的问题 #

1. 遍历法：
   步长过大时无法找到精确匹配，步长过小（如1e-8）会触发步数过多错误；添加容差筛选也未得到预期结果，代码如下：

   range2 <- seq(0.01,1, by = 0.001)
   range <- (1 - range2)
   # 精确匹配无结果
   range[which(range^(-1)*(w_L/delta)^delta == range2^(-1)*(w_H/alpha)^alpha)]
   # 容差尝试
   r <- range^(-1)*(w_L/delta)^delta - range2^(-1)*(w_H/alpha)^alpha
   rr <- ifelse(r<0.01, print(TRUE), print(FALSE))
   data.r <- data.frame(r, rr)
   length(unique(data.r$rr))
   

2. nleqslv包求解非线性方程组：
   错误地将mu_L和mu_H设为两个独立未知量，导致解超出[0,1]范围，代码如下：

   library(nleqslv)
   dslnex <- function(x) {
     y <- numeric(2)
     y[1] <- (x[1]^(-1)*(w_L/delta)^delta)
     y[2] <- (x[2]^(-1)*(w_H/alpha)^alpha)
     y 
   }
   xstart=c(0.5, 0.5)
   nleqslv(xstart, dslnex, control = list())
   

解决方案 #

核心思路：转化为单变量非线性方程求解 #

因为( \mu_L = 1 - \mu_H = 1 - I )，所有变量可统一用I表示，原条件转化为仅含I的单方程：

( (1-I)^{-1} \times (w_L/\delta)^\delta - I^{-1} \times (w_H/\alpha)^\alpha = 0 )
其中w_L和w_H本身也是I的函数，代入后用单变量求解器即可解决。

实现代码 #

方法1：使用R内置uniroot函数（推荐）

# 定义参数
alpha <- 0.55 
delta <- 0.83
L_L <- 0.8
L_H <- 0.2
Y <- 40000

# 定义目标函数：输入I，返回P1-P2的差值
target_func <- function(I) {
  # 计算当前I对应的工资
  w_H <- alpha * I * (Y / L_H)
  w_L <- delta * (1 - I) * (Y / L_L)
  # 计算生产价格差值
  P1 <- (1 - I)^(-1) * (w_L / delta)^delta
  P2 <- I^(-1) * (w_H / alpha)^alpha
  return(P1 - P2)
}

# 用uniroot求解，设定区间避免0/1处的奇点
result <- uniroot(target_func, interval = c(0.01, 0.99))
# 输出结果
cat("求解得到的边界I值为：", round(result$root, 4), "\n")

方法2：优化遍历法（粗搜+线性插值）

先粗遍历找到解所在区间，再用线性插值提高精度，避免步长过小的效率问题：

# 粗遍历，步长0.001
range_I <- seq(0.01, 0.99, by = 0.001)
# 计算每个I对应的P1-P2值
diff_vals <- sapply(range_I, target_func)
# 找到符号变化的区间（解在该区间内）
sign_change_idx <- which(diff(sign(diff_vals)) != 0)
# 提取区间端点
I_left <- range_I[sign_change_idx]
I_right <- range_I[sign_change_idx + 1]
diff_left <- diff_vals[sign_change_idx]
diff_right <- diff_vals[sign_change_idx + 1]
# 线性插值计算精确解
I_exact <- I_left - diff_left * (I_right - I_left) / (diff_right - diff_left)
cat("插值得到的边界I值为：", round(I_exact, 4), "\n")

关键说明 #

• 原问题是单变量联立问题：mu_L和mu_H由同一个I决定，不能作为独立未知量，这是之前nleqslv求解出错的核心原因。
• uniroot要求函数在区间端点符号相反，因此设定区间为(0.01, 0.99)，避开I=0或I=1时的分母为0奇点。
• 优化遍历法通过"粗搜区间+插值"平衡了效率与精度，适合对求解过程有直观需求的场景。

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