Программы, использующие случайные числа и рекурсию Задача 3.

Датчики случайных чисел находят интересное применение в вычислении пло­щадей и объемов: на этом применении основан известный метод Монте-Карло. Требуется, воспользовавшись датчиком случайных чисел, определить площадь фигуры.

Решение.

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

Такой метод называют методом статистических испытаний или методом Монте-Карло. Свое название он получил от названия небольшого городка, знаменитого игорными домами, в которых применяется рулетка - про­стейший механический прибор для получения случайных чисел. Создателями ме­тода считаются американские математики Дж. Нейман и С. Улам, а годом созда­ния - 1949 г. Хотя теоретическая основа метода известна давно, только с появле­нием ЭВМ этот метод нашел широкое применение и стал универсальным числен­ным методом.

Program Examp_29;

Uses crt, graph;

Var x,y:Integer;

ps,pf:LongInt;

z:String;

Procedure Graphinterface;

Var

driver, mode, error:Integer;

s:String;

Begin

driver:=detect;

s:='';

Initgraph(driver,mode,s);

error:=GraphResult;

if error<>GrOk then

begin

writeln(GraphErrorMsg(Error));

Halt(error)

end

end;

begin

Graphinterface;

SetTextStyle(DefaultFont,HorizDir,3);

OutTextXY(100,50,'Метод Монте-Карло');

SetTextStyle(DefaultFont,HorizDir,1);

OutTextXY(320,220,'Количество точек,

принадлежащих');

OutTextXY(320,220,'фигуре ');

OutTextXY(320,265,'квадрату :');

{оси координат и квадрат 200х200 точек}

Line(50,100,50,370);

Line(50,100,45,105);Line(50,100,55,105);

Line(300,350,30,350);

Line(300,350,295,345);Line(300,350,295,355);

OutTextXY(245,360,'1');OutTextXY(35,145,'1');

Rectangle(50,150,250,350);

{фигура: 1)эллипс; 2) прямоугольник; 3) круг}

SetColor(9); SetFillStyle(1,1);

{1} FillEllipse(140,240,40,80);

{2} {Bar(150,150,250,350);}

{3} {FillEllips(150,250,100,100);}

Randomize;ps:=0;pf:=0;

Repeat

x:=50+Random(201); {[50,250]}

y:=150+Random(201) {[150,350]}

SetViewPort(320,240,320,240,ClipOff);

Repeat

putpixel(round(r*x(t)),round(r*y(t)),11);

t:=t+0.001;

Until t>2*pi;

Readln; CloseGraph;

End.

Найдем, например, площадь эллипса, применяя метод Монте-Карло. Сторона) квадрата, условная единица измерения площади, в программе равна 200 пикселей. Используем формулу а+Random(b-а+1) для получения случайных координат точек, принадлежащих квадрату: х[50, 250], у[150, 350]. С помощью функции GetPixel(x,y:Integer), возвращающей код цвета пикселя с координатами (х,y), определим принадлежность точки либо только квадрату, либо одновременно квадрату и эллипсу (для этого точки внутри квадрата и внутри фигуры маркируются различными цветами). Ведя при этом раздельный подсчет точек (переменные рf и рs), найдем искомую площадь как отношение полученных чисел.

В качестве контрольного примера можно рассмотреть конфигурацию квадрата и вписанного в него круга. В этом случае известен точный ответ, равный /4, так как r²/(4r²)=/4. Если же полученный при этом приближенный ответ умножит на 4, то мы найдем еще одним способом приближенное значение константы n.

Задача 4. Доска Гальтона.

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

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

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

Решение.

Пусть т = 5 и введен коэффициент пропорциональности А. Тогда общее коли­чество шариков равно 16А, так как при т = 5 сумма чисел в соответствующей строке треугольника Паскаля равна 16:

Выберем систему координат с началом в точке (320,10) и введем условную единицу длины d. Для рисо­вания доски служит процедура Board. Вложенная про­цедура Half учитывает симметрию доски: изображает одну половину контура доски и вызывается дважды с противоположными по знаку параметрами. Другая внут­ренняя процедура Obstacle создана для рисования одно­го препятствия. Все 10 препятствий рисуются с помо­щью процедуры во вложенных циклах. Здесь же рису­ются и четыре вертикальных отрезка, которые разделя­ют отделения доски.

Процедуры Ball и PutBall создают и перемещают гра­фический образ (шарик) по доске. Они используют про­цедуры GetImage и PutImage модуля Graph, записываю­щие изображение в буфер и восстанавливающие его из буфера.

Процедуры Upper и Lower вызываются многократно и реализуют прохождение шариков по доске. В верхней части доски (от 0 до 16d) шарик четыре раза (2d, 6d, 10d, 14d) выбирает путь налево или направо: с помощью датчика случайных чисел в процедуре Choise опре­деляется приращение абсциссы. В нижней части доски (от 16d до 36d) шарики распределяются по отделениям и подсчитываются в каждом из них (соответствую­щие количества - элементы массива Q). Они укладываются по пять штук в ряд, что вместе с радиусом (константа r) используется в формуле определения текущей абсциссы:

