Численное интегрирование разрывной функции в нескольких измерениях

Question

У меня есть функция 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 или около того.

Я пробовал следующие библиотеки:

• Кубатурный код Стивена Джонсона: алгоритм hcubature работает, но ужасно медленен, он берет сотни миллионов оценок подынтегральных функций даже при n=2.
• HintLib: я попробовал адаптивную интеграцию с правилом Генца-Малика, кубатурными процедурами, VEGAS и MISER с твистером RNG Мерсенна. Для n=3 только первый вариант выглядит несколько жизнеспособным, но он снова требует сотни миллионов оценок подынтегральных функций для n=3 и relerr = 1e-2, что не внушает оптимизма.

Для области интеграции я испробовал оба подхода: интегрирование по [-200, 200]^n (т. Е. Область настолько большая, что она по существу захватывает большую часть интеграла) и подстановка x = sinh(t), которая кажется стандартный трюк.

У меня нет большого опыта в численном анализе, но, вероятно, трудность заключается в несплошности от знака (). Для n=2 и f(x)f(y)f(xy) имеются разрывы вдоль x=0, y=0, x=y. Они создают очень острый пик вокруг начала координат (с другим знаком в различных квадрантах) и своего рода "гребни" при x = 0, y = 0, x = y, вдоль которых подынтегральное выражение является большим по абсолютной величине и меняет знак как Вы пересекаете их. Так что, по крайней мере, я знаю, какие регионы важны. Я думал, что, возможно, я мог бы сделать Монте-Карло, но каким-то образом заранее "сказать" алгоритму, на чем сосредоточиться. Но я не совсем уверен, как это сделать.

Я был бы очень признателен, если бы у вас был какой-либо совет о том, как оценить интеграл с разумным количеством вычислительной мощности или как заставить мою идею Монте-Карло работать. Я застрял на этом какое-то время, поэтому любые отзывы будут приветствоваться. Заранее спасибо.

Answer 1

Одна вещь, которую вы можете сделать, это использовать направляющую функцию для вашей интеграции Монте-Карло: учитывая интеграл (я пишу его в 1D для простоты) из ∫ f(x) dx, запишите его как ∫ f(x) /g(x) г(х) дх, и использовать g(x) как распределение, из которого вы выбираете x,

поскольку g(x) произвольно, сконструируйте его так, чтобы (1) он имел пики там, где вы ожидаете f(x)и (2) такой, что вы можете попробовать x от g(x) (например, гауссов или 1/(1+x^2)).

В качестве альтернативы вы можете использовать цепь Маркова типа Metropolis. Он найдет соответствующие регионы подынтегральной функции (почти) сам по себе. Вот пара тривиальных примеров.