在R中绘制非线性函数F2=F1曲线的求解与绘图求助

核心问题原因 #

1. 定义域超出有效范围：原代码中F1_values取到100，但当F1接近R1=100时，方程无实数解，导致uniroot报错返回NA。
2. 单根求解限制：uniroot仅能找到区间内的一个根，无法处理一个F1对应多个F2解的场景。
3. 数值稳定性问题：当s=1时，原函数中的指数项会出现0次幂，直接计算易引发数值异常，需特殊简化处理。

修正后的代码实现 #

步骤1：加载必要包并定义常数 #

rm(list=ls())
library(ggplot2)
library(rootSolve) # 用于多根求解

# 固定常数
R1 <- 100
R2 <- 100
s <- 1
m <- 1

步骤2：定义目标函数（含特殊情况处理） #

针对s=1的场景单独简化函数，避免数值问题；s≠1时保留原形式并添加定义域校验：

# 定义关于F2的目标函数，输入为F2向量，返回函数值
target_func <- function(F2, F1, R1, R2, s, m) {
  if (s == 1) {
    # s=1时的简化形式
    term1 <- F1 / (F2^m)
    term2 <- m * ((R1 - F1) + (R2 - F2)) / (F1^m + F2^m)
    term1 - term2
  } else {
    # 通用形式，校验定义域有效性
    if (F1 >= R1 || any(F2 >= R2) || F1 <= 0 || any(F2 <= 0)) {
      return(rep(Inf, length(F2)))
    }
    term1 <- (F1 / (F2^m)) * ((R1 - F1)^((1 - s)/s))
    term2 <- m * ((R1 - F1)^(1/s) + (R2 - F2)^(1/s)) / (F1^m + F2^m)
    term1 - term2
  }
}

步骤3：求解每个F1对应的所有F2解 #

先确定F1的有效范围（以s=1,m=1为例，F1最大值约为41.42），再对每个F1用multiroot查找所有根：

# 确定F1的有效范围（避免无解方程）
max_F1 <- if (s ==1 && m==1) sqrt(200*R2) - R2 else R1*0.9 # 通用场景取R1的90%
F1_values <- seq(0.01, max_F1, by = 0.1)

# 收集所有有效解
solution_list <- lapply(F1_values, function(F1) {
  # 设定F2的初始猜测，覆盖可能的根区间
  initial_guesses <- c(0.1, R2*0.5)
  tryCatch({
    # 用multiroot求解多个根
    roots <- multiroot(f = target_func, start = initial_guesses, 
                       F1 = F1, R1 = R1, R2 = R2, s = s, m = m)$root
    # 过滤超出定义域的无效解
    valid_roots <- roots[roots > 0.01 & roots < R2 - 0.01]
    data.frame(F1 = rep(F1, length(valid_roots)), F2 = valid_roots)
  }, error = function(e) {
    # 无有效解时返回空数据框
    data.frame(F1 = numeric(0), F2 = numeric(0))
  })
})

# 合并所有解为数据框
data <- do.call(rbind, solution_list)

步骤4：绘制曲线 #

ggplot(data, aes(x = F1, y = F2)) +
  geom_line() +
  labs(title = "F2 as a Function of F1",
       x = "F1",
       y = "F2") +
  theme_minimal() +
  xlim(0, max_F1) +
  ylim(0, R2)

关键改进点 #

• 多根求解：使用rootSolve::multiroot替代uniroot，支持捕获一个F1对应的多个F2解。
• 定义域校验：提前过滤超出有效范围的F1/F2值，避免数值异常。
• 特殊场景优化：针对s=1的情况简化函数，提升数值稳定性。
• 有效范围缩小：根据方程特性调整F1的取值范围，减少无解方程的场景。

内容的提问来源于stack exchange，提问作者John M. Riveros