0

我有一个函数 f(x) = 1/(x + a+ b*I*sign(x)) 我想计算积分

dx dy dz f(x) f(y) f(z) f(x+y+z) f(xy - z)

在整个 R^3 上（b>0 和 a,- b 是统一的）。这只是一个有代表性的例子——实际上我有 n<7 个变量和 2n-1 个 f() 实例，其中 n 个涉及 n 个积分变量，其中 n-1 个涉及积分变量的一些线性组合。在这个阶段，我只对相对误差为 1e-3 左右的粗略估计感兴趣。

我尝试了以下库：

  • Steven Johnson 的 cubature 代码：hcubature 算法可以工作，但速度非常慢，即使 n=2 也需要数亿次被积函数。
  • HintLib：我尝试了与 Genz-Malik 规则、体积例程、VEGAS 和 MISER 与 Mersenne twister RNG 的自适应集成。对于 n = 3，只有第一个似乎是可行的选择，但对于 n = 3 和 relerr = 1e-2，它再次需要数亿次被积函数评估，这并不令人鼓舞。

对于积分区域，我尝试了两种方法：在 [-200, 200]^n 上进行积分（即一个区域如此之大以至于它基本上捕获了大部分积分）和替换 x = sinh(t) 这似乎是标准把戏。

我在数值分析方面没有太多经验，但大概困难在于 sign() 术语的不连续性。对于 n=2 和 f(x)f(y)f(xy)，沿 x=0、y=0、x=y 存在不连续性。这些在原点周围产生了一个非常尖锐的峰值（在各个象限中具有不同的符号），并且在 x=0,y=0,x=y 处形成了一种“脊”，其中被积函数的绝对值很大并且符号变化为你越过他们。所以至少我知道哪些地区很重要。我在想也许我可以做蒙特卡洛，但以某种方式提前“告诉”算法将注意力集中在哪里。但我不太确定该怎么做。

如果您对如何以合理的计算能力评估积分或如何使我的蒙特卡洛“想法”发挥作用有任何建议，我将不胜感激。我已经坚持了一段时间，所以欢迎任何意见。提前致谢。

4

1 回答 1

0

您可以做的一件事是使用指导函数进行蒙特卡罗积分：给定 ∫ f ( x ) dx 的积分（为了简单起见，我将其写成一维），将其写为 ∫ f ( x )/ g ( x ) g ( x ) dx，并g(x)用作您从中采样的分布x。

由于g(x)是任意的，因此构造它使得 (1) 它具有您期望它们在 中的峰值f(x)，并且 (2) 使得您可以从中采样x（g(x)例如，高斯或1/(1+x^2)）。

或者，您可以使用 Metropolis 型马尔可夫链 MC。它会（几乎）自己找到被积函数的相关区域。 这里有几个简单的例子。

于 2013-05-03T12:25:09.877 回答