多维整合 - 耦合极限

Question

我需要在 C++ 中计算高维积分的值。我找到了许多能够解决固定极限积分任务的库，

\int_{0}^{L} \int_{0}^{L} dx dy f(x,y) .

但是我正在查看的积分有可变限制，

\int_{0}^{L} \int_{x}^{L} dx dy f(x,y) .

为了澄清我的意思，这是一个简单的二维黎曼和实现，它返回所需的结果，

int steps = 100;
double integral = 0;
double dl = L/((double) steps);
double x[2] = {0};

for(int i = 0; i < steps; i ++){
    x[0] = dl*i;
    for(int j = i; j < steps; j ++){
        x[1] = dl*j;
        double val = f(x);
        integral += val*val*dl*dl;
    }
}

其中 f 是某个任意函数，L 是常见的积分上限。虽然此实现有效，但速度很慢，因此对于更高维度不切实际。

存在用于更高维度的有效算法，但据我所知，库实现（例如 Cuba）采用固定值向量作为限制参数，这使得它们对我的问题无用。

这有什么原因和/或有什么技巧可以规避这个问题吗？

Answer 1

您的集成顺序错误，应该是dy dx。

---

你正在对三角形进行积分

0

在正方形 [0,L]x[0,L] 内。这可以通过对整个正方形进行积分来模拟，其中被积函数 f 在三角形外部定义为 0。在很多情况下，当 f 定义在整个正方形上时，这可以通过将 f 与三角形的指示函数的乘积作为新的被积函数来实现。

Answer 2

在三角形区域（例如 0<=x<=y<=L）上积分时，可以利用对称性：在正方形 0<=x,y<=L 上积分 f(min(x,y),max(x,y))，然后将结果除以 2。这比将 f 扩展为零（方法LutzL 提到的），因为扩展函数是连续的，这提高了集成例程的性能。

我在 2x+y 对 0

扩展为零

f = @(x,y) (2*x+y).*(x<=y);
result = integral2(f, 0, 1, 0, 1);
fprintf('%.9f\n',result);

输出：

警告：已达到最大功能评估次数 (10000)。结果未通过全局错误测试。

0.666727294

对称延伸

g = @(x,y) (2*min(x,y)+max(x,y));
result2 = integral2(g, 0, 1, 0, 1)/2;
fprintf('%.9f\n',result2);

输出：

0.666666776

第二个结果比第一个结果准确 500 倍。

不幸的是，这种对称技巧不适用于一般领域；但是三角形上的积分经常出现，所以记住它很有用。

Answer 3

我对你的积分定义有点困惑，但从你的代码中我看到它是这样的：

刚刚做了一些测试，所以这是你的代码：

//---------------------------------------------------------------------------
double f(double *x) { return (x[0]+x[1]); }
void integral0()
    {
    double L=10.0;
    int steps = 10000;
    double integral = 0;
    double dl = L/((double) steps);
    double x[2] = {0};
    for(int i = 0; i < steps; i ++){
        x[0] = dl*i;
        for(int j = i; j < steps; j ++){
            x[1] = dl*j;
            double val = f(x);
            integral += val*val*dl*dl;
        }
    }
    }
//---------------------------------------------------------------------------

这是优化后的代码：

//---------------------------------------------------------------------------
void integral1()
    {
    double L=10.0;
    int i0,i1,steps = 10000;
    double x[2]={0.0,0.0};
    double integral,val,dl=L/((double)steps);
    #define f(x) (x[0]+x[1])
    integral=0.0;
    for(x[0]= 0.0,i0= 0;i0<steps;i0++,x[0]+=dl)
    for(x[1]=x[0],i1=i0;i1<steps;i1++,x[1]+=dl)
        {
        val=f(x);
        integral+=val*val;
        }
    integral*=dl*dl;
    #undef f
    }
//---------------------------------------------------------------------------

结果：

[ 452.639 ms] integral0
[ 336.268 ms] integral1

• 因此速度提升约 1.3 倍（在 WOW64 AMD 3.2GHz 上的 32 位应用程序上）
• 对于更高的维度，它会成倍增加
• 但我仍然认为这种方法很慢
• 我能想到的唯一可以降低复杂性的方法是通过代数简化事物
  • 通过积分表或拉普拉斯变换或 Z 变换
  • 但为此必须知道 f(*x) ...
• 当然可以不断减少时间
  • 通过使用多线程
  • 和/或 GPU 使用情况
  • 这可以给你N倍的速度提升
  • 因为这一切都可以直接并行化