научная статья по теме АЛГОРИТМ МЕТОДА МОНТЕ-КАРЛО ДЛЯ ОЦЕНКИ УГЛОВОГО РАСПРЕДЕЛЕНИЯ РАССЕЯННОГО ПОЛЯРИЗОВАННОГО ИЗЛУЧЕНИЯ НА ОСНОВЕ ОРТОГОНАЛЬНОГО РАЗЛОЖЕНИЯ Математика

Текст научной статьи на тему «АЛГОРИТМ МЕТОДА МОНТЕ-КАРЛО ДЛЯ ОЦЕНКИ УГЛОВОГО РАСПРЕДЕЛЕНИЯ РАССЕЯННОГО ПОЛЯРИЗОВАННОГО ИЗЛУЧЕНИЯ НА ОСНОВЕ ОРТОГОНАЛЬНОГО РАЗЛОЖЕНИЯ»

ДОКЛАДЫ АКАДЕМИИ НАУК, 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 = о

где а,- = яН ^ у, (р.у). При этом Ефп (х) = ф„(х) +

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 рублей.

Показать целиком

Пoхожие научные работыпо теме «Математика»