非恒定间隔下如何用有限差分计算向量的二阶导数？

你的判断没错——原代码的误差根源在于一阶导数dydx对应的是相邻x的中点，而后续计算时用了错位的x序列，导致二阶导数的计算逻辑完全不匹配非均匀网格的特性。针对非均匀间隔的x，必须使用基于拉格朗日插值推导的三点公式来计算二阶导数，才能保证精度。

解决方案：三点公式直接计算二阶导数 #

对于非均匀网格中的每个内部点x[i]（2 ≤ i ≤ n-1），二阶导数的计算公式由三点拉格朗日插值多项式求导推导而来：

d2ydx2[i] = 2 * [ y[i+1]/(h2*(h1+h2)) - y[i]/(h1*h2) + y[i-1]/(h1*(h1+h2)) ]

其中h1 = x[i] - x[i-1]，h2 = x[i+1] - x[i]。

向量化实现（高效适合大样本） #

用向量化操作替代循环，避免性能损耗：

set.seed(123) # 固定随机种子方便复现
x <- rnorm(2000)
y <- x^2

# 排序x和对应y（保持x递增）
ord <- order(x)
x <- x[ord]
y <- y[ord]

n <- length(x)
d2ydx2 <- rep(NA, n)

# 向量化计算各部分项
h1 <- diff(x)[-length(diff(x))]  # x[i]-x[i-1]，对应i=2到n-1
h2 <- diff(x)[-1]                # x[i+1]-x[i]，对应i=2到n-1
term1 <- y[-c(1, 2)] / (h2 * (h1 + h2))
term2 <- -y[-c(1, n)] / (h1 * h2)
term3 <- y[-c(n-1, n)] / (h1 * (h1 + h2))

# 填充内部点的二阶导数
d2ydx2[2:(n-1)] <- 2 * (term1 + term2 + term3)

# 绘图验证（理论二阶导数为2）
plot(x, d2ydx2, main = "非均匀网格下的二阶导数", ylab = "d2y/dx²")
abline(h = 2, col = "red", lwd = 2)

循环实现（逻辑更直观） #

如果需要更清晰的逻辑展示，循环版本也能满足需求：

set.seed(123)
x <- rnorm(2000)
y <- x^2

ord <- order(x)
x <- x[ord]
y <- y[ord]

n <- length(x)
d2ydx2 <- rep(NA, n)

for (i in 2:(n-1)) {
  h1 <- x[i] - x[i-1]
  h2 <- x[i+1] - x[i]
  d2ydx2[i] <- 2 * (y[i+1]/(h2*(h1+h2)) - y[i]/(h1*h2) + y[i-1]/(h1*(h1+h2)))
}

plot(x, d2ydx2, main = "非均匀网格下的二阶导数", ylab = "d2y/dx²")
abline(h = 2, col = "red", lwd = 2)

为什么原代码会出错？ #

原代码中dydx的每个值对应x[i]和x[i+1]的中点，而diff(x[-1])对应的是x[2]到x[n]的间隔，这两个序列的位置完全错位。当x间隔变化较大时，这种错位计算会产生极大的误差，完全不符合二阶导数的数学定义。

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