ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 06.04.2021
Просмотров: 2196
Скачиваний: 4
УДК 517.958
МЕТОД МАКСИМАЛЬНОГО СЕЧЕНИЯ
С ИСПОЛЬЗОВАНИЕМ ВЕТВЯЩИХСЯ
МАРКОВСКИХ ЦЕПЕЙ ДЛЯ РЕШЕНИЯ
УРАВНЕНИЯ ПЕРЕНОСА
А.С. Жуплев
Институт прикладной математики ДВО РАН
Россия, 690041, Владивосток, Радио 7
Дальневосточный федеральный университет
Россия, 690091, Владивосток, Суханова, 8
E-mail:
zhuplev@gmail.com
Ключевые слова:
Алгоритмы Монте-Карло, уравнение переноса излуче-
ния, метод максимального сечения
В работе рассматриваются модификация метода максимального сечения
для решения стационарного уравнения переноса излучения в трехмерной
неоднородной среде. Модификация основана на применении метода Монте-
Карло к суммированию ряда Неймана для решения уравнения переноса.
Идея модификации основана на использовании ветвящихся цепей Маркова.
Проводится численное сравнение этих алгоритмов.
Введение
Не смотря на многолетнюю историю, использование метода максимального се-
чения продолжает привлекать внимание специалистов, как с точки зрения его обос-
нования, так и с точки зрения его приложений [1-3]. Причем, это касается не только
теории переноса излучения, но и теории управления, массового обслуживания и др..
Одним из универсальных способов уменьшения дисперсии в задачах суммирования
рядов Неймана, в том числе и для решения уравнения переноса, является использо-
вание ветвящихся марковских цепей [4]. В настоящее время такой подход активно
применяется в различных задачах вычислительной математики и математической
физики [5-7]. В данной работе проводится численное сравнение метода максималь-
ного сечения с использованием ветвящихся траекторий с ветвлением в узлах фик-
тивного и истинного рассеяния частиц с традиционными методами максимального
сечения для стационарного уравнения переноса излучения.
1.
Постановка задачи
Рассматривается стационарное монохроматическое уравнение переноса излуче-
ния следующего вида [8-10]
ω
· ∇
r
f
(
r, ω
) +
µ
(
r
)
f
(
r, ω
) =
µ
s
(
r
)
Z
Ω
P
(
r, ω
·
ω
0
)
f
(
r, ω
0
)
dω
0
+
J
(
r, ω
)
,
(1)
Сборник материалов XXXVII Дальневосточной Математической Школы-Семинара
имени академика Е.В. Золотова, Владивосток, 8 – 14 сентября 2013 г.
86
где
r
= (
r
1
, r
2
, r
3
)
∈
G
,
G
– выпуклая ограниченная область в пространстве
R
3
,
ω
∈
Ω =
{
ω
∈
R
3
:
|
ω
|
= 1
}
,
Ω
– единичная сфера в
R
3
. В уравнении (1) функ-
ция
f
(
r, ω
)
интерпретируется как плотность потока частиц в точке
r
в направлении
ω
, функция
J
(
r, ω
)
описывает распределение внутренних источников излучения в
среде,
µ
(
r
)
,
µ
s
(
r
)
– называются коэффициентами полного ослабления и рассеяния,
а функция
P
(
r, ω
·
ω
0
)
– индикатрисой рассеяния. Для характеристики неоднород-
ности среды
G
, в которой изучается процесс переноса излучения, обычно вводится
в рассмотрение некоторое разбиение
G
0
области
G
[10]. Множество
G
0
является
объединением конечного числа областей
G
0
=
p
[
i
=1
G
i
,
G
i
∩
G
j
=
∅
,
i
6
=
j,
таких, что
G
0
=
G
. Области
G
i
можно интерпретировать как некоторые части неод-
нородной среды
G
, заполненные
i
-м веществом. Относительно геометрии множе-
ства
G
0
традиционно предполагается выполнения свойства обобщенной выпуклости
[13,14]. Согласно этому свойству любой луч
L
r,ω
=
{
r
+
tω, t
>
0
}
, исходящий из
точки
r
∈
G
0
в направлением
ω
∈
Ω
пересекает
∂G
0
в конечном числе точек. Все
функции
µ
,
µ
s
,
J, P
, – неотрицательные,
µ
s
(
r
)
6
µ
(
r
)
, функция
P
(
r, ω
·
ω
0
)
удовле-
творяет условию нормировки:
Z
Ω
P
(
r, ω
·
ω
0
)
dω
0
= 1
.
Обозначим
Γ
±
=
{
(
z, ω
)
∈
∂G
×
Ω :
L
z,
∓
ω
∩
G
0
6
=
∅}
и присоединим к уравнению (1) граничное условие
f
(
ξ, ω
) =
h
(
ξ, ω
)
,
(
ξ, ω
)
∈
Γ
−
.
(2)
Неотрицательная функция
h
(
ξ, ω
)
∈
C
b
(Γ
−
)
описывает входящий в среду
G
поток
излучения. Введем обозначения
d
(
r, ω
)
есть расстояние от точки
r
∈
G
до границы
∂G
=
G
\
G
в направлении
ω
(
Lf
)(
r, ω
) =
ω
· ∇
r
f
(
r, ω
) +
µf
(
r, ω
)
,
µ
= sup
r
∈
G
0
µ
(
r
)
,
(
Sf
)(
r, ω
) =
µ
s
(
r
)
Z
Ω
P
(
r, ω
·
ω
0
)
f
(
r, ω
0
)
dω
0
+ (
µ
−
µ
(
r
))
f
(
r, ω
)
.
Тогда уравнение (1) можно переписать в следующем виде
Lf
=
Sf
+
J,
(3)
Под решением прямой задачи (1),(2) будем понимать функцию
f
(
r, ω
)
, которая удо-
влетворяет соотношениям
(
Lf
)(
r, ω
) = (
Sf
)(
r, ω
) +
J
(
r, ω
)
,
(
r, ω
)
∈
G
0
×
Ω
,
f
(
ξ, ω
) =
h
(
ξ, ω
)
,
(
ξ, ω
)
∈
Γ
−
.
Сборник материалов XXXVII Дальневосточной Математической Школы-Семинара
имени академика Е.В. Золотова, Владивосток, 8 – 14 сентября 2013 г.
87
2.
Ряд Неймана для нахождения решения
краевой задачи
Введём функции, называемые оптическими расстояниями [9, 10]
τ
(
r, ω
) =
d
(
r,
−
ω
)
Z
0
µ
(
r
−
tω
)
dt,
τ
(
r, ω, t
) =
t
Z
0
µ
(
r
−
t
0
ω
)
dt
0
,
Поскольку
µ
s
(
r
)
6
µ
(
r
)
, то решение уравнения (2) может быть найдено в виде ряда
Неймана
f
=
∞
X
m
=0
(
AS
)
m
f
0
,
(4)
где
f
0
(
r, ω
) =
h
(
r
−
d
(
r,
−
ω
)
ω, ω
) exp(
−
µd
(
r,
−
ω
)) + (
AJ
)(
r, ω
)
(
Aφ
)(
r, ω
) =
d
(
r,
−
ω
)
Z
0
exp(
−
µt
)
φ
(
r
−
tω, ω
)
dt.
Ряд (4) равномерно сходится со скоростью
q
m
, где
q
=
k
AS
k
6
sup
(
r,ω
)
∈
G
0
×
Ω
µ
−
µ
(
r
) +
µ
s
(
r
)
µ
1
−
e
−
µd
(
r,
−
ω
)
3.
Метод Монте-Карло
Опишем модификацию вычисления суммы ряда (4), которая соответствуют ме-
тоду максимально сечения с использованием ветвящихся марковских цепей. Пусть
M
– число членов усеченного ряда Неймана (4), а
N
– число моделируемых траек-
торий, тогда для решения
f
можно использовать следующую оценку
f
(
r, ω
)
≈
f
M
(
r, ω
)
,
(5)
f
msb
(
r, ω
) =
1
µ
(1
−
exp(
−
µd
(
r,
−
ω
)))
×
×
1
N
N
X
n
=1
[
µ
s
(
r
n
)
f
m
−
1
(
r
n
,
e
ω
n
) + (
µ
−
µ
(
r
n
))
f
m
−
1
(
r
n
, ω
n
)]
,
(6)
где
m
= 1
, . . . , M
и
f
0
(
r, ω
) =
f
0
(
r, ω
)
. В (19) точки
r
n
определяются следующей
формулой
r
n
=
r
−
ωt
n
,
(7)
где
t
n
— реализация случайной величины, распределенной на отрезке
[0
, d
(
r,
−
ω
)]
с
плотностью
µ
exp(
−
µt
)
1
−
exp(
−
µd
(
r,
−
ω
))
.
(8)
Для определения вектора
ω
используем следующие расчетные формулы. Пусть
α
—
равномерно распределенная случайная величина на
[0
,
1]
. Тогда, при
α
>
µ
−
µ
(
r
)
µ
−
µ
(
r
) +
µ
s
(
r
)
(9)
e
ω
=
Qv
(10)
Сборник материалов XXXVII Дальневосточной Математической Школы-Семинара
имени академика Е.В. Золотова, Владивосток, 8 – 14 сентября 2013 г.
88
где
v
=
cos(
ϕ
)
√
1
−
ν
2
sin(
ϕ
)
√
1
−
ν
2
ν
(11)
Q
=
ω
x
ω
z
√
1
−
ω
z
−
ω
y
√
1
−
ω
z
ω
x
ω
y
ω
z
√
1
−
ω
z
ω
x
√
1
−
ω
z
ω
y
−
√
1
−
ω
z
0
ω
z
,
ω
z
6
= 1
1
0
0
0
1
0
0
0
1
,
ω
z
= 1
(12)
а при
α
n,i
<
µ
−
µ
(
r
)
µ
−
µ
(
r
) +
µ
s
(
r
)
(13)
ω
n,i
=
ω.
(14)
В матрице (12) угол
ϕ
есть реализация равномерно распределенной на промежутке
[0
,
2
π
)
случайной величины, а
ν
реализация случайной величины распределенной
плотностью вероятности
P
(
r, ν
)
,
ν
∈
[
−
1
,
1]
.
4.
Численные эксперименты
Оценка работоспособности и точности работы алгоритмов производилась на
ряде численных экспериментов. Для эксперимента были реализованы три метода
Монте-Карло. Первый
(
ms
)
–– стандартный метод максимального сечения, также
описан в [12], Второй способ
(
msb
)
–– метод максимального сечения, с использовани-
ем ветвящихся марковских цепей., третий
(
s
)
— стандартный способ, описан в [12],
второй и третий — методы Монте-Карло с использованием алгоритма максимально-
го сечения. Во всех экспериментах область
G
представляет собой шар единичного ра-
диуса с центром в начале координат. Решение уравнение переноса
f
(
r, ω
)
находится
в точке
r
= (0
.
5
,
0
.
5
,
0
.
5)
,
ω
= (1
.
0
,
1
.
0
,
0
.
0)
. Разбиение
G
0
области
G
состоит из объ-
единения конечного числа шаров (включений).
G
i
=
r
:
|
r
−
r
i
|
< R
i
,
i
= 2
, . . . , p
и
множества постоянными характеристиками
µ
1
,
µ
s
1
, и множества
G
1
=
G
\
p
S
i
=2
G
i
некоторым количеством включений (также в виде шара) с постоянными характери-
стиками
µ
k
,
µ
sk
. Коэффициенты
µ
,
µ
s
внутри
G
i
постоянны и равны соответственно
µ
i
и
µ
s,i
,
i
= 1
, . . . , p
. Количество членов ряда Неймана
M
в каждом эксперимен-
те варьируется с целью анализа качества работы. Тестирование производилось на
аналитических примерах и вычисления осуществлялись на персональном компью-
тере с процессором Intel Core i7 (3.40 GHz) и 8Гб оперативной памяти. В каждом
эксперименте находились относительная погрешность, дисперсия оценки и время (в
микросекундах), затраченное на выполнение алгоритмов. Опуская зависимость от
переменых
(
r, ω, N, M
)
через
f
s
,
f
ms
и
f
msb
будем обозначать статистические оценки
для суммы
f
рядов Неймана для соответствущего метода решения уравнения пере-
носа. Через
ε
m
s
,
Df
ms
,
s
(
f
ms
)
,
ε
msb
,
Df
msb
,
s
(
f
msb
)
,
ε
s
,
Df
s
,
s
(
f
s
)
обозначим относи-
тельную погрешность приближенного решения, дисперсию оценки и трудоемкость
соответствующего метода Монте-Карло. Согласно [11], произведение
Dξ
случайной
величины
ξ
на среднее время
t
расчета одного выборочного значения
ξ
называется
Сборник материалов XXXVII Дальневосточной Математической Школы-Семинара
имени академика Е.В. Золотова, Владивосток, 8 – 14 сентября 2013 г.
89
трудоемкостью
и является критерием оценки качесва алгоритма. Сравнение алго-
ритмов проводилось на аналитическом примере, который соответствует релеевско-
му рассеянию в среде. Результаты численного эксперимента для
10
-ти включений
в области приведены в таблице 1. В этом эксперименте коэффициенты ослабления
и рассеяния имели следующие значения:
µ
(
r
) = 4
,
µ
s
(
r
) = 2
в основной среде, и
µ
(
r
) = 2
,
µ
s
(
r
) = 1
– во включениях. Количество траекторий предполагалось рав-
ным
N
= 10
6
.
Таблица 1.
таблица 1
M
ε
ms
Df
m
s
s
(
f
ms
)
ε
msb
Df
msb
s
(
f
msb
)
ε
s
Df
s
s
(
f
s
)
3
10
.
6620
0
.
2265
4
.
0990
10
.
6658
0
.
2266
4
.
1709
10
.
6641
0
.
2266
24
.
8816
4
4
.
9321
0
.
1016
2
.
4253
4
.
9436
0
.
1019
2
.
5116
4
.
9369
0
.
1017
14
.
9792
5
2
.
1570
0
.
0753
2
.
2660
2
.
1536
0
.
0755
2
.
3066
2
.
1641
0
.
0754
13
.
8580
6
0
.
8698
0
.
0725
2
.
6121
0
.
8775
0
.
0725
2
.
6476
0
.
8651
0
.
0723
15
.
8910
7
0
.
3311
0
.
0740
3
.
0932
0
.
3199
0
.
0738
3
.
1556
0
.
3154
0
.
0737
18
.
8537
8
0
.
1101
0
.
0752
3
.
5767
0
.
1053
0
.
0752
3
.
6965
0
.
1226
0
.
0753
22
.
0282
9
0
.
0437
0
.
0759
4
.
0746
0
.
0307
0
.
0758
4
.
1839
0
.
0332
0
.
0761
25
.
0274
10
0
.
0102
0
.
0761
4
.
5371
0
.
0016
0
.
0761
4
.
7024
0
.
0066
0
.
0760
27
.
7672
5.
Заключение
В работе проведен сравнительный анализ алгоритмов. Анализ таблицы 1 пока-
зывает, что при одинаковом количестве членов ряда Неймана методы максимально-
го сечения проигрывают в точности. Данная проблема широко известна, поскольку
скорость сходимости ряда Неймана для методов максимального сечения меньше.
Данная проблема исчезает при уменьшении вариации коэффициента
µ
(
r
)
в области
G
. При одном и том же
M
методы максимального сечения, в целом выигрывают по
трудоёмкости (являются более экономичными). Это особенно проявляется в ситу-
ации, когда структура неоднородной среды усложняется, и количество включения
увеличивается. В частности, уже при
10
-ти включениях методы максимального се-
чения выигрывают более чем в
2
раза.
Список литературы
1.
Михайлов Г.А., Аверина Т.А. // Алгоритм "максимального сечения"в методе Монте-
Карло / Доклады Академии наук. 2009. Т. 428, Вып. 2. С. 163-165.
2.
Аверина Т.А. // Методы статистического моделирования неоднородного пуассонов-
ского ансамбля / Сиб. журн. вычисл. матем. 2009. Т. 12, Вып. 4. С. 361-374.
3.
Аверина Т.А., Михайлов Г.А. // Алгоритмы точного и приближенного статистиче-
ского моделирования пуассоновских ансамблей / Журнал вычисл. матем. и матем.
физики. 2010. Т. 50, Вып. 6. С. 1005-1016.
4.
Ермаков С.М. // Метод Монте-Карло в вычислительной математике / Санкт-
Петербург. Бином. Лаборатория знаний. 2009.
5.
Бреднихин С.А., Медведев И.Н., Михайлов Г.А. // Оценка параметров критичности
ветвящихся процессов методом Монте-Карло / Журнал вычислительной математики
и математической физики. 2010. Т. 50, Вып. 2. С. 362-374.
6.
Медведев И.Н., Михайлов Г.А. // Исследование весовых алгоритмов метода Монте-
Карло с ветвлением / Журнал вычислительной математики и математической фи-
зики. 2009. Т. 49, Вып. 3. С. 441-452
Сборник материалов XXXVII Дальневосточной Математической Школы-Семинара
имени академика Е.В. Золотова, Владивосток, 8 – 14 сентября 2013 г.