Компьютерра - Компьютерра PDA N71 (06.11.2010-13.11.2010)

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

Поэтому большинство расчетов по методу Монте-Карло производится с использованием генераторов псевдослучайных чисел.

- Как устроены генераторы псевдослучайных чисел?

- Большинство таких генераторов основаны на применении так называемого метода вычетов и его модификаций. Идея довольно проста. Берется дробное число α i с большим "хвостом" (то есть с длинной мантиссой), умножается на большое целое число M , в результате получается большое целое плюс дробная часть. Потом целую часть результата убирают, а дробную берут в качестве следующего числа:

α i +1= { M α i }

Оказывается, если множитель M взять достаточно большим (например, в современных генераторах используются множители порядка M = 5 100109) получается, что "хвосты" α i +1ведут себя как настоящие стандартные случайные числа α.

На самом деле "настоящее" (теоретическое) значение стандартного случайного числа получить невозможно, так как α представляет собой дробь с бесконечной мантиссой , состоящей из нулей и единиц (такую дробь в принципе воспроизвести нельзя). Здесь ситуация похожа на проблему воспроизведения вещественных (в частности, иррациональных) чисел на компьютере.

На практике в методе вычетов при представлении чисел α i берут "длинные" мантиссы (например, в современных генераторах используется T = 128 разрядов мантиссы).

В методе вычетов имеется также проблема периодичности : не позднее, чем через 2 T шагов произойдет "зацикливание" генератора. В расчетах по методу Монте-Карло не рекомендуется использование более чем L /2 обращений к генератору; здесь L - длина периода, равная числу шагов метода вычетов, после которого начинается повторение последовательности α i . При удачном подборе множителя M можно получить величину периода L = 2 T - 2(это едва ли не "рекордный" результат). Для T = 128 величина L /2 равна 2 125, этого вполне хватает для широкого класса современных задач, решаемых с помощью численного статистического моделирования.

Решение проблем конечности мантиссы (периодичности) не гарантирует качества получаемых чисел α i . Требуется проведение тестов, показывающих, что эти числа по свойствам близки к настоящим (теоретическим) стандартным случайным числам α (тем, что имеют бесконечную мантиссу). Здесь используют широкий спектр критериев и методик математической статистики.

Тестом можно считать и любую задачу с известным ответом, решаемую методом Монте-Карло. В этом смысле процесс проверки генераторов псевдослучайных чисел неограничен. Более того, для любого генератора, основанного на методе вычетов, можно найти "тяжелую" задачу, с которой он "не справится" (то есть правильный ответ не получится). Для такой задачи придется проводить усовершенствование метода вычетов. В частности, можно увеличить длину контролируемой мантиссы T и множитель M .

Следует, однако, учитывать, что увеличение этих величин ведет к росту компьютерных затрат при обращении к подпрограммам типа RAND и RANDOM (если в этих подпрограммах "запаян" метод вычетов). Вообще следует отметить, что обращение к генератору случайных чисел - достаточно дорогостоящая компьютерная операция (по сравнению, например, с простым сложением или умножением чисел). Поэтому считается, что тот алгоритм метода Монте-Карло будет работать эффективнее (быстрее), который использует меньше обращений к генератору псевдослучайных чисел.

- А какие задачи решаются методом Монте-Карло?

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

Излучение трактовалось как поток малых частиц ("фотонов"), пролетающих сквозь слой защиты, в котором хаотично расположены крупные частицы. "Фотон", сталкиваясь с крупной частицей, либо захватывается ею ("поглощается"), либо рассеивается по некоторому вероятностному закону. Можно проследить (реализовав на компьютере) траектории фотонов и подсчитать, какая доля фотонов поглощается в слое защиты. Если эта доля близка к единице, то защита может считаться хорошей.

При реализации траектории "фотона" до поглощения нужны выборочные значения случайных величин с различными законами распределения. Для получения таких значений используют преобразования стандартных случайных чисел α j .

Далее выяснилось, что с описанным случайным процессом движения "фотонов" можно соотнести определенное уравнение (интегральное уравнение Фредгольма второго рода), на основе которого можно строить так называемые весовые оценки для вычисления требуемых характеристик (функционалов) физического процесса. Введение весов позволяет в ряде случаев упростить компьютерную реализацию траекторий "фотонов".

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

Вообще в литературе метод Монте-Карло обычно представляется как специальный способ вычисления многократных интегралов. Часто для иллюстрации рисуют такую картинку.

Численное интегрирование функции методом Монте-Карло (график из "Википедии")

Предположим, нам нужно вычислить интеграл, равный площади S под кривой, изображенной на рисунке. Для этого поместим ее в прямоугольник с известной площадью U , и будем кидать в него равномерно распределенные случайные точки. Понятно, что вероятность P попадания случайной точки в интересующую нас область равна отношению площади этой области к площади прямоугольника: P = S / U . Реализуем большое количество точек N , и подсчитаем, какое количество точек K попадет под кривую. Частота K / N попадания случайных точек под кривую приближает вероятность P , и поэтому S / U ≈ K / N , а искомый интеграл приближенно равен S ≈ KU / N .