метод Монте-Карло для вычисления многомерного интеграла

[Nobody]

Всем здравствуйте!
Возникла следующая задачка.
Имеется многомерный интеграл, скажем восьмимерный, в пределах от нуля до бесконечности. Проблема такая - какой из существующих математических пакетов может его вычислить? (инеграл явно не "берется").
и в каком источнике довольно понятно можно узнать про сам метод Монте-Карло?
спасибо заранее.

[ruan]

У нас о нем рассказывают в курсе теории вероятности. Можно поискать что-то и по численным методам.

[Руст]

Вряд ли методом Монте Карло можно вычислить этот интеграл эффективно, если не учитывать вид функции. Например, ваш интеграл скорее всего имеет вид преобразования Лапласа (раз от 0 до бесконечности). Соответственно эффективность вычисления существенно улучшится, если и сами точки будут разбросаны по Пуассоновому распределению, а интегральная функция заменена с учётом этого.

[maxal]

пакетов/библиотек много всяких.
GSL тот же: http://www.gnu.org/software/gsl/manual/gsl-ref_23.html

[Oleg S.]

Проблема сформулирована нечетко -
1) неизвестно есть ли особенности у интегрируемых функций ,
2) какова должна быть точность вычислений,
3) есть ли ограничения на время вычислений.

Я бы рекомендовал Вам посмотреть такие книги:
Соболь И.М. — Многомерные квадратные формулы и функции Хаара
Калиткин Н.Н. — Численные методы

Есть еще хорошая книга Соболя по методам Монте-Карло.

[sexstant]

Соболь И.М. Численные методы Монте Карло. М.:Наука. 1973.

Метод конечно хорош именно для многомерных интегралов, и был очень популярным на старых ЭВМ. Все упирается в датчик случайных чисел. Раньше это была электронная лампа. Так что ищите эффективный датчик случайных чисел

[maxal]

Вот еще одна библиотека Cuba - a library for multidimensional numerical integration кроме прочего предлагающая несколько разных RNG на выбор.

[Nobody]

Спасибо всем за ответы!

подынтегралная функция - это экспонента, причем в показателе экспоненты - сумма произведений логарифмов всех величин.
подскажите как сюда запостить вид функции, что то не получается.

особенности : при стремлении величин к бесконечности, функция бесконечно мала, и при стремлении величин к нулю, фунция бесконечно растет.

точность интегрирования - устроит 6 знаков после запятой, время в принципе не важно, главное чтоб результат был. пытался вычислить его в Математике, но не получается пока.

Впрочем, есть еще одна идея: а возможно ли вычислить восьмимерный интеграл по гиперобъемам Дарбу? сложить объемы и получить интеграл или нет?

Заранее спасибо.

[LynxGAV]

Никогда не поздно что-то полезное для себя выудить, зайдя в "левую" тему. Спасибо за ссылки, maxal. Какой-нибудь Максим Александрович.. Такие смелые заявки, даже не знаю.
Случайно никто не занимается методом Монте-Карло в применении к физике элементарных частиц (симуляция распадов заряженных и нейтральных частиц, термализации нейтронов, симуляция детекторов (калориметров в частности), взаимодействия излучения с веществом, симуляция столкновений частиц (пакет PYTHIA, JETSET), cчитает диаграммы Фейнмана (CompHEP))? В применении к физике конденсированного вещества (магнитные явления)? Кто-то работал с netlogo (моделирование сложных систем (в целом) с графическим интерфейсом)? С какими ЦЕРНовскими программами последних лет 10?

PS Надо будет посмотреть на это искусство. Из Вашего много чего хочется покликать.

[Николай]

