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

Существует много способов вычисления числа Пи. Самым простым и понятным является численный метод Монте-Карло, суть которого сводится к простейшему перебору точек на площади. Суть расчета заключается в том, что мы берем квадрат со стороной a = 2 R, вписываем в него круг радиусом R. И начинаем наугад ставить точки внутри квадрата. Геометрически, вероятность P1 того, чтот точка попадет в круг, равна отношению площадей круга и квадрата:
P1=Sкруг / Sквадрата = πR2 / a 2 = πR2 / (2 R ) 2= πR2 / (2 R) 2 = π / 4 (1)
Выглядит это так:

image

Вероятность попадания точки в круг можно также посчитать после численного эксперимента ещё проще: посчитать количество точек, попавших в круг, и поделить их на общее количество поставленных точек:
P2=Nпопавших в круг / Nточек; (2)
Так, при большом количестве точек в численном эксперименте вероятности должны вести себя cледующим образом:
lim(Nточек→∞)⁡(P2-P1)=0; (3)
Следовательно:
π / 4 = Nпопавших в круг / Nточек; (4)
π =4 Nпопавших в круг / Nточек; (5)
НО! При моделировании мы применяем псевдослучайные числа, которые не являются случайным процессом.
Поэтому, выражение (5), к сожалению, строго не выполняется. Логичны вопросы, каковы оптимальные размеры квадрата и как много нужно применить точек?
Чтобы это выяснить, я написал такую программу:
#include<stdio.h> 
#include<math.h> 

#define limit_Nmax 1e7 //Максимальное количество точек
#define limit_a 1e6 //Максиальный радиус круга
#define min_a 100 //Начальный радиус

double circle(double, double); //Выдает квадрат y в зависимости от координаты Х и радиуса круга.

int main() 
{ 

double x,y,Pi; 
long long int a=min_a//сторона квадарата
i=0;//Счетчик 
double Ncirc=0;//Количество точек, попавших в круг 
double Nmax=a; //Общее количество точек
while (a<limit_a)  //Перебор  значений радиуса
{ 
Nmax=a; 
 while (Nmax<=limit_Nmax) // Перебор значения количества точек
 { 
 Ncirc=0; i=0; //обнуляторы
    while (i<Nmax) 
    { 
    x = (random() % (a * 1000))/1000;  //Рандомный Х с 3 знаками после запятой
    y = (random() % (a * 1000))/1000;  //Рандомный Y с 3 знаками после запятой
        if (y*y<=circle(x,(a/2)))  //Условие принадлежности точки к кругу
        { 
        Ncirc++; 
        } 
    i++; 
    } 

 Pi=(Ncirc/Nmax)*4; 
 Nmax *= 2; 

 printf("\n%lld,%.0f,%f",a,Nmax,Pi); 
 } 
a*=2; 
} 

} 

double circle(double x, double radius) 
{ 
double y=radius*radius-x*x; 
return y; 
}

Программа выводит значения числа Пи в зависимости от радиуса и количества точек. Единственное, что остается читателю, это скомпилировать её самостоятельно и запустить с параметрами, которые желает он.

Приведу лишь одну таблицу с полученными значениями:
Радиус
Nточек
Pi
102400
204800
3,145664
102400
409600
3,137188
102400
819200
3,139326
102400
1638400
3,144478
102400
3276800
3,139875
102400
6553600
3,142611
102400
13107200
3,140872
102400
26214400
3,141644
102400
52428800
3,141217
102400
1,05E+08
3,141324
102400
2,1E+08
3,141615
102400
4,19E+08
3,141665
102400
8,39E+08
3,141724
102400
1,68E+09
3,141682


