R语言中如何对长度不等的两个向量做数值除法求解目标积分

变量定义 #

• $\Delta G(z)$：300行2列数值矩阵，第一列为距离z的采样值，第二列为对应距离下的自由能
• $D(z)$：30行2列数值矩阵，第一列为距离z的采样值，第二列为对应距离下的扩散系数

原有分步拟合计算的方案误差极大，核心问题是两次独立拟合会引入双重拟合偏差，若未对齐采样区间、使用了不合适的拟合方法（如高次多项式），误差会在除法、积分步骤被进一步放大，完全偏离真实值。

正确求解流程（无额外拟合误差） #

核心思路：放弃分步拟合的思路，先通过保单调插值将两个离散数据集映射到统一的连续函数上，直接构造被积函数做数值积分，全程避免不必要的拟合引入的误差。

1. 预处理确定有效积分区间 #

积分仅在两个数据集都有实测值的z区间内计算，禁止无依据外推：

# 输入数据命名约定：
# DG：300*2矩阵，列1为z值，列2为对应ΔG值
# Dz：30*2矩阵，列1为z值，列2为对应D值
colnames(DG) <- c("z", "G")
colnames(Dz) <- c("z", "D")

# 取两个数据集z值的公共区间
z_min <- max(min(DG[, "z"]), min(Dz[, "z"]))
z_max <- min(max(DG[, "z"]), max(Dz[, "z"]))

2. 构造保单调连续插值函数 #

使用R自带stats包的保单调三次样条插值，既保证函数光滑可积，又不会出现高次多项式拟合的龙格震荡问题：

# 注意替换beta为你实际单位、温度下的β值，以下为298K、能量单位kJ/mol的示例
beta <- 1/(8.314e-3 * 298)

# 生成D(z)的插值函数
interp_D <- splinefun(
  x = Dz[, "z"], 
  y = Dz[, "D"], 
  method = "monoH.FC" # 保单调插值，避免非物理震荡
)

# 生成指数项exp(-βΔG(z))的插值函数
interp_expG <- splinefun(
  x = DG[, "z"],
  y = exp(-beta * DG[, "G"]),
  method = "monoH.FC"
)

3. 直接构造被积函数计算积分 #

不需要额外做函数拟合、除法运算，直接基于插值得到的连续函数做数值积分：

# 定义被积函数
integrand <- function(z) {
  interp_expG(z) / interp_D(z)
}

# 计算定积分，调高收敛阈值保证精度
calc_result <- integrate(
  integrand,
  lower = z_min,
  upper = z_max,
  rel.tol = 1e-8,
  subdivisions = 1000L
)

# 输出积分结果
print(calc_result$value)

关键注意事项 #

• 若积分区间内存在D(z)趋近于0的点，先核实数据有效性，避开数值奇点区间再计算，避免除以0导致的结果失真
• 不要使用普通线性插值：线性插值的函数一阶导数不连续，会降低积分精度
• 不要使用高次多项式全局拟合：全局拟合很容易在采样点间隙出现非物理的震荡，直接导致积分结果错误

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