Добрый день. Я сам совсем недавно (2 дня назад ) начал изучать метод Монте-Карло (понадобилось для решения некоторых диффуров из МСС), поэтому знаю и умею оччень немного. Но сегодня утром написал програмку на Фортране для вычисления примитивного одномерного интеграла от sin(x), x=[0,Pi/2]. Он равен 1, что очень удобно для оценки погрешности. Скажу, что этот интеграл я посчитал на Matlab и на Fortran. Так вот, когда я задал в Matlab 100000 "испытаний", он ушел в себя и н вернулся, пока я пил чай . А Fortran справился и с 10000000 испытаниями всего за 26 секунд (засекал секундомером ). Все основные элементарные функции в фортране есть, поэтому, полагаю, вашу функцию там можно будет задать, хотя, понятно, что обобщить мою програмку на восьмимерный случай будет непросто. Ну, надеюсь, хоть чем-то вам помогу.

Листинг:
PROGRAM INTEGRAL
REAL, ALLOCATABLE :: RANX(:), RANY(:), X(:), Y(:), S(:)
INTEGER :: NMAX
NMAX=10000000
ALLOCATE (RANX(NMAX))
ALLOCATE (RANY(NMAX))
ALLOCATE (S(NMAX))
ALLOCATE (X(NMAX))
ALLOCATE (Y(NMAX))
CALL RANDOM_SEED()
CALL RANDOM_NUMBER(RANX)
CALL RANDOM_NUMBER(RANY)
XMIN=0.
YMIN=0.
XMAX=3.14159/2
YMAX=1.5
SUM=0.
DO I = 1,NMAX
X(I)=XMIN+(XMAX-XMIN)*RANX(I)
Y(I)=YMIN+(YMAX-YMIN)*RANY(I)
IF (SIN(X(I)).GE.Y(I)) THEN
S(I)=1
ELSE
S(I)=0
ENDIF
SUM=SUM+S(I)
END DO
PL=0.
PL=(XMAX-XMIN)*(YMAX-YMIN)*SUM/NMAX
WRITE(8,3)
3 FORMAT(3X,'ПЛОЩАДЬ ФИГУРЫ')
WRITE(8,22)
22 FORMAT(' ')
WRITE(8,21) PL
21 FORMAT((2X,F12.8))
END PROGRAM INTEGRAL

Я получал результаты 1.000***. То есть, абсолютная погрешность была меньше одной тысячной. Естественно, количество испытаний (NMAX) стоит уменьшить.

[Николай]

В листинге, естественно, вместо " " было двоеточие и скобка - одномерный динамический массив

[Руст]

Встречал людей, считающих методом Монте-Карло бесконечномерные интегралы (интегралы Фейнмана по траекториям). Однако, сам отношусь к этому методу скептический. Достаточно посмотреть равномерное случайное распределение 100 точек на плоскости. На глаз видно, что есть сгустки и пустоты. В восьмимерном пространстве 10000 точек будет иметь равномерность как 3-4 точки на размерность. Хотя для достаточно хороших функций ошибки между собой компенсируются, нет гарантии, что даже для некоторых достаточно гладких нелинейных функций это имеет место.

[Николай]

Ну, если бы равномерное распределение с.в. было бы действительно равномерным - без сгустков, оно уже не было бы случайным . Когда я был школьником и ходил на курсы, к нам однажды в аудиторию по ошибке забрёл какой-то забавный товарищ (какой-то ученый, занимавшийся поиском нефти). Он привёл интересный пример распределения с.в. - нужно взять лист бумаги, сильно его смять и проткнуть иглой. Потом развернуть лист и убедиться, что дырки распределились сгустками. Это не по теме, конечно, а так - интересный пример

[Руст]

Любая случайная величина, используемая в компьютере, не является случайной. Равномерность распределения для вычисления интеграла важнее случайности. Например, если сгустки окажутся там, где функция принимает максимальное значение, а просветы, там, где у неё минимальное значение, то этот метод даст сильно завышенное значение.

[Николай]

Да, вероятно, Вы правы. Но ведь это не проблема - задать равномерное
распределение на некотором (даже N-мерном) прямоугольнике. Я немного переделал свой пример - задал 1000 точек по х и 1000 по у - т.е. миллион точек, равномерно распределённых на прямоугольнике
0<x<Pi/2
0<y<1
Результат очень неплохой 0.99957080.
Полагаю,я воспользовался каким-то общеизвестным методом?