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

Пусть всего было испытано N точек, и из них M попало в круг. Рассмотрим отношение количества точек, попавших в круг, к общему количеству точек (M/N). Очевидно, что чем больше случайных точек мы испытаем, тем это отношение будет ближе к отношению площадей четверти круга и квадрата. Таким образом, имеем, что, для достаточно больших N, верно равенство:

.                    

Из полученного равенства:   .                       (13)

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

 

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

           Program pi;

VAR

K,       {количество испытаний}

N,      {количество точек}

i, j: word;     {для циклов}

s,      {сумма всех Пи}

P: real;    {среднеарифметическое значение Пи}

{функция возвращает число  Пи}

FUNCTION raschet: real;

VAR

x, y: word;    {координаты точек}

M: word;     {число точек попавших в окружность}

BEGIN

     M:=0;

     for i:=1 to N do

begin

  x:=random(2);   {x, y – случайные числа}

  y:=random(2);

  if sqr(x)+sqr(y)<=1 then inc(M); {точка с координатами x, y попала в круг}

end;

     raschet:=4*M/N;   {из формулы [13]}

END;

BEGIN

    write('Введите количество испытаний: ');

    readln(K);

    write('Введите количество испытываемых точек: ');

    readln(N);

    randomize;

    s:=0;

    for j:=1 to K do s:=s+raschet;

    P:=s/K;

    writeln('Число Пи, рассчитанное методом Монте-Карло равно:');

    writeln(P:1:6);

    writeln;

    writeln('Точное число Пи равно:');

    writeln(Pi:1:6);

     readln;

END.

Итак, с помощью этой программы я проверил верность формулы [13]. Я получил число Пи равное: 3.000808, при количестве испытаний 500 раз с количеством точек 5000. Точное число Пи равно: 3.141593.Как и говорилось выше более точный ответ можно получить при очень большом количестве проведенных опытов, при испытании большего количества точек и при использовании качественного генератора псевдослучайных чисел.

   

2.5  Пример решения задачи методом статистического моделирования

На пластину из алюминия толщиной D = 75 см перпендикулярно ее поверхности падает пучок моноэнергетических тепловых нейтронов. Вероятность свободного пробега нейтрона в веществе до взаимодействия равна,

где μ = 0.0981 см^-1

 

В результате взаимодействия с веществом нейтрон либо поглощается  с вероятностью p_a = 0.141, либо рассеивается с вероятностью p_s = 0.859 . Считать, что задача двумерна и рассеяние нейтронов происходит изотропно в плоскости рисунка.

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

Рисунок 12 - На пластину перпендикулярно ее поверхности падает пучок нейтронов.

 

Решение

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

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

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

1. Инициализируем начальное  положение нейтрона при входе  в пластину: x = 0 , y = 0 , φ = 0 (угол  φ отсчитывается от оси x );

2. Моделируем движение  нейтрона без взаимодействия (свободный  пробег). Так как плотность распределения  длины свободного пробега нейтрона  вдоль заданного направления   φ описывается экспоненциальной  функцией распределения  , то длина пробега до взаимодействия Определяем новые координаты положения нейтрона: x = x + l cos φ , y = y + l sin φ , φ =  φ;

3. Проверяем, не вышел  ли нейтрон за границы пластины ? Если вышел ( x > D ), то вносим в гистограмму точку с углом φ  и переходим к розыгрышу траектории следующего нейтрона на п.1. Если нейтрон не вышел за границы пластины ( x ≤ D ), то переходим к следующему пункту;

4. Розыгрыш типа взаимодействия  нейтрона с веществом. Определяем, какой процесс произошел −  поглощение или рассеяние. Для этого получаем очередное случайное число γ и сравниваем его со значением вероятности поглощения нейтрона p_a . Если γ ≤ p_a , то это означает, что произошло поглощение нейтрона, поэтому процесс прослеживания текущей траектории заканчивается и переходим к п.1 для розыгрыша движения следующего нейтрона. Если γ > p_a ,то произошло рассеяние нейтрона и переходим к пункту 5;

5. Определение параметров  рассеяния нейтрона. Считая рассеяние  нейтрона изотропным, определяем  угол рассеяния по формуле  φ = 2πγ , так как рассеяние во все стороны равновероятно. Теперь, определив новое направление движения нейтрона, возвращаемся к п.2 для розыгрыша движения нейтрона в этом направлении до следующего взаимодействия;

