Многомерная интеграция — связанные ограничения

Мне нужно вычислить значение многомерного интеграла в С++. Я нашел множество библиотек, способных решить эту задачу для фиксированных предельных интегралов,

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

Однако интегралы, на которые я смотрю, имеют переменные пределы,

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

Чтобы прояснить, что я имею в виду, вот наивная реализация 2D-суммы Римана в 2D, которая возвращает желаемый результат,

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 — общий верхний предел интегрирования. Хотя эта реализация работает, она медленная и, следовательно, непрактичная для более высоких измерений.

Существуют эффективные алгоритмы для более высоких измерений, но, насколько мне известно, реализации библиотек (например, Куба) принимают фиксированный вектор значений в качестве предельного аргумента, что делает их бесполезными для моей проблемы.

Есть ли какая-то причина для этого и / или есть ли какой-нибудь трюк, чтобы обойти проблему?


person emher    schedule 30.05.2014    source источник
comment
И у вас тоже был вопрос? (кроме ОТ)   -  person πάντα ῥεῖ    schedule 30.05.2014
comment
Да, мой вопрос в том, есть ли какие-либо эффективные библиотечные процедуры и/или какие-либо хитрости для адаптации стандартных методов к моей проблеме. Реализация, представленная в ответе, работает, но непрактична для более высоких измерений. Я немного перефразировал вопрос, чтобы прояснить этот момент.   -  person emher    schedule 30.05.2014
comment
Какое отношение тег [cuba] имеет к чему-либо?   -  person Tom Fenech    schedule 30.05.2014
comment
Именно это я имел в виду под помимо Off Topic ( см. 5.)   -  person πάντα ῥεῖ    schedule 30.05.2014
comment
Куба является примером библиотеки, которая решает проблему для фиксированных пределов интеграции (но, насколько мне известно, не для связанных пределов). feynarts.de/cuba   -  person emher    schedule 30.05.2014
comment
@user22496 Тег cuba здесь гласит: Cuba — это микрофреймворк Ruby для веб-разработки.   -  person πάντα ῥεῖ    schedule 30.05.2014
comment
Ах, я этого не видел; Я удалил тег Куба. Также я обновил вопрос, чтобы направить его в сторону от ОТ, как определено в (5).   -  person emher    schedule 30.05.2014
comment
Вы пытались оптимизировать свой код? Я вижу слишком много ненужных умножений, даже смешанных fp/int, которые часто медленные, также вызов f(x) медленнее, чем просто скопировать его код внутрь...   -  person Spektre    schedule 30.05.2014


Ответы (3)


Ваш порядок интеграции неверен, должен быть dy dx.


Вы интегрируете по треугольнику

0 <= x <= y <= L

внутри квадрата [0,L]x[0,L]. Это можно смоделировать путем интегрирования по полному квадрату, где подынтегральная функция f определяется как 0 вне треугольника. Во многих случаях, когда f определено на полном квадрате, это можно сделать, взяв произведение f на индикаторную функцию треугольника в качестве нового интеграла.

person Lutz Lehmann    schedule 07.06.2014

При интегрировании по треугольной области, такой как 0<=x<=y<=L, можно воспользоваться симметрией: проинтегрировать f(min(x,y),max(x,y)) по квадрату 0<=x,y<=L и разделить результат на 2. Это имеет преимущество перед расширением f на ноль (метод, упомянутый Лутцем) в том, что расширенная функция является непрерывным, что повышает производительность процедуры интегрирования.

Я сравнил их на примере интеграла от 2x+y по 0‹=x‹=y‹=1. Истинное значение интеграла равно 2/3. Давайте сравним производительность; для демонстрационных целей я использую подпрограмму Matlab, но это не относится к используемому языку или библиотеке.

Расширение на ноль

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 раз точнее первого.

К сожалению, этот трюк с симметрией недоступен для общих доменов; но интегрирование по треугольнику встречается достаточно часто, поэтому полезно помнить об этом.

person Community    schedule 01.03.2015

Я был немного смущен вашим интегральным определением, но из вашего кода я вижу это так:

интегральный

только что провел некоторое тестирование, вот ваш код:

//---------------------------------------------------------------------------
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 раза (на 32битном приложении на WOW64 AMD 3.2ГГц)
  • для более высоких измерений он будет умножаться
  • но все же я думаю, что этот подход медленный
  • The only thing to reduce complexity I can think of is algebraically simplify things
    • either by integration tables or by Laplace or Z transforms
    • но для этого нужно знать f(*x)...
  • constant time reduction can of course be done
    • by the use of multi-threading
    • и/или использование графического процессора
    • это может дать вам увеличение скорости в N раз
    • потому что это все напрямую распараллеливается
person Spektre    schedule 30.05.2014