Если что, значение числа Пи можно посмотреть с точностью до определенного знака здесь.
Источник картинки — википедия.
Поделиться публикацией
Комментарии 37
    +6
    Вопрос, зачем округлять рандомные значения до 3 знаков после запятой? Сколько точек не бери, в пределе вычислится не π, а площадь некоей многоугольной фигуры.
      0
      Хотя я засомневался что-то, наверное всё же будет стремиться к π. Но какой смысл в округлении?
      –2
      Это не округление. Дело в том, что рандомизатор выдает только целые числа. Поэтому, чтобы избежать только целых значений Х и Y, я сделал 3 знака после запятой. Можно и больше знаков сделать. Просто не вижу в этом смысла пока.
        +9
        «Округлять» нельзя — в этом то и смысл Монте-Карло. С помощью округления точки могут попасть лишь на пересечении линий, которые отстают друг от друга на растояние dx — которое в Вашем случае равно 0.001. А теперь представьте себе вырожденный случай, когда в квадрате не ~100000 таких пересечений, а например лишь 4 либо 9 — и станет понятно к чему приводит округление… Чем больше знаков после запятой — тем точнее результат — при условии, что генератор действительно случайный:)
          –2
          Мы с Вами знаем, что любая программа использует совсем даже не случайный генератор. Поэтому я хочу заняться этой темой. Не посоветуете какой-нибудь хороший генератор псевдослучайных чисел?
            +3
            34
            +9
            4
              0
              Если вы работаете под линуксом, посмотрите на разницу между /dev/random и /dev/urandom, в любом случае потребуется «прогрев».

              В Windows XP генератор псевдослучайный чисел мягко говоря безобразный, в Висте и 7ке лучше. Почитайте про CryptoAPI и, в частности, CryptGenRandom.
                0
                Так вот зачем Intel встраивает аппаратный генератор случайных чисел! Чтобы skobeltsyn смог точно вычислить площадь многоугольника близкого к дуге окружности. На самом деле подлинность случайности при выборе точек в данном алгоритме не важна.
                  0
                  А что важно?
                    +2
                    разрядность действительных чисел, количество итераций
                  0
                  есть хороший генератор на основе LPtau последовательностей
                    0
                    Спасибо! Я обязательно его опробую!
              +2
              Не хватает колонок с погрешностями.
                –2
                А что бы дали эти погрешности?
                –5
                Согласен, но чем больше я беру радиус, тем меньше будет погрешность от округления. Заковыка в поиске оптимального размера фигур и количества точек. То есть, когда я беру радиус 100 точек, то число пи получается порядка 3.18… не слишком удачный результат. Но чем больше круг(от 1000 и более), тем точнее становится расчёт. Минус метода заключается в том, что я не контролирую точность из-за того, что генератор псевдослучайный.
                  0
                  это с какой это стати?

                  относительные погрешности одинаковы, вне зависимости от радиуса. Вы же только умножаете или делите на него
                    0
                    … а поскольку абсолютное значение величины фиксировано, то, стало быть, из одинаковости относительных вытекает и одинаковость абсолютных
                      0
                      У меня ограничена разрядность. Из-за этого, чем больший радиус я покрою, тем большее количество разнообразных точек получу, тем больше процесс будет похож на случайный. И тем большая точность.
                  +9
                  это задача школьной программы. Зачем писать об этом на Хабре?
                    –5
                    Задачка интересная. Почему бы её и не вспомнить?
                      +4
                      Я полагаю, что вы бы могли рассказать не столько и не только о вычислении числа Пи, сколько о применении метода Монте-Карло для решения различных геометрических задач.
                        –4
                        Я только начал заниматься этим вопросом. В перспективе вычисление этим методом диффузии веществ в микропористых средах. Решение геометрических задач методом Монте-Карло — это очень известная тема, о которой, считаю, можно не упоминать.
                          +6
                          Решение геометрических задач методом Монте-Карло — это очень известная тема, о которой, считаю, можно не упоминать

                          А вот вычисление числа Пи — это тайные знания масонов, да?

                          Решение геометрических задач — это более-менее реальный способ применения метода Монте-Карло, в отличие от вычисления числа Пи (которое, по сути, лишь является примитивной иллюстрацией самого алгоритма).
                            –2
                            Это и есть решение геометрической задачаи. Если разобраться повнимательнее в алгоритме вычисления пи, то можно сразу заметить, что тут вычисляется площадь круга методом Монте-Карло. И по идее она же должна делиться на площадь квадрата. Так что косяк в алгоритме детектед. И он и не в разрядности, и не в генераторе.
                              –1
                              Простите, возможно я не совсем точно выразился. Я имел в виду то, что одно дело — «научить» хабраюзеров находить площадь круга, а другое — рассказать, как же методом Монте-Карло вычислять, например, многомерные интегралы или находить корни системы линейных алгебраических уравнений.

                              Вас же минусуют не только за, как вы выразились, «косяк в алгоритме», но и за то, что вы нас недооцениваете :).
                                +1
                                И в чём же это проявляется?
                    +8
                    i=0;//Счетчик 
                    
                    Ncirc=0; i=0; //обнуляторы
                    

                    Сделайте меня развидеть это :(
                      –1
                      А что тут такого?
                        0
                        Слово то какое… заморское!
                          0
                          нулифайеры!
                            0
                            Ncirc=(Ncirc-Ncirc)*0; //Брутальный обнулятор
                      +3
                      И всё ради этого?
                      #include <boost/random/normal_distribution.hpp>
                      #include <boost/random/mersenne_twister.hpp>
                      
                      bool inCircle(double x, double y) {
                      	return x*x + y*y < 1;
                      }
                      
                      double secretFunction(size_t n) {
                          boost::uniform_01<boost::mt19937> gen((boost::mt19937()));
                          size_t inside = 0;
                        
                          for (size_t i = 0; i < n; i++)
                              inside += inCircle(gen(), gen());
                          
                          return 4.0 * inside / n;
                      }

                      Заодно приделал вихрь Мерсенна, который для метода Монте-Карло и создавался.
                        +1
                        Скажите, а какая точность у вашей реализации алгоритма?
                          +1
                          Как уже было сказано выше, никого не интересует таблица с полученными значениями, ибо рандом по определению. Я прогнал по 20 раз для каждого значения n от 1 до 20 миллионов (итого 400 замеров).

                          Наверняка я что-то накосячил со статистикой, но вот что получилось (все графики в зависимости от числа миллионов случайных точек):

                          Синяя линяя — стандартное отклонение.
                          Желтая линия — стандартное отклонение / pi.
                          Красная линия — предполагаемое отклонение. Вот тут я мог напутать. По моим подсчётам это (4/pi-1)/sqrt(n). По другим источникам — 1/sqrt(n). Что общего в формулах — так это то, что погрешность от double-арифметики при N < 2^32 всегда будет всегда намного меньше погрешности метода Монте-Карло.
                            +2
                            Погрешность стремится приблизительно к 0.0001. Это значит, что 4й знак после запятой — не верный. Этот же вывод можно сделать из таблички в статье. Причём это заметно уже и при 5млн замерах. Метод расчета погрешности мне понравился. Взял на вооружение.

                      Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                      Самое читаемое