多维不连续函数的数值积分



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

的积分

dxdy dz f (x) (y) f (z) f (x + y + z) f (x - y - z)

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

我尝试了以下库:

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

对于积分区域,我尝试了两种方法:对[- 200,200]^n进行积分(即一个非常大的区域,它基本上捕获了大部分积分)和替换x = sinh(t),这似乎是一个标准的技巧。

我没有太多的数值分析经验,但大概困难在于符号()项的不连续。对于n=2和f(x)f(y)f(x-y)在x=0, y=0, x=y处存在不连续点。它们在原点周围形成了一个非常尖锐的峰值(在不同的象限上有不同的符号),在x=0,y=0,x=y处形成了一些"脊",沿着这些"脊",被积函数的绝对值很大,当你穿过它们时,符号就会改变。至少我知道哪些区域是重要的。我在想,也许我可以用蒙特卡罗算法,但以某种方式提前"告诉"算法聚焦在哪里。但我不太确定该怎么做。

如果您有任何关于如何用合理的计算能力评估积分或如何使我的蒙特卡洛"想法"工作的建议,我将非常感激。我已经在这个问题上纠结了一段时间,所以欢迎任何意见。

你可以做的一件事是为你的蒙特卡罗积分使用一个引导函数:给定一个积分(为了简单起见,我用1D来写)∫f(x) dx,写成∫f <我> (<我> x ) g/<我> (<我> x ) g <我> (<我> x ) dx,并使用g(x)作为一个分布的样本x

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

或者,你可以使用Metropolis-type Markov chain MC。它会(几乎)自己找到被积函数的相关区域。下面是几个简单的例子:

最新更新