6. Повторяем описанный  процесс необходимое число раз,  соответствующее числу нейтронов, упавших на пластину.

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

 

uses Gisto,Graph;

const

    MashX = 2; { Масштабный коэффициент по X для вывода на экран }

   MashY = 1.2; { Масштабный коэффициент по Y для вывода на экран }

    NGist = 60; { Число интервалов гистограммы }

    Mu = 0.0981; { Mu - коэффициент поглощения нейтрона }

    Pa = 0.0141; { Pa - вероятность поглощения нейтрона }

    D = 75.; { D - толщина пластины }

var

  Gd, Gm, I, N : Integer;

  Xscr, Yscr, Xscr0, Yscr0 : Integer;

  OldPattern : FillPatternType;

  x, y, Lpr, Fi : real;

label P1, P2, P3, P4, P5, OutGist;

begin

         writeln('Введите число нейтронов  N = ');

          readln(N);

{ Переходим в  графический режим }

          Gd := Detect; InitGraph(Gd, Gm, 'с:\tp\bgi');

          if GraphResult <> grOk then begin

          writeln (' Graph error or BGI driver not found. ');

          Halt;

end;

{ Рисуем пластину  в центре экрана }

   GetFillPattern(OldPattern);

   SetFillPattern(OldPattern,7); { стиль и цвет для вывода пластины }

   Bar(320-round(D*MashX/2), 0, 320+round(D*MashX/2), 479);

   Randomize; { подкрутка Random }

   i := 0; { очистка счетчика числа нейтронов }

   InitGist(NGist,-Pi/2,Pi/2); { инициализация гистограммы в диапазоне углов от -Pi/2 до Pi/2 }

 P1:

{Инициализируем начальные координаты и угол движения нейтрона }

   i := i + 1;

   if (i > N ) then goto OutGist; { проверка на исчерпание числа нейтронов }

   x := 0; y := 0; Fi := 0;

   Xscr0 := 320-round(D*MashX/2); Yscr0 := 240;

   MoveTo(Xscr0,Yscr0); { установка точки на экране в позицию x=0, y=0 }

   SetColor (i); { установка цвета вывода траектории }  

 P2:

{ Определяем  длину свободного пробега и  новые координаты нейтрона}

   Lpr := - ln(Random)/Mu; { длина пробега }

   x := x + Lpr*cos(Fi); { новая координата x }

   y := y + Lpr*sin(Fi); { новая координата y }

{ Определяем  экранные координаты и рисуем  траекторию }

   Xscr := Xscr0 + round(x*MashX);

   Yscr := Yscr0 + round(y*MashY);

   LineTo(Xscr,Yscr);

 P3:

{ Проверка выхода  за границы пластины }

   if (x < 0 ) then GoTo P1; { выход за левую границу пластины,

переходим к  следующему нейтрону }

   if (x > D) then

begin { выход за  правую границу }

   if (Fi > Pi) then Fi := Fi - 2*Pi; { коррекция диапазона по углам }

   InGist(Fi); { заносим точку в массив гистограмм }

   GoTo P1;

end;

P4:

{ Проверка на  поглощение или рассеяние нейтрона }

   if (Random <= Pa) then GoTo P1; { поглощение нейтрона }

P5:

{ Рассеяние нейтрона, розыгрыш угла рассеяния }

   Fi := 2*Pi*Random;

   GoTo P2;

   OutGist:

   write(' Вывод гистограммы <Enter> ? ');

   readln; { ожидание <Enter> }

   ClearDevice; { очистка экрана перед выводом гистограммы }

   Gistogram(13); { вывод гистограммы }

   readln; { ожидание <Enter> }

   CloseGraph;

   end.               

 

 

Рисунок 13 - Нарисован набор траекторий движения нейтронов в пластине

 

{Пакет процедур для построения гистограммы распределений}

{ Файл: gisto.pas }

UNIT Gisto;

INTERFACE

uses Graph;

var

   I, NGist : Integer; { NGist-число интервалов гистограммы }

   OldPattern : FillPatternType;

   Xstart,Xfin,deltaX : real;

   MGist : array[1..580] of integer; { массив гистограммы }

   procedure InitGist(Ng : integer; Xl,Xh : real);

   procedure InGist(Value : real);

   procedure Gistogram(Color : Byte);

IMPLEMENTATION

{ Процедура инициализации  для построения гистограммы }

