ДОКЛАДЫ АКАДЕМИИ НАУК, 2015, том 464, № 4, с. 401-405
= МАТЕМАТИКА
УДК 519.676
АЛГОРИТМ МЕТОДА МОНТЕ-КАРЛО ДЛЯ ОЦЕНКИ УГЛОВОГО РАСПРЕДЕЛЕНИЯ РАССЕЯННОГО ПОЛЯРИЗОВАННОГО ИЗЛУЧЕНИЯ НА ОСНОВЕ ОРТОГОНАЛЬНОГО РАЗЛОЖЕНИЯ
© 2015 г. Член-корреспондент РАН Г. А. Михайлов, Н. В. Трачева, С. А. Ухинов
Поступило 18.05.2015 г.
Разработаны алгоритмы метода Монте-Карло для исследования угловой зависимости интенсивности прошедшего через слой вещества излучения на основе разложения плотности соответствующего углового распределения квантов по полиномам, ортонормированным с ламбертовским весом; оптимизация осуществляется по критерию трудоемкости вычислений. Алгоритм позволяет достаточно точно оценить даже малое влияние поляризации, а также отклонение изучаемого углового распределения от ламбертовского.
DOI: 10.7868/S0869565215280063
В настоящей работе рассматривается математическая модель односкоростного процесса переноса частиц (квантов излучения), в основе которой лежит цепь Маркова "столкновений", разделенных "свободными пробегами", которые имеют неоднородное экспоненциальное распределение с коэффициентом <(г), г е В3. Предполагается, что ст(г) = <5(г) + сте(г), где ст^(-) — коэффициент рассеяния, а сте(-) — коэффициент поглощения, причем вне среды < = <с > 0 для нормировки распределения свободного пробега. Распределение направления ю после рассеяния определяется индикатрисой ^(ю', ю) (см., например, [1]). Для учета поляризации с частицей ассоциируется вектор Стокса (I, 0, и, V), первая компонента которого определяет интенсивность излучения; при этом индикатриса заменяется на матрицу рассеяния. В скалярной модели 0 = и = V = 0 [1].
Для определенности будем рассматривать перенос частиц через плоский слой 0 < г < Н от расположенного на границе г = 0, направленного вдоль Ог источника. Практически важным является численное исследование интенсивности Ф(ц, Н) проходящего излучения, где ц = юг. Известно [1], что
Ф(ц, Н) = ц' Н , причем в скалярном случае ц
Институт вычислительной математики и математической геофизики Сибирского отделения Российской Академии наук, Новосибирск
Новосибирский государственный университет E-mail: gam@sscc.ru; sau@sscc.ru
F(^, H) — плотность распределения по ц вылетающих из слоя частиц. Введем обозначения: PH =
JF(x
, H)dx, fx) = fx, H) = PH F(x, H), ф(х) = f(x)
x
Работа посвящена решению задачи минимизации трудоемкости статистического моделирования для оценки углового распределения интенсивности проходящего излучения, т.е. функции f( х H)
--—-, 0 < х < 1, с целью анализа изменения этой х
функции при увеличении H, а также при введении поляризации в математическую радиационную модель, как указано выше.
Рассмотренный далее алгоритм основан на методе Н.Н. Ченцова для оценки неизвестной плотности вероятностей [2], причем используется модификация этого метода, предложенная в [3]. Пусть у,(х) — система ортонормированных с ве-
1
сомр(х) на (0, 1) полиномов, т.е. |p (х)у((х)уу(х)йгх =
0
= djj, и для плотности вероятностейДх), согласно [3], выполняется соотношение
f( x) = p (x) ^ a^i (x),
i = 0
где
at = Jf(x )v;(x) dx = EVi(^),
0
да
2
401
ц — случайная величина, распределенная с плотностью /(х).
Равенства (1) определяют алгоритм метода Монте-Карло для оценки величин а, и функции/х). Они фактически являются исходным пунктом метода Н.Н. Ченцова, а точнее, его модификации, так как в
да
[2] рассматривается представление /х) = ^ а1 у,(х)
I = о
и, соответственно, а, = Е(р(ц)у;(ц)). Как видно из дальнейшего, модификация, основанная на (1) является эффективной для оценки и сравнительного анализа функции Ф(х, Н) при Н > 1.
Известно (см., например, [4]), что для больших значений Н при слабом поглощении частиц в среде плотность /(х) близка к плотности Ламберта
/0(х) = 2х. Поэтому, с учетом равенства ф(х) = ,
х
здесь целесообразно положить р(х) = х. При этом полиномы у,-(х) можно получить нормировкой полиномов у1 (х), полученных заменой переменных у = 2х — 1 из полиномов уI (у), ортогональных с весом у + 1 на отрезке (—1, 1), которые представляют собой частный случай полиномов Якоби (см., например, [5]).
Для полиномов у,-(х) получен следующий явный вид:
щенная оценка коэффициента a, имеет вид
(x) = Jli + 2 X
k = 0
(-1) k ( 2 i + 1 - k) ! , - k ( i - k)!k!(i + 1 - k)! '
Введем обозначения: N — объем выборки случайных траекторий частиц, ^ — случайное число частиц, достигших границы г = Н, 8Н — средняя трудоемкость моделирования траектории. Случайная оценка функции ф(х, Н), согласно (1), строится следующим образом:
Ф(х, Н) - ^ а у,(х) =
I = о
п
= Фп( х )-ф п (х) = ^ а, у, (х),
I = о
Nн
где а,- = яН ^ у, (р.у). При этом Ефп (х) = ф„(х) +
1 = о 1
+ О^-1) и |хфп (х)йх = 1. При весовом моделиро-
о
вании процесса переноса излучения, в том числе при учете поляризации, асимптотически несме-
N„
X Qj )
а, = ■
N„
где Qj — вес частицы (первой
j = о
компоненты вектора Стокса) при вылете из рассматриваемого слоя в направлении
Согласно [2] с точностью до величины О^-1) имеем следующее выражение для среднеквадратичной с весом р(х) погрешности оценки фп:
ElФ - ф J2 = ¡Р(x)Офn(x)dx + ||ф - ф
' = A n + A„
С другой стороны, трудоемкость реализации фn определяется величиной S = N(SH + snPH), где s — это фактически трудоемкость операции умножения для вычисления очередной степени ц' + 1 = ц • ц', так как определение значений а, можно осуществить в конце моделирования; следовательно s « 1. Трудоемкость оценок величин Da, здесь не учитывается, так как их точность не играет роли в последующем анализе.
Далее рассматривается следующая задача минимизации относительной трудоемкости S численно-статистической оценки функции ф(х, H):
S = N(SH + snPH)^ min, An + An = 52,
n, N
или (n0, N0) = arg minS, при условии An + An = 52.
n, N
Теорема. Пусть выполняются следующие условия:
1) производная ф(г)(х) порядка r непрерывна на [0, 1];
2) 0 < f1 < ¿Щ < f2 < +да.
Р (x)
Тогда
1) равенство An + An = 52 достигается, если
n = n1 = (2C1)1/2r5-1/r, N = N1 = 2C0n1 pH 5-2,
где fi < C0 < /2 (причем fi, /2 — f1, f2 при 5 ^ 0), а C1 определяется из соотношения An < C1n-2r.
При этом S = S1 = (SH + sn1PH) • 2C0n1 PH 5-2;
2) асимптотически при H ^ да оптимальные значения величин n, N, S равны
П0 = (2r + 1)1/2rCj/2r 5-1/r, N0 = 2r+1 CnpH 5-2,
2r
S
S0 = SHN0, причем 1 < — =
2r + 1
2 r + 1 2 r
2r < 2.
0
n
0
Доказательство. Соотношение Дп < С1п—2г, как и оценка \а,-\ < Сгг, ( < п, следует из общей теории аппроксимации гладких функций суммами ортогональных многочленов (см. [6, глава IV, § 7]). Далее по аналогии с [2] получается выражение
~ Р1 " (11
Дп = -Н- ^ о + о(-) , причем здесь
Таблица 1. Влияние поляризации на поток P
н
D у (и) = jfx )V2 (x) dx - a2
1 f
0
С учетом соотношений jp (x) у2 (x)dx = 1 и a2 < Ci
= fv^P(x)V2(x)dx - a2' Jp( x)
отсюда получаем A„
ConPH
N
соответственно
условию 2) теоремы.
Таким образом, задача оптимизации относительной трудоемкости приобретает вид S = N(SH +
+ snPH) ^ min, при условии
n, N
Co nPH + C = S2.
N
2 r
(2)
Приравнивая слагаемые суммы в (2) к — , получаем первое утверждение теоремы.
Далее при H ^ да имеем PH ^ 0 и величина SH + snPH асимптотически заменяется на SH, после чего задача оптимизации S формулируется в виде N ^ min при условии (2), из которого получаем
P hN =
Con
C
52 - Cn2
-1 с 2 ^ - 2r - 1
n о -C1n
Минимизируя это выражение по п, после дополнительных вычислений получаем второе утверждение теоремы, завершая таким образом доказательство.
Заметим, что иногда практически использу-
ется оценка ф(х)
Фh (х) = , где fh (х) - ги-
x
стограмма с шагом к для выборки {цк}, к = 1, 2, ..., Ын. Однако соответствующие оценки практически важных функционалов |ф (x)g(x)dx, g(x) > б > 0
имеют бесконечную дисперсию. Возможно также использование "локальных оценок" [1], т.е. оценок вклада в значения {ф^,)} с шагом к непосредственно от каждого столкновения. Если гладкое
H AP PH aph ) PH
2 -8.3 • 10-5 5.4 ■ 10-7 0.69596
5 1.3 • 10-4 1.7 ■ 10-6 0.60416
10 2.5 ■ 10-4 4.8 ■ 10-6 0.41314
20 2.9 ■ 10-4 1.6 ■ 10-5 0.25008
2r
восполнение получаемой таблицы {ф ^^ осу-
п
ществлять на основе разложения ф^) = ^ а1 у ((г),
г = 1
п = к-1, то получаем соотношение, близкое к (2). Однако величина Б = Бн + «Рнк-1 при этом заменяется на Б = Бн + Б'нк-1, где Б'н > Бн, т.е. "локальные оценки" не эффективны для н > 1.
При проведении численных расчетов использовалась матрица аэрозольного рассеяния, рассчитанная по теории Ми для следующих параметров аэрозольной среды [7]: коэффициент преломления частиц п = 1.331 — I -1.3 • 10-4 (вода), распределение частиц по размерам логнормально
с плотностью /(г) = 1 ехр (—— 1п2 ( -)) , г е (0, 10),
г К 2ст] Ч"
г& = 0.12, ст = 0.5, длина волны излучения равна 0.65 мкм. Средний косинус угла рассеяния = = 0.7292.
Проводилось исследование влияния поляризации на интегральный поток Рн, проходящего через слой излучения, и на коэффициенты разложения а1 из (1). Проведенные прецизионные расчеты (моделировалось 1010 траекторий) показали, что для слоев оптической толщины от н = 2 до н = 20 влияние поляризации на Рн возрастает с ростом толщины слоя; кроме того, подтверждается известное для случая ст = ст8 асимптотическое соотношение Рн = 0(н—1) [4]. В табл. 1 приведены статистически значимые оценки относительных (по отношению к Рн) разностей ДР потоков с учетом поляризации и без, а также их среднеквадратичные отклонения стя{ ДРн).
В табл. 2 приведены значения относительной разности Да,- оценок коэффициентов с учетом поляризации и без, а также оценки соответствующих среднеквадратичных отклонений ст^Да,-). Для слоя н = 20 влияние поляризации на коэффициенты а-оказалось статистически незначимым.
В табл. 3 приведены статистически значимые оценки коэффициентов а1 для н = 2, 5, 10. В связи с быстрым убыванием этих оценок и возрастанием относительной погрешности оценок Да,- по (
0
2
Таблица 2. Влияние поляризации на коэффициенты а, (а0 = а0 = ^2 )
Н = 2 Н = 5 Н = 10
- Аа, Аа, МАа,-) Аа,-
1 8.5 ■ 10-3 4.1 ■
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.