0

我知道我之前问过同样的问题，但是由于我在这里很新，所以这个问题被问得很糟糕而且无法重现。因此，我尝试在这里做得更好。（如果我只编辑旧的可能没有人会读它）

我有这个想要积分的双积分：这是一张图片

在此处输入图像描述

ff<-function(g,t) exp((16)*g)*exp(-8*t-(-t-0.01458757)^2/(0.0001126501))

integrate(Vectorize(function(t) integrate(function(g) 
                                          ff(g,t), -2.5,0)$value), -2, 2)

在 R 中运行它会给我错误：

  the integral is probably divergent

当我尝试在 Wolfram 中运行 sam 函数时，它给了我一个合适的值：（我不得不切换 g=x 和 t=y）

关联：

http://www.wolframalpha.com/input/?i=integration+[%2F%2Fmath%3Aexp%28%2816%29*x%29*exp%28-8*y-%28-y-0.01458757%29 ^2%2F%280.0001126501%29%29%2F%2F]+[%2F%2Fmath%3Adx+dy%2F%2F]+for+x+from+[%2F%2Fmath%3A-2.5%2F%2F] +to+[%2F%2Fmath%3A0%2F%2F]+for+y+从+[%2F%2Fmath%3A-2%2F%2F]+to+[%2F%2Fmath%3A2%2F%2F]

如您所见，它得到了一个有限的结果，有人可以在这里帮助我吗？

在此处输入图像描述

我在定义的区域上绘制了函数，但找不到奇点问题。看：

library('Plot3D')
x <- seq(-2.5,0, by = 0.01) #to see the peak change to: seq(-0.2,0, by = 0.001)
y <- seq(-2,2, by = 0.01) #"": seq(-0.1,0.1, by = 0.001)
grid <- mesh(x,y) 
z <- with(grid,exp((16)*x)*
  exp(-8*y-(-0.013615734-y-0.001+0.5*0.007505^2*1)^2/(2*0.007505^2)))
persp3D(z = z, x = x, y = y)

感谢您的帮助，我希望问题的结构比旧问题更好。

4

3 回答 3

2

还值得注意的是，在Integrate.c源文件中，错误信息的描述是

error messages
...
ier = 5 the integral is probably divergent, or
    slowly convergent. it must be noted that
    divergence can occur with any other value of ier.

因此，尽管消息说“可能发散”，但您的代码似乎更有可能缓慢收敛。

此外，您可以在收到此消息时继续运行并提取错误，如果您设置stop.on.error=FALSE

r <- integrate(Vectorize(function(t) 
    integrate(function(g) ff(g,t), -2.5,0)$value
), -2, 2, stop.on.error=FALSE); 
r$value

R 并不像 Mathematica 等 Wolfram 产品那样声称是一个花哨的数学求解器。它没有对积分进行任何象征性的简化，而这正是 Wolfram 多年来一直在完善的东西。如果您只是想对一堆双积分进行数值求解，那么像 Mathematica 或 Maple 这样的程序可能是更好的选择。这似乎并不是 R 花费大量开发资源的地方。

于 2014-09-08T05:48:06.117 回答
2

您的被积函数仅在 y=0 附近的小范围内显着非零。从?integrate

当在无限区间积分时，明确地这样做，而不是仅仅使用一个大数作为端点。这增加了正确答案的机会——任何在无限区间上的积分是有限的函数在该区间的大部分时间里都必须接近于零。

虽然您没有在无限区间上严格积分，但同样的数值问题也适用。事实上：

ff <- function(x, y)
exp(16*x - 8*y - (-y - 0.01458757)^2/0.0001126501)

f <- function(y)
integrate(ff, lower=-2.5, upper=0, y=y)$value

integrate(Vectorize(f), lower=-Inf, upper=Inf)
0.001323689 with absolute error < 4.4e-08

有趣的是，答案与从 Wolfram Alpha 获得的不同。我不确定在这里可以信任谁；一方面，我integrate多次使用 R 并且没有遇到问题（我可以说）；但是正如@MrFlick 所说，R 不是像 Wolfram Alpha 那样的专用数学求解器。

您还可以将rel.tol收敛参数设置为更严格的值，例如 1e-7 或 1e-8。这在内部积分中比外部积分更重要，因为前者的错误会传播到后者。在这种情况下，它对最终结果没有影响。

于 2014-09-08T06:31:42.583 回答
1

对于双积分，最好使用cubature包。

library(cubature)
f <- function(x){
  exp(16*x[1] - 8*x[2] - (x[2] + 0.01458757)^2/0.0001126501)
}

该hcubature函数给出的结果在减小容差时不稳定：

> hcubature(f, c(-2.5, -2), c(0,2))$integral
[1] 0.001285129
> hcubature(f, c(-2.5, -2), c(0,2), tol=1e-10)$integral
[1] 0.001293842

与 的结果相反pcubature，它是稳定的：

> pcubature(f, c(-2.5, -2), c(0,2))$integral
[1] 0.001323689
> pcubature(f, c(-2.5, -2), c(0,2), tol=1e-10)$integral
[1] 0.001323689

p 自适应版本 pcubature 反复将正交规则的度数加倍，直到实现收敛，并且基于 Clenshaw-Curtis 正交规则的张量积。对于少数 (<=3) 维的平滑被积函数，该算法通常优于 h 自适应积分，但对于更高维度或非平滑被积函数来说，这是一个糟糕的选择。

其次，RcppNumerical提供了强大的多重集成。

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include <RcppNumerical.h>
#include <cmath>
using namespace Numer;

class ValegardIntegrand: public MFunc
{
public:    
  double operator()(Constvec& x)
  {
    return exp(16*x[0] - 8*x[1] - pow(-x[1] - 0.01458757,2)/0.0001126501);
  }
};    

// [[Rcpp::export]]
Rcpp::List Valegard()
{
  ValegardIntegrand f;  
  Eigen::VectorXd lower(2);
  lower << -2.5, -2;
  Eigen::VectorXd upper(2);
  upper << 0.0, 2.0;
  double err_est;
  int err_code;
  double res = integrate(f, lower, upper, err_est, err_code, 
                         10000000, 1e-14, 1e-14);
  return Rcpp::List::create(
    Rcpp::Named("approximate") = res,
    Rcpp::Named("error_estimate") = err_est,
    Rcpp::Named("error_code") = err_code
  );
}

它给出了相同的结果pcubature：

> Valegard()
$approximate
[1] 0.001323689

$error_estimate
[1] 9.893371e-14

$error_code
[1] 0

顺便说一句，与 Mathematica 11 的精确集成也提供了这个结果：

Integrate[E^(16 x - 8877.04 (0.0145876 + y)^2 - 8 y), {y, -2, 2}, {x, -2.5, 0}]
0.00132369
于 2017-09-15T12:38:47.620 回答