{ Ng - число интервалов  гистограммы;

Xl,Xh - начальное  и конечное значение гистограммы по оси X }

procedure InitGist(Ng:integer; Xl,Xh:real);

var

i : integer;

begin

   if Ng <= 580 then Ngist:=Ng

   else begin

   writeln('Число интервалов  гистограммы > 580 !');

   Halt;

end;

   for i := 1 to NGist do MGist[i] := 0; { очистка массива гистограммы }

   Xstart := Xl; { начальное значение гистограммы по оси X }

   Xfin:=Xh; { конечное значение гистограммы по оси X }

   deltaX := (Xh-Xl) / NGist; { интервал построения гистограммы по оси X }

end;

{ Процедура занесения  значения в массив гистограммы }

procedure InGist(Value : real); { Value - вносимое значение }

var

   i : integer;

begin

   i := trunc((Value - Xstart)/deltaX)+1; { номер элемента по x}

   if i < 1 then i := 1; { защита от выхода }

   if i > NGist then i := NGist; { за границы массива }

   MGist[i] := MGist[i] + 1; { добавление точки в массив }

end;

{ Процедура вывода  гистограммы на экран }

  procedure Gistogram(Color : Byte); { Color - цвет вывода гистограммы }

const

   MaxY = 420; { максимальный размер гистограммы по Y }

   X0 = 50; Y0 = 450; { координаты начальной точки вывода }

var

   max,x,y,i,shir : integer;

   mash : real;

   s : string;

   begin

{ Определение  максимума массива и масштаба  по Y }

   max := MGist[1];

   for i := 2 to NGist do

   if MGist[i] > max then max := MGist[i];

   if max <> 0 then mash := MaxY / max;

   if max = 0 then mash:=0;

{ Определение  ширины прямоугольника для вывода  интервала }

   shir := trunc((630-X0) / Ngist);

   GetFillPattern(OldPattern); {Цвет и стиль вывода гистограммы}

   SetFillPattern(OldPattern,Color);

{ Вывод гистограммы  на экран }

   for i := 1 to NGist do begin

  y := Y0 - round(MGist[i]*mash);

  x := X0 + round(Shir*(i-1));

   bar(x+1, Y0, x + shir-1, y);

end;

{ Построение  осей и разметки гистограммы  }

SetColor(white);

line(X0-10,Y0,X0+NGist*shir+9,Y0);

line(X0,Y0+10,X0,Y0-MaxY-10);

line(X0-10,Y0-MaxY,X0+10,Y0-MaxY);

line(X0+NGist*shir,Y0+10,X0+NGist*shir,Y0-10);

 

{ Вывод значений  максимума гистограммы, начального и конечного значений по оси X }

   Str(0,s);

   MoveTo(X0-22,Y0-4);

   OutText(s);

   Str(max,s);

   MoveTo(0,Y0-MaxY-4);

   OutText(s);

   Str(Xstart:8:2,s);

   MoveTo(X0,Y0+14);

   OutText(s);

   Str(Xfin:8:2,s);

   MoveTo(X0+NGist*shir-64,Y0+14);

   OutText(s);

end;

end. 

 

Рисунок 14 - Гистограмма распределений полученная на программе Тurbo Pascal

 

Заключение о практике

 

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

  За время прохождения производственной практики:

1. получил практический опыт работы с  технической документацией;
2. ознакомился с принципами моделирования физических явлений на ЭВМ;
3. изучил теоретические основы методов статистического моделирования;
4. при помощи метода Монте-Карло и языка программирования Turbo Pascal     научился построению  некоторых трехмерных объектов, вычислению  определенного интеграла, числа π, и решения задач методом статистического моделирования;
5. пoзнакомился с рабочим коллективом и узнал о деятельности предприятия, на котором  проходил практику.

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

 

 

Список использованной литературы

 

1. Соболь И.М. «Метод Монте-Карло», М., 1985.
2. Соболь И.М. «Численные методы Монте-Карло». М.: "Наука", 1973.
3. Ермаков С. М. «Методы Монте-Карло и смежные вопросы».М.:“Наука”, 1971.
4. Фаронов В. В.. «Основы Турбо Паскаля». М., «НВТУ-Фесто Дидактик», 1992.
5. Фаронов  В. В «Turbo Pascal», – CПб.: БВХ- Петербург, 2004.
6. Ермаков, С.М. Михайлов. Г.А. «Курс статистического  моделирования». М.:"Наука", 1976.