R:如何获得用于分析积分的插值样条的分段系数?
R: How to get piecewise coefficients of an interpolation spline for analytical integration?
动机
我正在对一个深度嵌套的多重积分进行数值计算。在每个嵌套级别,我得到一个位于下面级别的积分向量,将其乘以密度函数向量,得到该级别的被积函数向量 y
。 x
值间隔不均匀。
被积函数是弯曲的,梯形积分不够准确,所以我想做一个允许曲率的积分。 Simpson 规则不适用,因为横坐标不是均匀分布的。所以我建议做一个三次样条插值,然后通过解析计算每个段的三次样条的积分来计算样条函数的积分
问题
我一直在研究 spline
和 splinefun
以及 splines2
包中的函数。但是我找不到任何可以告诉我一系列三次多项式的系数的信息 - 节点之间的每段一个。
如果有人能给我指点一个进行样条插值并提供三次系数数组的函数,我将不胜感激。
谢谢。
这是在这里扩展我的新答案的好机会:有了分段参数化,计算分段积分就很容易了。
这是一个用于计算的(矢量化)函数:
## a function for integration on a piece
piecewise_int <- function (hi, yi, bi, ci, di) {
yi * hi + bi * hi ^ 2 / 2 + ci * hi ^ 3 / 3 + di * hi ^ 4 / 4
}
下面我将以那个帖子里的小例子,展示如何整合spline。
## the small example in the linked thread
set.seed(0)
xk <- c(0, 1, 2)
yk <- round(runif(3), 2)
f <- splinefun(xk, yk, "natural") ## natural cubic spline
construction_info <- environment(f)$z
## information for integration
int_info <- with(construction_info,
list(h = diff(x), y = y[-n], b = b[-n], c = c[-n], d = d[-n])
)
## cubic spline integration on all pieces
integral <- sum(do.call(piecewise_int, int_info))
#[1] 0.81375
我们也可以进行数值积分来验证这个结果。
integrate(f, 0, 2)
#0.81375 with absolute error < 9e-15
动机
我正在对一个深度嵌套的多重积分进行数值计算。在每个嵌套级别,我得到一个位于下面级别的积分向量,将其乘以密度函数向量,得到该级别的被积函数向量 y
。 x
值间隔不均匀。
被积函数是弯曲的,梯形积分不够准确,所以我想做一个允许曲率的积分。 Simpson 规则不适用,因为横坐标不是均匀分布的。所以我建议做一个三次样条插值,然后通过解析计算每个段的三次样条的积分来计算样条函数的积分
问题
我一直在研究 spline
和 splinefun
以及 splines2
包中的函数。但是我找不到任何可以告诉我一系列三次多项式的系数的信息 - 节点之间的每段一个。
如果有人能给我指点一个进行样条插值并提供三次系数数组的函数,我将不胜感激。
谢谢。
这是在这里扩展我的新答案的好机会:
这是一个用于计算的(矢量化)函数:
## a function for integration on a piece
piecewise_int <- function (hi, yi, bi, ci, di) {
yi * hi + bi * hi ^ 2 / 2 + ci * hi ^ 3 / 3 + di * hi ^ 4 / 4
}
下面我将以那个帖子里的小例子,展示如何整合spline。
## the small example in the linked thread
set.seed(0)
xk <- c(0, 1, 2)
yk <- round(runif(3), 2)
f <- splinefun(xk, yk, "natural") ## natural cubic spline
construction_info <- environment(f)$z
## information for integration
int_info <- with(construction_info,
list(h = diff(x), y = y[-n], b = b[-n], c = c[-n], d = d[-n])
)
## cubic spline integration on all pieces
integral <- sum(do.call(piecewise_int, int_info))
#[1] 0.81375
我们也可以进行数值积分来验证这个结果。
integrate(f, 0, 2)
#0.81375 with absolute error < 9e-15