x=x-d+(Q(num) mod 5)*(3r+1).

Program Examp_30;

Uses crt, graph;

Const d=12; {13}

r=2;

del=#219#219#219;

Var i,j,x,y,k,n:Integer;

Q:Array[1..5] of Integer;

b:Pointer;

z:String;

Procedure Graphinterface;

Var

driver, mode, error:Integer;

s:String;

Begin

driver:=detect;

s:=' ';

Initgraph(driver,mode,s);

error:=GraphResult;

if error<>GrOk then

begin

writeln(GraphErrorMsg(Error));

Halt(error)

end

end;

Procedure Init;

begin

ClrScr;

SetViewPort(GetMaxX div 2,10,

GetMaxX,GetMaxY,ClipOff);

SetColor(2);

SetLineStyle(SolidLn,0,ThickWidth);

Line(-6,0,-25*d,0); Line(6,0,25*d,0);

OutTextXY(3*d,-10,'(осталось)');

SetColor(10);

OutTextXY(-10*d-180,38*d,'Числа Паскаля ->');

SetTextJustify(CenterText, TopText);

Str(6*k,z); OutTextXY(0,38*d,z);

Str(4*k,z); OutTextXY(4*d,38*d,z);

OutTextXY(-4*d,38*d,z);

Str(k,z);OutTextXY(8*d,38*d,z);

OutTextXY(-8*d,38*d,z);

SetFillStyle(9,2);

Randomize;

for i:=1 to 5 do Q[i]:=0;

n:=0;

end;

Procedure Board;

Procedure Half(a:Integer);

begin

MoveTo(a,0);

for i:=1 to 4 do

begin

LineRel(0,abs(a));LineRel(a,abs(a))

end;

LineRel(0,abs(10*a)); LineRel(-5*a,0);

end;

Procedure Obstacle(x,y:Integer);

begin

MoveTo(x,y);

LineRel(d,-d); LineRel(0,-2*d); LineRel(-d,-d);

LineRel(-d,d); LineRel(0,2*d); LineRel(d,d);

FloodFill(x,y-d,10);

end;

begin

Half(2*d); {правая половина}

Half(-2*d); {левая половина}

for i:=1 to 4 do

for j:=1 to i do

begin

Obstacle(2*(2*i-j-4)*d,(23-4*j)*d);

{препятствия}

if j=1 then LineRel(0,17*d)

{перегородки}

end;

end; {board}

Procedure Ball;

var size:Word;

begin

SetColor(13); SetFillStyle(1,13);

FillEllipse(3,3,r,r);

size:=ImageSize(0,0,6,6);

GetMem(b,size); GetImage(0,0,6,6,b^);

SetColor(0);FillEllipse(3,3,r,r);

FreeMem(b,size);

end; {Ball}

Procedure PutBall;

const v=5;

{константа, влияющая на скорость движения

шарика}

begin

PutImage(x,y,B^,XorPut); Delay(v);

PutImage(x,y,B^,XorPut);

Inc(y);

end; {PutBall}

Procedure Upper;

{прохождение верхней части доски}

var dx:-1..1;

Procedure Choise; {выбор пути}

begin

dx:=-1+2*Random(2);

{приращение абсциссы:1 или -1 }

Sound(40);Delay(1); NoSound;

end; {Choise}

begin

Inc(n); Str(n,z);

SetColor(0); OutTextXY(0,-10,del);

SetColor(13); OutTextXY(0,-10,z);

Str(16*k-n,z);

SetColor(0); OutTextXY(12*d,-10,del);

SetColor(2); OutTextXY(12*d,-10,z);

x:=0; y:=0;

for j:=1 to 4 do

begin

for i:=1 to 2*d do PutBall; {вниз}

Choise; {налево или направо?}

for i:=1 to 2*d do

begin

PutBall; Inc(x,dx)

end; {наискось}

end;

end; {Upper}

Procedure Lower;

{прохождение нижней части доски}

const botton=36*d-2*r-2;

{ордината центров шаров нижнего ряда}

var top:Word;

num:1..5;

begin

num:=x div(4*d)+3;

{определение номера отделения}

top:=botton-(Q[num] div 5)*3*r;

if top<19*d then top:=19*d;

Repeat

Putball; {Inc(y)}

Until y=top;

Sound(3000); Delay(2); NoSound;

SetColor(0);OutTextXY(x,37*d,del);

Str(Q[num]+1,z);

SetColor(13);OutTextXY(x,37*d,z);

x:=x-d+(Q[num] mod 5)*(3*r+1);

{абсцисса фиксируемого}

FillEllipse(x,y,r,r); {в отделении шарика}

Inc(Q[num]) {подсчет шариков в отделениях}

end;

begin

ClrScr; Writeln('Всего шариков:16*k');

Write('Введите коэффициент (от 1 до 60)

k = ?'#8);

Readln(k);

Graphinterface;

Init; Board; Ball;

Repeat Upper;

Lower Until (n=16*k) or KeyPressed;

Write(#7);

Readln; CloseGraph;

End.