代码之家 › 专栏 › 技术社区 › jay.sf

如何将梯形积分法调整为自定义零点?

numerical-integration r

jay.sf · 技术社区 · 7 年前

我要计算序列向量的积分。由于没有可用的功能,我使用 trapezoidal 方法 ¹ .

iglTzm <- function(x, y) sum(diff(x) * (head(y, -1) + tail(y, -1))) / 2

序列的第一个元素应该是零点,所以原理是:如果序列的值主要低于第一个值,那么积分应该是负的,否则是正的,或者是0。

考虑矩阵 m1 :

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    6    7    8    8    6    8   10
[2,]    9    9    8    9    9    8    9
[3,]    9   10   10    9    9    9    9
[4,]    9    8    8    8    6    8    9
[5,]   10   10   10    9   10    8    0
[6,]    9    8    9   10    9    9    9

与这些原始值的集成很可能导致不一致的值:

> setNames(apply(m1, 1, iglTzm, 0:6), 1:6)
  1   2   3   4   5   6 
 15   2  -2   7 -52   0

所以我调整第一个值(第1列)上的序列(行),以便设置正确的符号,并得到矩阵 m2 :

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    1    2    2    0    2    4
[2,]    0    0   -1    0    0   -1    0
[3,]    0    1    1    0    0    0    0
[4,]    0   -1   -1   -1   -3   -1    0
[5,]    0    0    0   -1    0   -2  -10
[6,]    0   -1    0    1    0    0    0

从逻辑上讲,这不会改变任何价值观 iglTzm() 抛出,因为 diff() 是一样的:

> setNames(apply(m2, 1, iglTzm, 0:6), 1:6)
  1   2   3   4   5   6 
 15   2  -2   7 -52   0

不管怎样,因为我不能简单地缩放或反转它,所以我还没有一个聪明的想法,如何调整函数以获得正确的符号,这是很明显的:

#  1   2   3   4   5   6 
# 15  -2   2  -7 -52   0

有人知道怎么适应吗 IGLTZ() 得到带正确符号的积分?

情节 m2 应该更详细地说明这个原理:

数据

m1 <- matrix(c(6, 7, 8, 8, 6, 8, 10,
                9, 9, 8, 9, 9, 8, 9,
                9, 10, 10, 9, 9, 9, 9,
                9, 8, 8, 8, 6, 8, 9, 
                10, 10, 10, 9, 10, 8, 0, 
                9, 8, 9, 10, 9, 9, 9), 6, byrow=TRUE)

m2 <- t(apply(m1, 1, function(x) scale(x, center=x[1], scale=FALSE)))

# plot
par(mfrow=c(2, 3))
lapply(1:nrow(m2), function(x) {
  plot(m2[x, ], type="l", main=x)
  abline(h=m2[x, 1], col="red", lty=2)
})

1 回复 | 直到 7 年前

Julius Vainora 7 年前

首先,还有一个小问题,但更重要,尽管在解决它之后,您的问题仍然有效。我的意思是 x 和 y 由于函数在 apply .

但这还不够,现在我们回到你的问题上来。为此,让我们回顾一下通常的集成:_f(x)dx(从a到b的限制)将集成f以下的区域,这就是您已经成功完成的功能。现在你想要的是调整它的水平。但是如果我们从A到B整合,这与(f(x)-f(a))dx=f(x)dx-(b-a)f(a)相同,这导致

iglTzm <- function(y, x) sum(diff(x) * (head(y, -1) + tail(y, -1))) / 2 - y[1] * (max(x) - min(x))
setNames(apply(m1, 1, iglTzm, 0:6), 1:6)
#  1  2  3  4  5  6 
#  9 -2  2 -7 -8  0

只有两个绝对值不同于 X 和 Y 是颠倒的。我们来看看第一个函数:应该是9还是15?我们有2*2/2+1*2+1*2/2+2*4/2=9,所以我们确实想反转 X 和 Y .

编写函数的另一种方法是

iglTzm <- function(y, x) sum(diff(x) * (head(y - y[1], -1) + tail(y - y[1], -1))) / 2
setNames(apply(m1, 1, iglTzm, 0:6), 1:6)
#  1  2  3  4  5  6 
#  9 -2  2 -7 -8  0

编辑 :通过反转,我只表示函数定义中的顺序或您如何使用它 应用 ;函数本身的y(函数值)和x(网格值)是可以的。