Материал: Имитационное и статистическое моделирование — копия
лонке реализуется случайное направление вектора перемещения частицы, во второй – абсцисса конца вектора (имеющего единичную длину), в третьей – его ордината. Последняя строка таблички – A3:C3 – формирует следующее положение частицы (в момент времени t = 1), определяемое расчётом смещения частицы относительно предыдущего положения (сформированного в ячейках B2 и C2). Эту строку A3:C3 можно копировать вниз. Скопируйте её до 1001-й строки листа. В результате в колонках B и C появляется последовательность из 1000 положений (координат) частицы, совершающей Броуновское движение на плоскости.
Выделите диапазон B2:C1001 и постройте по нему диаграмму типа Точечная вида Диаграмма со значениями, соединёнными отрезками без маркеров. Реализуйте несколько независимых Броуновских движений посредством нажатия клавиши F9. На рис. 1.5 показаны три реализации Броуновского процесса, полученные описанным способом.
|
|
40 |
|
|
|
|
40 |
|
|
|
|
40 |
|
|
|
|
30 |
|
|
|
|
30 |
|
|
|
|
30 |
|
|
|
|
20 |
|
|
|
|
20 |
|
|
|
|
20 |
|
|
|
|
10 |
|
|
|
|
10 |
|
|
|
|
10 |
|
|
|
|
0 |
|
|
|
|
0 |
|
|
|
|
0 |
|
|
-40 |
-20 |
0 |
20 |
40 |
-40 |
-20 |
0 |
20 |
40 |
-40 |
-20 |
0 |
20 |
40 |
|
|
-10 |
|
|
|
|
-10 |
|
|
|
|
-10 |
|
|
|
|
-20 |
|
|
|
|
-20 |
|
|
|
|
-20 |
|
|
|
|
-30 |
|
|
|
|
-30 |
|
|
|
|
-30 |
|
|
|
|
-40 |
|
|
|
|
-40 |
|
|
|
|
-40 |
|
|
|
|
Рис. 1.5. Реализации Броуновского движения на плоскости |
|
|
||||||||||
Оцените на глаз расстояние, на которое частица удаляется от своего начального положения за 1000 шагов. Теоретически оно в среднем должно быть примерно равно квадратному корню из числа шагов (т. е. составлять около 32 единиц).
3.3.2. Броуновское движение в замкнутой области
При исследовании Броуновского движения, в котором нужно учитывать какие-либо ограничивающие факторы, аналитические методы математики оказываются малоэффективными. Статистическое моделирование позволяет учитывать подобные факторы без особого труда.
Пример 1 0 . На рис. 1.6 приведена траектория Броуновского движения
вкольце, полученного путём небольшой модификации модели, рассмотренной
вп. 3.2.1. Справа на рис. 1.6 показан увеличенный фрагмент этой же траектории частицы. Модель Броуновского движения в кольце отличается от модели Броуновского движения на неограниченной плоскости только тем, что посредством условного оператора =ЕСЛИ… выполняется проверка, не выходит ли новое случайное положение частицы за пределы кольца, и если выходит, то расчётные координаты определяются как координаты «отражения» частицы от достигнутой границы кольца. Выполните самостоятельно указанную несложную модификацию модели Броуновского движения на плоскости.
21
Рис. 1.6. Реализация Броуновского движения в кольце
При необходимости нетрудно ввести в модель и другие дополнительные факторы – мембраны и фильтры, влияющие на перемещение частицы, анизотропию направлений, удлиняющую смещение от центра по сравнению со смещением к центру (этим способом можно промоделировать вращающиеся цилиндры сепаратора) и т. д.
3.3.3. Фрактальные пейзажи
В математическом эссе [25] приводится обзор результатов относительно новой математической теории, основанной на теории множеств и теории размерностей, которая позволяет простыми средствами описывать и исследовать большое число пространственных (геометрических) закономерностей структурирования живой и неживой материи. В основе этой теории лежит понятие фрактала – множества точек, обладающего свойством наглядного или абстрактного (рис. 1.7) «самоподобия» (масштабной инвариантности) и характеризуемого наряду с обычной целой (Евклидовой) размерностью дробной размерностью Хаусдорфа – Безиковича.
Броуновское движение представляет собой одну из разновидностей фракталов.
Рис. 1.7. Случайные пейзажи (графики построены в Ms Excel)
22
На рис. 1.7 приведены две диаграммы, на каждой из которых верхний график (лес) представляет собой зависимость абсолютной величины ординаты Броуновской частицы на плоскости от времени (номера шага). Нижний график (отражение в озере) получается умножением верхнего графика на (–0,5).
Постройте такие графики самостоятельно, это несложно. Интересной фрактальной особенностью таких графиков является то, что при растягивании рисунка с графиком по горизонтали в любое число раз остаются неизменными все статистические характеристики «кромки леса», воспринимаемые визуально. Этот эффект вызывает впечатление, что происходит не растягивание рисунка, а многократное приближение к реальному горизонту пейзажа. При растягивании рисунка между верхушками деревьев появляются всё новые и новые верхушки. Указанное фрактальное свойство графика проявляется уже при использовании приближения, построенного на одной – десяти тысячах перемещений блуждающей частицы (точное математическое Броуновское движение складывается из бесконечного множества бесконечно малых перемещений). Нажимая клавишу F9, можно за несколько секунд сгенерировать десятки таких фрактальных пейзажей.
4.Расчёт интегралов методом Монте-Карло
4.1.Основная идея
Основные вероятностные характеристики непрерывной с. в. x определяются интегралами, связанными с её п. р. в. f(t). К их числу относятся следующие характеристики:
– вероятность попадания x в интервал (a,b)
P{a ≤ x ≤b} = ∫b |
f (t)dt ; |
(1.31) |
a |
|
|
– математическое ожидание x |
|
|
M (x) = ∞∫t f (t)dt ; |
(1.32) |
|
−∞ |
|
|
– математическое ожидание любой функции φ от с. в. x
M[ϕ(x)] = |
∞∫ϕ (t) f (t)dt |
(1.33) |
|
−∞ |
|
и, в частности, математическое ожидание xr (r-й начальный момент)
M (xr ) = |
∞∫t r f (t)dt . |
(1.34) |
|
−∞ |
|
Часто эти интегралы сложны для вычисления. При статистическом моделировании вместо точных значений (1.31)–(1.34) можно вычислять и использовать соответствующие приближённые эмпирические оценки, легко рассчитываемые по выборке x1, x2, …, xn:
23
– оценку вероятности попадания x в интервал (a,b) |
|
|||||
ˆ |
|
1 |
n |
|
||
P{a ≤ x |
≤b} = |
|
n |
I [a ≤ xi ≤b] , |
(1.35) |
|
|
|
|
∑i=1 |
|
||
где функция I [A] – индикатор события A – равна 1, когда событие A имеет ме- |
||||||
сто, и равна нулю в противном случае; |
|
|
|
|
|
|
– оценку математического ожидания x |
|
|||||
ˆ |
|
1 |
|
n |
|
|
M (x) = |
n |
xi ; |
(1.36) |
|||
|
|
∑i=1 |
|
|||
– оценку математического ожидания любой функции φ(x) от с. в. x
ˆ |
1 |
n |
|
|
M [ϕ(x)] = |
|
ϕ (xi ) , |
(1.37) |
|
n |
||||
|
∑i=1 |
|
и, в частности, оценку математического ожидания xr (r-го начального момента)
ˆ |
r |
) = |
1 |
n |
r |
(1.38) |
M (x |
|
n |
∑i=1 |
xi . |
||
|
|
|
|
|
Оценки (1.35)–(1.38) легко вычисляются, и погрешности, с которыми они приближают точные значения характеристик, легко контролируются. Действительно, как суммы большого числа независимых с. в. они имеют нормальное распределение. При этом м. о. оценок характеристик совпадают с точными значениями характеристик, а фактические отклонения оценок от точных значений характеристик подчиняются известному правилу «трёх сигм» [16] и могут быть сделаны сколь угодно малыми путём увеличения объёма выборки n.
Более подробные сведения о статистических оценках, включая теорию состоятельности, несмещённости и эффективности оценок, можно найти в учебниках по математической статистике [23, 24].
Основная идея метода Монте-Карло состоит в том, чтобы использовать расчёт эмпирических оценок в качестве способа вычисления определённых интегралов вида (1.31)–(1.34). Поскольку к такому виду легко сводятся практически любые определённые интегралы, метод Монте-Карло оказывается достаточно универсальным методом.
Преимущества метода Монте-Карло и условия, в которых они проявляются, можно увидеть, рассмотрев некоторые примеры его применения.
4.2. Сведение интеграла к вероятности
Наиболее простой способ вычисления определённого интеграла методом Монте-Карло – это сведение интеграла к вероятности. Суть этого способа легко понять из следующего примера.
Пример 1 1 . Рассчитаем методом Монте-Карло определённый интеграл
S = ∫1 |
1−t 2 dt , |
(1.39) |
0 |
|
|
24
который представляет собой площадь чет- |
|
y(t) = |
1 − t 2 |
|||
верти круга единичного радиуса с цент-ром в |
1 |
|||||
начале координат (рис. 1.8). Если в еди- |
|
|
||||
|
|
|
||||
ничный квадрат, ограниченный на рис. 1.8 |
|
|
|
|||
координатными осями и штриховой линией, |
|
S |
|
|||
бросить случайную точку, то она попадёт |
|
|
||||
в круг с вероятностью p = S/S0, где S0 – пло- |
|
|
t |
|||
щадь квадрата – известна. Отсюда следует, |
|
|
||||
|
|
|
||||
что |
искомое |
значение |
интеграла |
|
0 |
1 |
S = pS0 однозначно выражается через веро- |
Рис. 1.8. К расчету интеграла |
|||||
ятность p. В данном случае S0 = 1, значит |
|
|
|
|||
S = p. Вероятность p легко можно оценить |
|
|
|
|||
методом Монте-Карло, выбрасывая в квадрат достаточно большое число n слу-
ˆ |
|
|
|||
чайных точек. Оценка p вероятности p вычисляется при этом по формуле |
|
||||
|
1 |
n |
|
|
|
pˆ = |
∑i=1 |
I (Ai ) , |
(1.40) |
||
n |
|||||
где событие Ai – это событие попадания точки в круг.
Случайные точки должны иметь равномерное распределение внутри квадрата. Такие точки мы уже «выбрасывали» в примере 4 (см. рис. 1.2). Их координаты записаны в файле «СлВеличины.xls» на листе Выборка в диапазоне A1:B10000. Откройте этот файл или заполните на любом новом листе Ms Excel указанный диапазон значениями БСВ с помощью функции =СЛЧИС(). Чтобы вычислить оценку (1.40), введите в ячейку С1 формулу =ЕСЛИ((A1^2+B1^2<1);1;0). Эта формула выдаёт значение 1, если сумма квадратов координат случайной точки не превосходит 1, и выдаёт 0 – в противном случае. Таким образом, введённая формула вычисляет индикатор попадания точки в круг. Скопируйте её в диапазон C1:C10000 двойным щелчком по маркеру заполнения выделенной ячейки C1. В ячейку D1 введите формулу =СУММ(C1:C10000)/10000, вычисляющую оценку (1.40). В этой ячейке завершается расчёт приближённого значения искомого определённого интеграла S методом Монте-Карло.
В ячейку E1 введите точное значение величины S, составляющее, очевидно, S = π/4 = 0,785398, в ячейку F1 – формулу =D1-E1, вычисляющую фактическое отклонение приближённой оценки от точной. После завершения пересчёта листа несколько раз нажмите функциональный ключ F9 и понаблюдайте за величиной погрешности в ячейке F1.
Убедитесь, что погрешность не выходит за пределы, определённые правилом «трёх сигм». Согласно этому правилу, с вероятностью 1–0,0027 погреш-
|
ˆ |
лежит в пределах |
||||
ность эмпирической оценки M (x) |
||||||
|
|
= |
|
ˆ |
|
≤ 3σn , |
|
|
|
|
|||
|
|
|
M (x) − M (x) |
|
||
|
|
|
|
|
|
|
где σn = 
D(x) / n – среднее квадратичное отклонение оценки
D(x) – дисперсия усредняемой выборочной с. в. x; n – объём выборки.
(1.41)
ˆ
M (x);
25