首页 > 解决方案 > 计算R中函数乘积的双重积分的问题

问题描述

我正在尝试使用 R 包“pracma”中的函数“integral2”计算双积分。我在计算时遇到问题

integral2(function(x,y){ X(x)*R(x,y)*X(y) }, 0, 10, 0, 10)

在哪里

X <- function(t) {
  -0.4*sqrt(2)*sin(pi*1*t)+0.016*sqrt(2)*sin(pi*2*t)-0.01*sqrt(2)*sin(pi*3*t)
}

R <- function(x,y){(1/2*(x^2-x+1/6))*(1/2*(y^2-y+1/6))-
           (1/24*((abs(x-y)^4)-2*(abs(x-y)^3)+(abs(x-y)^2)-1/30))}.

我对 r 中的双积分的结果是

integral2(function(x,y){ X(x)*R(x,y)*X(y) }, 0, 10, 0, 10)$Q = 80.77929, 

但如果我在 Maple 中计算相同的积分,结果是 87.911。

标签: r

解决方案


以下将起作用。

X <- function(t){-0.4*sqrt(2)*sin(pi*1*t)+0.016*sqrt(2)*sin(pi*2*t)-
    0.01*sqrt(2)*sin(pi*3*t)}

R <- function(x,y){(1/2*(x^2-x+1/6))*(1/2*(y^2-y+1/6))-
    (1/24*((abs(x-y)^4)-2*(abs(x-y)^3)+(abs(x-y)^2)-1/30))}

f <- function(x, y){X(x)*R(x, y)*X(y)}

integral2b <- function(f, lower, upper){
  integrate( function(y) {
    sapply(y, function(y) {
      integrate(function(x) f(x,y), lower[1], upper[1])$value
    })
  }, lower[2], upper[2])
}

integral2b(f, c(0, 0), c(10, 10))
#84.94517 with absolute error < 0.0081

请参阅R-Help,此答案改编自该线程。


推荐阅读