Отчет по практике ТОО «Институте Высоких Технологий» АО «НАК «Казатомпром»

ξi+1 ={M * ξ_i} , ξ₀ = 2^{−m ,}

где M − достаточно большое целое число, фигурные скобки обозначают дробную часть, а m − число двоичных разрядов в мантиссе чисел в ЭВМ.

Методы выбора значений M , ξ0 и m разнятся для разных вариантов реализаций данного метода (это своя собственная "наука") и определяют основные свойства датчика случайных чисел (соответствие статистическим критериям, длину периода повторения последовательности и т.п.).

В составе математического  обеспечения системы Turbo Pascal также имеется датчик псевдослучайных чисел − функция Random() . Если эта функция используется без параметра, то ее значением будет равномерно распределенное реальное число в интервале от 0 до 1. Если же она вызывается с параметром (целым числом n ), то результатом будет равномерно распределенное целое число в интервале от 0 до n . Вызов функции в программе осуществляется следующим образом:

x := Random; { 0 < x < 1 }

или

m := Random(n); { 0 ≤  m < n } .

       Поскольку получение псевдослучайных чисел осуществляется по рекуррентной формуле, то при запуске программы функция Random будет всегда начинать последовательность с одного и того же первого числа (аналог ξ₀ в приведенных выше примерах). Это означает, что после компиляции и запуска программы вы все время будете получать одну и ту же наперед заданную последовательность чисел. С одной стороны, это иногда удобно для отладки программы, а, с другой стороны, может мешать при моделировании. Чтобы исключить такую ситуацию, существует специальная процедура Randomize , которая инициализирует случайным значением (текущим системным временем) генератор псевдослучайных чисел, используемый функцией Random.

       Таким  образом, обращение к процедуре Randomize позволяет смоделировать начальную "случайность" программного датчика псевдослучайных чисел.

 

         2.2 Моделирование дискретных случайных величин

Рассмотрим дискретную случайную величину ξ , принимающую n значений X₁ , X₂ , ..., X _n с вероятностями P₁ , P₂ ,..., P_n . Эта величина задается таблицей распределения

Для моделирования такой  дискретной случайной величины разбивают отрезок [0,1] на n последовательных отрезков Δ₁ , Δ₂ ,..., Δ_n , длины которых равны соответствующим вероятностям P₁ , P₂ ,..., P_n . Получают случайную величину γ , равномерно распределенную в интервале (0,1), и полагают ξ = X_k , если γ Δ_k .

Аналогичным образом  осуществляется моделирование случайных  событий. Пусть необходимо смоделировать  реализацию (осуществление или неосуществление) какого-либо события S_k из ряда событий S₁ , S₂ ,..., S_n , образующих полную группу событий. Свяжем с этой системой событий случайную величину ξ , равную номеру события:

 где P_k − вероятность наступления S_k

 

Тогда полагаем, что произошло  случайное событие S_k , если случайная величина ξ приняла значение, равное k.

 

 

 

          2.3 Моделирование непрерывных случайных величин

Пусть нам нужно получать значения случайной величины ξ, распределенной в интервале (a,b) с плотностью вероятности f (x) . Стандартный метод моделирования основан на том, что интегральная функция распределения любой непрерывной случайной величины равномерно распределена в интервале (0,1), т.е. для любой случайной величины x с плотностью распределения f (x) случайная величина       

 

(1) 

                                                       

равномерно распределена на интервале (0,1).

Тогда случайную величину ξ с произвольной плотностью распределения f (x)

можно найти следующим  образом:

1. Получаем случайную величину γ , равномерно распределенную в интервале (0,1).     

  

Такой способ получения  случайных величин называется методом обратных функций.

В качестве примера рассмотрим получение случайной величины τ , имеющей экспоненциальное распределение  с плотностью  p(τ) = λ * e^−λτ , где λ − параметр распределения. Такое распределение широко используется при моделировании различных физических явлений: это и длина свободного пробега ионизирующих частиц в веществе, и распределение интервалов времени между моментами попадания ионизирующих частиц в регистрирующий прибор и т.д.

Следуя указанной методике, получаем:

                                                        

                                                                                        (2)

Отсюда      Так как величина 1− γ распределена

точно так же, как и  γ, то последнюю формулу можно  переписать в виде

                                                                                                    (3)

 

 

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

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

                                                                                                (4)

 

Разработано много разных методов вычисления интегралов с  помощью статистического моделирования. Рассмотрим здесь лишь два наиболее простых.

Первый метод аналогичен способу определения площади  произвольной плоской фигуры, описанному выше.

 

 

Рисунок 9 - Задана некоторая функция y = f (x) по которой нужно вычислить определенный инеграл от этой функций  на интервале (a,b)

 

     Рассмотрим прямоугольник со сторонами x = a , x = b , y = 0 и y = M , где M − число, равное или превышающее по величине максимум функции на интервале (a,b) , так что 0 ≤ f (x) ≤ M на всем интервале.

Будем заполнять этот прямоугольник случайными  точками (ξ,η) , координаты которых равномерно распределены в интервалах a < ξ < b , 0 < η< M.

     Тогда геометрически очевидно, что приближенное значение интеграла (заштрихованная на рисунке площадь) будет равно площади S₀ прямоугольника, умноженной на отношение числа точек N′ , попавших под

кривую y = f (x) , к общему числу точек N . Следовательно:

 

                                                             (5)

 

На практике это означает, что для каждой разыгранной точки ( ξ_i ,η_i ), где i = 1,2,..,N , проверяется условие η_i < f (ξ_i ) , и, если оно выполнено, то в счетчик для N′ добавляется единица, в противном случае − ничего не добавляется. После проведения N испытаний по приведенной формуле вычисляется приближенное значение интеграла.

