在R中通过插值推导近似函数并求解累积和最大化的N值

嘿，我来帮你解决这个问题！首先得提一句：你给的示例数据里x有6个值，但y只有5个，这会直接导致拟合失败，我先补个合理的y值（比如最后加个1200）来演示，你替换成自己的完整数据就行。

1. 实现你要的findfun()：生成五次多项式拟合函数 #

在R里，五次多项式拟合其实可以用线性模型lm()来做，因为多项式的系数求解本质是线性回归问题。我给你封装好一个findfun()函数，完全符合你的需求：

findfun <- function(x, y, pow = 5) {
  # 先检查x和y长度是否匹配，避免报错
  if (length(x) != length(y)) {
    stop("喂，x和y的长度得一样啊！")
  }
  # 拟合多项式回归模型，raw=TRUE表示用原始的x^k项，不是正交多项式
  model <- lm(y ~ poly(x, pow, raw = TRUE))
  # 提取系数：截距项是x^0的系数，然后依次是x^1到x^5的系数
  coeffs <- coef(model)
  # 生成拟合函数，输入新x就能返回拟合值
  fit_fun <- function(x_new) {
    # 计算每一项的取值然后求和
    terms <- sapply(0:pow, function(i) coeffs[i+1] * (x_new^i))
    rowSums(terms)
  }
  return(fit_fun)
}

用你的数据测试一下： #

# 修正后的示例数据（x和y长度一致）
x <- c(0, 1, 2, 3, 4, 5)
y <- c(200, 320, 455, 612, 899, 1200)

# 生成五次多项式函数
approxfun <- findfun(x, y, pow = 5)

# 试一下x=2的时候，应该接近455
approxfun(2)

这个approxfun就是你要的五次多项式函数，形式就是f(x) = a*x^5 + b*x^4 + c*x^3 + d*x^2 + e*x + f（注意这里的截距项是公式里的最后那个f）。

2. 找最大N：让累积消耗不超过总苹果数 #

接下来是核心需求：给定总苹果数，找到最多能吃多少天，也就是找到最大的N，使得sum(f(1), f(2), ..., f(N)) ≤ 总苹果数。

思路解析： #

多项式的累积和是有公式的，比如对于f(x) = a0 + a1*x + a2*x² + ... + a5*x⁵，前N项的和可以用求和公式直接计算，不用循环累加（效率更高）。然后我们用二分法来找最大的N，避免一个个试的麻烦。

完整实现代码： #

首先调整一下findfun，让它同时返回拟合函数和系数，这样计算累积和更准确：

findfun <- function(x, y, pow = 5) {
  if (length(x) != length(y)) {
    stop("x和y长度必须一致！")
  }
  model <- lm(y ~ poly(x, pow, raw = TRUE))
  coeffs <- coef(model)
  fit_fun <- function(x_new) {
    terms <- sapply(0:pow, function(i) coeffs[i+1] * (x_new^i))
    rowSums(terms)
  }
  # 返回包含函数和系数的列表
  return(list(fun = fit_fun, coeffs = coeffs))
}

# 计算前N天的累积消耗量
calc_cumulative <- function(N, coeffs) {
  a0 <- coeffs[1]  # x^0的系数（截距）
  a1 <- coeffs[2]  # x^1的系数
  a2 <- coeffs[3]  # x^2的系数
  a3 <- coeffs[4]  # x^3的系数
  a4 <- coeffs[5]  # x^4的系数
  a5 <- coeffs[6]  # x^5的系数
  
  # 用多项式求和公式计算累积和
  sum_total <- a0*N + 
               a1*(N*(N+1)/2) + 
               a2*(N*(N+1)*(2*N+1)/6) + 
               a3*(N^2*(N+1)^2/4) + 
               a4*(N*(N+1)*(2*N+1)*(3*N^2+3*N-1)/30) + 
               a5*(N^2*(N+1)^2*(2*N^2+2*N-1)/12)
  return(sum_total)
}

# 找最大的N：累积消耗不超过总苹果数
find_max_days <- function(total_apples, fit_result) {
  # 先找一个足够大的上界，确保累积和超过总苹果数
  upper <- 1
  while (calc_cumulative(upper, fit_result$coeffs) < total_apples) {
    upper <- upper * 2
  }
  # 二分法查找最大N
  lower <- 0
  while (upper - lower > 1) {
    mid <- floor((lower + upper)/2)
    current_sum <- calc_cumulative(mid, fit_result$coeffs)
    if (current_sum <= total_apples) {
      lower <- mid
    } else {
      upper <- mid
    }
  }
  return(lower)
}

测试示例： #

假设总共有5000个苹果：

total_apples <- 5000
# 先得到拟合结果
fit_result <- findfun(x, y, pow=5)
# 找最大天数
max_days <- find_max_days(total_apples, fit_result)

# 输出结果
cat("苹果最多能维持", max_days, "天\n")
cat("前", max_days, "天总共吃了", round(calc_cumulative(max_days, fit_result$coeffs), 2), "个苹果\n")
cat("第", max_days+1, "天要吃", round(fit_result$fun(max_days+1), 2), "个，加上之前的就超量啦！\n")

几个小提醒： #

• 如果你每天的消耗量是从x=1开始的（x=0不算一天），那建议把你的x数据改成从1开始（比如x <- c(1,2,3,4,5,6)），避免x=0的初始值干扰结果。
• 五次多项式拟合需要至少6个数据点（因为有6个系数要解），如果你的数据点少于6个，会出问题，这时候可以降低多项式的次数，或者用样条插值（不过你明确要五次多项式，所以注意数据量）。

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