Аналогично могут быть вычислены и кратные интегралы. Например, для двойного интеграла  от функции f (x, y) по области G = {a < x < b, c < y < d } будем иметь

                                                        (6)

 

где V₀ − объем параллепипеда со сторонами (b − a) и (d − c) , и с высотой M , равной или превышающей по величине максимум функции в области G , так что во всей области                 0 ≤ f (x, y) ≤ M . В данном случае уже не прямоугольник, а параллепипед заполняется случайными точками (ξ,η,χ), координаты которых имеют равномерное распределение в интервалах a < ξ < b ,  c < η < d ,  0 < χ < M , и, если для произвольной точки ( ξi ,ηi ,χi) выполняется условие χi < f (ξi ,ηi ) , то в счетчик для N′ добавляется единица. В общем случае формула вычисления для k -кратных интегралов будет

иметь вид 

 

                                                                                (7)

 

где Vk − k-мерный объем области интегрирования.

Чем больше точек, тем  точнее значение интеграла, но следует  помнить, что точность вычисления пропорциональна 1 N , и для повышения точности в 10 раз нужно увеличить N в 100 раз и т.д.

Другой способ вычисления интегралов с помощью метода Монте-Карло  основан на вычислении среднего значения функции. Пусть ξ – случайная величина, равномерно распределенная в интервале (a,b) с плотностью вероятности

 

 

Так как любая функция f (ξ) от случайной величины ξ также будет случайной, то ее математическое ожидание равно

 

(8)

Отсюда видно, что интеграл I может быть вычислен через математическое ожидание или среднее значение случайной величины f (ξ) , являющейся функцией равномерно        распределенной случайной величины ξ . Используя оценку для среднего значения по      выборочным значениям из N испытаний, получаем

 

                                            (9)

 

Описанный способ может  также рассматриваться как модификация  метода прямоугольников для вычисления интегралов. Правда, в данном случае область интегрирования не разбивается  на множество интервалов с шагом h , которые и служат основанием прямоугольников, аппроксимирующих при суммировании значение интеграла. Наоборот, основание прямоугольников, использующихся в данном методе, равно полному интервалу интегрирования (b − a) , а вот высота этих прямоугольников случайна и равна значению функции f (x_i ) в случайной точке x_i интервала (a,b).

Будем "набрасывать" на область интегрирования такие  случайные прямоугольники, а площадь каждого прямоугольника принимать за значение интеграла на каждом этапе. Из-за случайности высоты f (x_i ) прямоугольников оценки значения интеграла будут также носить случайный характер.

 

Рисунок 10 - Показан метод “набрасывания”прямоугольников

 

Проведя N вычислений со случайными, равномерно распределенными в интервале (a,b) узлами x_i , после усреднения получим  

                                                                              (10)

 

Аналогично могут быть вычислены и кратные интегралы. Например, для двойного интеграла  формула вычислений имеет вид

                                                          (11)

где независимые случайные  величины ξ , η равномерно распределены в

интервалах a < ξ < b , c < η < d .

Для интегралов произвольной кратности k аналогично будем иметь

                                          (12)

где V_k − k -мерный объем интегрирования, а x_1i , x_2i ,,,, x_ki – случайная реализация значений независимых, равномерно распределенных в соответствующих интервалах интегрирования случайных величин

ξ₁ , ξ ₂,,,, ξ_k .

 

Найти методом Монте-Карло  значение интеграла и сравнить с  его точным значением:

 

 

Листинг программы

 

  Program mmk;

  Uses crt;

  Var

  F : double;

  N : longint; {N – число испытаний  }

  k,x,g,Integral :double ;

  i,a,b : integer;

  BEGIN

       Randomize;

       writeln('Введите промежуток интегрирования (a;b):');

       readln(a);

       readln(b);

       writeln('Введите количество случайных значений(число испытаний) N:');

       readln(N);

       k:=(b-a); {Переменной “k” присвоим значение длины промежутка интегрирования }

       writeln('k=',k);

       for i:= 1 to N do

   begin {Проведем N испытаний }

     g:=random; {g  переменная вещественного типа ,случайная величина из промежутка  [0;1]}

     x:= a + g*(b-a); {По этой формуле получается произвольная величина из  [a; b] }

     F:= F + sqrt(1-x*x);  { Подынтегральная функция }

     delay(0); {задержка,чтобы произвольные значения не повторялись }

     writeln('F=',F); {Сумма функций для N произвольных значений}

     Integral:=(1/N)*k*F ;

     writeln('Интеграл=',Integral);

   end; {Конец испытаний}

readln;

  END.     

                                            Протокол работы программы

Введите промежуток интегрирования (a;b):

0

1

Введеите количество случайных значений (Число  испытаний ) N:

1000

 

Интеграл = 7.848161045677E-0.01

 

 

Точное значение интеграла  π/4: =7.8539816339744E-0.01

 

                      

Определить  с помощью метода Монте-Карло значение числа π

Указание: Для вычисления числа π используйте определение площади четверти окружности, вписанной в квадрат со стороной 1.

Вычисление числа Пи методом Монте-Карло. Для этого рассмотрим четверть круга единичного радиуса (рис.11). Площадь круга равна , очевидно, площадь четверти круга равна:

.

Зная, что радиус круга  равен 1, получим:

 

 

 

 

 

Рисунок 11 - Определение площади четверти окружности, вписанной в квадрат со стороной 1.

 

Площадь же всего единичного квадрата OABC равна 1. Будем случайным образом выбирать точки внутри квадрата OABC. Координаты точек должны быть и . Теперь подсчитаем количество точек таких, что , т.е. те точки которые попадают внутрь круга.