Расчет параметров временной асимптотики выходящего из полубесконечного слоя поляризованного излучения Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 13, Специальный выпуск 4, 2008

Расчет параметров временной асимптотики выходящего из полубесконечного слоя поляризованного излучения*

Н.В. Трачева, С. А. Ухинов, A.C. Чимаева Институт вычислительной математики и математической геофизики, Новосибирск, Россия e-mail: tnv@osmf.sscc .ru, sau@sscc.ru, chimaevaa@gmail.com

Based on the method of parametric differentiation of a special representation for a solution of the nonstationary vector transfer equation, the constants of exponential and power temporal asymptotes are estimated using the Monte Carlo method.

Введение

Существует ряд задач теории переноса, при решении которых исследователей интересует асимптотическое поведение потоков излучения на больших временах в еветораеее-ивающих и поглощающих средах. Известно, что для неполяризованного излучения эта асимптотика при выполнении довольно общих условий экспоненциальная. Параметром экспоненциальной асимптотики при этом является ведущее характеристическое число А* соответствующего однородного стационарного уравнения переноса [1].

В работе [2] решаются вопросы, связанные с распространением этого утверждения на случай поляризованного излучения. Для однородной бесконечной среды, т. е. для полного пространства, это достигается довольно просто на основе теории весового моделирования процесса переноса, причем оказывается, что А* = -acv независимо от типа поляризации. Здесь ас — сечение поглощения, v — скорость частиц ("фотонов").

Специалистам хорошо известно также, что А* = — асv и для полупространства. Эвристически это достаточно ясно из соображений непрерывности величины А* как функции оптической толщины плоского слоя.

Для ограниченной среды теоретически нерешенным оказался вопрос о связи значений А^ и AV, т. е. значений параметра асимптотики соответственно для скалярного и векторного вариантов с одинаковой индикатрисой рассеяния Кц{и,и'). Это фактически вопрос о соотношении скоростей деполяризации потока частиц и его перехода к экспоненциальной асимптотике.

В работе [2] с помощью расчетов методом Монте-Карло показано, что для ограниченных сред в случае молекулярного и аэрозольного рассеяния А^ = А^, т.е. деполяризация потока частиц несколько запаздывает относительно перехода к асимптотике.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 06.01-00046-а, № 07.01-00673-а, № 09.01-00035-а) и INTAS (Ref. N 05-1000008-8024).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

В настоящей работе рассматривается задача определения параметров временной асимптотики потока 3(г) поляризованного излучения, выходящего из полубесконечного слоя рассеивающего и поглощающего вещества при освещении его внешним плоскопараллельным источником. По аналогии со скалярным вариантом задачи [3] в предположении, что

3(г) х Гаех*(0.1)

строится алгоритм метода Монте-Карло для нахождения параметров степенной а и экспоненциальной А* асимптотик в векторном случае,

1. Система уравнений переноса с учетом поляризации. Временная постоянная

Для описания поляризационных свойств света используется вектор Стокса I = (I ,(,и,У^ , компоненты которого определяют интенсивность, степень и плоскость поляризации, а также степень эллиптичности излучения. Математическая модель переноса поляризованного излучения строится на основе феноменологического предположения о том, что в результате рассеяния ассоциируемый с "фотоном" вектор Стокса преобразуется заданной матрицей рассеяния.

Рассматривается интегродифференциальное уравнение переноса излучения с поляризацией:

иУФ + аФ = ! а8(т)Р(и' ,и,т)Ф(т,и') ди' + ?о(т,и), (1.1)

п

где Ф = (Ф1, Ф2, Фз, Ф4)Т _ вектор-функция плотности потока частиц ("фотонов"); П — пространство единичных векторов направления; и € П, т € О С К3; Р(и',и,т) — матричная функция рассеяния; а = а(т) — полное сечение, а = а8 + ас, ас — сечение поглощения, а8 — сечение рассеяния; ^ = (/о^,/о2),/о3),/о4))т _ вектор-функция плотности распределения источника частиц,

А

иУФх(т,и) + а Ф\(т,и) = ! а8(т)Р (и' ,и,т)Ф\(т,и') ¿и', (1,2)

п

где а = а8 + а0 ас = ас + А/у.

Вектор-функция плотности столкновений р(т,и) с вектор-функцией плотности потока связана соотношением

^ = (аФ1, аФ2, аФз, аФ4)т = (ф1, ф2, фз, ф4)Т.

Известно, что она удовлетворяет интегральному уравнению переноса с учетом поляризации:

ф(х) = J К(х',х)р(х') дх' + ¥(х) (1,3)

X

или

фг(х) = / кц (х',х)ф (х') дх' + /г(х), % =\,. . ., 4,

¿=1 X

где x = (r, и) G R3 х П; f = (f1, f2, f3, f4)T — вектор-функция плотности распределения начального столкновения.

Здесь

, q(r/)exp(-r(r',r))a(r)P(u',u,r') ( r - r' \ K(x,x) =-——-, (1,4)

где

|r'-r|

a(r) f q(r) = s, , r(r',r;cj)= / a(r' + tu)dt. ^(r) J

0

В нестационарном случае вводится дополнительная фазовая временная координата t и ядро (1.4) домножаетея на 6 (t — (r — r')/v).

После подстановки диагональной матрицы PR = diag(Rn, 0, 0, 0) вместо P в систему уравнений переноса с поляризацией получаем уравнение переноса излучения без поляризации, В этом случае наибольшее вещественное характеристическое значение Л* задачи (1.2) называется временной постоянной, так как она определяет асимптотику потока частиц I(r, и) ~ C(r, и) ехр(ЛЧ) [1],

2. Решение системы интегральных уравнений методом Монте-Карло

Методом Монте-Карло обычно оценивают линейные функционалы от решения рассматриваемого интегрального уравнения. Ниже приводится общий алгоритм метода Монте-Карло для оценки таких функционалов в случае системы интегральных уравнений второго рода [4-7],

Для вычисления искомых функционалов используется "оценка по столкновениям", В случае процесса переноса излучения с учетом поляризации будем оценивать линейные функционалы вида

j=ол h ) = £ у

<Pi(x)hi(x).

i=i x

Здесь Ь — вектор-функция с абсолютно ограниченными компонентами, т. е, ЛДх) € Ь^ при г = 1,..., 4.

Рассмотрим цепь Маркова с начальной плотностью п(х) и субстохастической переходной плотностью р(х',х) такими, что

п(х) = 0, если ¥(х) = 0, и р(х', х) = 0, если К(ж',ж) = 0,

для которой N — случайный номер последнего состояния.

Введем также вспомогательный случайный вектор "весов" Qn:

л{г) _ МХо) п(г) _ кч(Хп-1}Хп)

п(г) _ 0) n{i) _ ST^ n(j) ^

Vo - / ч> Чп - / у Уга-Г

п(хо) ' П ^ П р(хп-1,хп)

Аналогично тому, как это делается для одного интегрального уравнения, можно показать [5, 6], что

N Г 4

еЕ EQni)hi(xn) = (р, Ь) = 3 = ^, <р*). (2.1)

n=0

i=1

Соотношение (2,1) описывает алгоритм метода Монте-Карло для оценки величины 3, При обосновании этого соотношения используется разложение решения системы уравнений в ряд Неймана [4, 6]. Почленное осреднение первой компоненты суммы из (2,1) допустимо в силу ее неотрицательности, а остальных компонент — вследствие мажорантного свойства первой компоненты [5, 6], В (2,1) — это решение сопряженной к (1,3) системы.

Таким образом, получаем "оценку по столкновениям" [4-6]

N £

п=0

N

ЬЫ, Ее = 3, (2.2)

п=0

.1=1

а также оценки параметрических производных функционала 3:

Е^М = з М.

N

Ат) = = \ -

4 Р)\т ¿_

дЛг

п=0

\ ^ д €)п , , . .¿=1

Здесь предполагается, что Ь те зависит от Л,

Для исследования вопроса о дисперсии используемых оценок величину £ представим

N ~ ~ ~ Кт (Хп-1,Хп) ~

в виде £ = где £хо = ^ С}пЬ, С}п = —-—-—<50 = — единичная

п=0 р(хп-1,хп)

матрица.

Известно, что Е£Ж0 = ф*(х0) и при выполнении определенных условий ковариационная матрица Ф(х) = Е(£х, ^) удовлетворяет матричному интегральному уравнению

Т/ ч ,/ ч [ КТ(х,у)Ф(у)К(х,у) , , ч

Ф(х) = А(х) + / -1 УУ{ 1 ,уиу, 2.3)

] Р(х,У)

или в операторном виде Ф = А + КРФ, где А = Ьр*т + <^*ЬТ — ЬЬТ. Дисперсия исходной оценки £ связана с ковариационной матрицей Ф следующим соотношением:

= Ех0 (ЯТ ФЫО)) — 32,

а следовательно, из ограниченности матрицы Ф(х) вытекает конечность дисперсии оценки £.

В работах [5, 7] доказано, что если спектральный радиус р(Кр) < 1, то для любой вектор-функции Ь(х) € Ьж, ковариационная матрица Ф(х) удовлетворяет уравнению (2.3), < В [7] также фактически показано, что при выполнении условия р(Кр) < 1 справедливы соотношения Е£(т) = 3(т) и (т) <

В работе [8] на основе теории положительных операторов вычислена величина спектрального радиуса р(Бр) для оператора переноса поляризованного излучения в бесконечной однородной среде. Показано, что величина р(Кр) приближенно равна произведению значения р(Бр) па спектральный радиус оператора, соответствующего переносу излучения без поляризации, и доказано неравенство

2

р(Кр) < р(Бр) 8ПР / --(II, (2.4)

Т,Ш 0 Ч1(Г + и1) ([■ Г,Ш)

где рх(1; г, и) = а (г + и1)ехр(—а (г + и1)) — плотность распределения длины свободного пробега I в точке г в направлении и; рХ? (I; г, и) — плотность распределения длины

пробега в моделируемой цепи Маркова; ^1(г) — вероятность выживания цепи Маркова в точке г,

В случае молекулярного рассеяния [4] в работе [8] вычислено значение р(Бр) = 1.178, а следовательно, дисперсия оценки £ в реальной среде конечна только при достаточно большом поглощении. Например, для "физического" моделирования из (2,4) легко получить оценку

р(Кр) < (1 - р)р(Бр), (2.5)

где р — нижняя граница вероятности поглощения в среде, а значит, р(Кр) < 1 и дисперсия конечна при р > 0.151, Для использованной в [8] модели аэрозольного рассеяния р(Бр) и 1.02077 и р(Кр) < 1 при р > 0.02034.

Если же процесс переноса модифицируется путем замены о — о3, ос — 0 с учетом поглощения весовым множителем в-ас1 [4, 6, 9], то из (2,4) следует, что

р(К„) < ^Р(^) (2.6)

и р(Кр) < 1 при р > 0.082 в случае молекулярного рассеяния и при р > 0.0103 в случае аэрозольного рассеяния, т. е. при относительно меньшем поглощении в среде.

Полученные из неравенств (2.4)-(2.6) выводы о конечности дисперсий оценок функционалов и их производных являются основанием для применимости метода Монте-Карло в том или ином случае.

3. Оценка параметров временной асимптотики методом Монте-Карло

Величину Л* можно вычислять, используя специальный итерационный метод [5|:

4т)( Ао)

Ит/'^ ^=А*-А0. (3.1)

Здесь Л0 — некоторое начальное значение; 3д — произвольный линейный функционал от плотности потока, соответствующего уравнению (1.2); 3(т)(Л0) — производная от 3 по параметру Л в точке Л = Л0,

Данное в [5] обоснование этого метода для скалярного варианта легко распространяется на векторный случай [2]. Однако этот метод не позволяет находить параметр асимптотики а, и мы будем использовать метод, основанный на вычислении параметрических производных по времени [3, 9].

Рассмотрим нестационарный процесс переноса от источника, вообще говоря, поляризованного излучения, причем ¥(г,и,Ь) = f (г, и,£)10, где 10 = (10^0,и0,У0). Справедливо соотношение

г

3 = У У (г, и, ¿)Ь(г, =УУ У f (г0,и0,т )30(г0,и0,£ — т )^г0^и0^т,

я п я п 0

где

30(г0, = У У (г, и, ¿; Г0, ^)Ь(г, и)<1г<1и. яп

Здесь р0(х; то,ио) — векторная плотность столкновений (по аргументу х) от одного столкновения в точке (то,ио, 0), т. е. для функции источника 5(т — то)5(и —

Обозначим через п(то,ио) "оценку по столкновениям" (см, уравнение (2,2)) для функционала

(то,ио) = ^ ^ У (т,и,т; то,ио)Ь(т,и)^и^т. я по

Справедливо следующее утверждение.

Теорема. Пусть точка (то,ио) распределена для, го = 0 с плотностью /х(т,и), причем функция /(п х) (х) абсолютно непрерывна по г во вся,ком конечном интервале времени V (т,и) € Я х П и | (х) | < С/х(т, и) для почти в сех х, т = 0,... ,п. Пусть также УхЕ| п |2 € Ьх(Я х П) и |Зо(х)| < С < Тогда, выполняются соотношения,

J{Ш) (£) = где

= ЕГМ (3-2)

причем 0£(ш) < если р(Кр) < 1.

Теорема доказывается по аналогии со скалярным вариантом, который изучен в [9]. Рассмотрим теперь оценку параметра экспоненциальной временной асимптотики. Как указано выше, в скалярном случае при выполнении довольно общих условий имеет место асимптотическое соотношение

Jо(т, и, г) ~ С(т, и)вхч, г ^ (3.3)

Эти условия, в частности, имеют место для односкоростного процесса переноса в ограниченной среде с функцией источника, достаточно быстро убывающей по времени [1]. Поэтому мы предполагаем, что если /(т, и,г) ехр(—Л*г) —> 0 V (т,и), то выполняется

соотношение (3,3) и

J (г) = Сехч[1 + е(г)], е(г) —► 0, (3.4)

также и в случае поляризованного излучения.

Так же как в [9], можно показать, что функция J' (г) обладает аналогичным свойством, т.е. если /'(т,и,г) ехр(—Л*г) —> 0 V (т,и), то

J'(г) = С\*вхч[1 + £х(г)], £х(г) —► 0. (3.5)

Вследствие соотношений (3.4) и (3.5) величина J'(t)/J(г) для достаточно большого значения г дает оценку временной константы Л*.

В [2] доказано, что главное характеристическое число Л* векторного уравнения (1.2) в пространственно-однородной среде равно — асV. Если сделать очевидное предположение, что в полубесконечном слое это также выполняется, то можно второй параметр а асимптотики (0.1) находить из следующего соотношения:

(1п = *»1.

То есть а ~ (—Л — асv)t, где Л = .]'(г)/-1 (г) = Е££/Е£4 — приближенное значение Л*,

г

4. Результаты численных экспериментов

Рассмотрим полупространство г > 0, заполненное рассеивающим и поглощающим свет веществом, источник излучения находится в точке с координатами (0, 0, 0) и излучает "фотоны" в направлении и = (0, 0,1), причем /0 = (1, 0, 0, 0)т. Имеет место одноекороет-ной (V = 1) процесс переноса частиц, коэффициенты рассеяния и поглощения заданы, полное сечение ослабления о = 1,

В расчетах использовались матрица молекулярного рассеяния и матрица аэрозольного рассеяния, рассчитанная по теории Ми для следующих параметров среды. Коэффициент преломления частиц п = 1.331 — ¿1.3 • 10-4 (вода), распределение частиц по размерам логнормально с функцией распределения

/(г) = - ехр r

V,

2s2 V Гд

r £ (0,10 мкм), rg = 0.12 мкм, s = 0.5,

длина волны излучения равна 0,65 мкм.

Плоскопараллельный приемник регистрирует поток J(t) излучения, выходящего из полубесконечного слоя (освещенность границы z = 0), т. е, определяет среднее значение первой компоненты вектора Стокса для вылетающих "фотонов",

В расчетах методом Монте-Карло использован стандартный алгоритм моделирования траекторий "векторного фотона" с учетом "веса". Использовалась модификация локальной оценки по столкновениям в состоянии после рассеяния с моделированием цепи Маркова без поглощения, т.е. с заменой a — as, ac — 0, Также использовалась модификация "моделирования без вылета" [4, 6],

В качестве временной составляющей плотности распределения первых столкновений взята функция вида f (t) = te-t, t > 0, Соответствующие "оценки" (3,2) для вычисления потока и его производной имеют следующий вид:

N

£(t) = £ A-(wn)At(t„)(t - t„)e-(i-in+1"*>Q«,

n=0

N

£'(t) = £ A-(^n)At(tn)(1 - (t - tn))e-(i-in+in^Q«,

n=0

где Qn определяется выражением

4 P

Q« = exp{—/га<тс}(1 - exp{-/>s}) V -^Q^, Qo = (1, 0, 0, 0)T

Здесь l^ — расстояние от точки рассеяния до границы слоя в направлении wn; ln — длина п-го пробега между столкновениями; tn = (Ln + in)/v, Ln — полная длина пробега до n-го столкновения; A-(w) — индикатор пересечения границы слоя направлением w; At(r) — индикатор интервала (0,t); R11 — индикатриса рассеяния; Piyj — (i,j)-a компонента матрицы рассеяния P(w', w, r), Траектория обрывается при Ln/v > t.

Для проведения расчетов создана параллельная реализация алгоритма. Код программы написан на языке Intel Fortran с использованием MVAPICH на коммуникационной среде InfiniBand [10]. Также использовалась библиотека Intel MKL, Вычисления

выполнялись на многопроцессорной системе НКС-160 (ССКЦ) [11]. Число реализованных траекторий с использованием 100 процессоров составило N = 1010,

Расчеты проводились с помощью двух различных генераторов псевдослучайных чисел: 128-битного конгруэнтного генератора [12] и генератора МТ2203 из библиотеки \ IКI. [13]. Отметим, что результаты расчетов с разными генераторами псевдослучайных чисел статистически не различались.

4.1. Расчет параметров асимптотики

Результаты расчетов параметров асимптотики Л и а для а,3 = 0.9, ас = 0.1 и аэрозольной матрицы рассеяния приведены в табл. 1 и 2 соответственно. Также в таблицах приведены значения среднеквадратичной погрешности результатов.

Заметим, что значения параметров временной асимптотики в случаях с учетом поляризации и без ее учета совпадают с точностью до статистической погрешности.

В табл. 3 и 4 приведены результаты расчетов, полученные для сред с параметрами а3 = 0.95, ас = 0.05 и а3 = 1.0, ас = 0.0 соответственно. Для первой среды погрешности вычисления Л те превосходят 0.22% при моделировании без учета поляризации, 1 % при моделировании с поляризацией. Для второй среды погрешности вычисления 2%

Таблица 1. Параметр Л экспоненциальной асимптотики освещенности, вычисленный методом параметрических временных производных для а3 = 0.9, ас = 0.1

г Без учета поляризации С учетом поляризации

21 -0.14850806 ±0.00002813 -0.14857793 ± 0.00002874

41 -0.13246373 ±0.00004165 -0.13239391 ± 0.00004521

61 -0.12291465 ± 0.00005455 -0.12291154 ± 0.00006531

81 -0.11766099 ±0.00006639 -0.11750441 ± 0.00009052

101 -0.11431486 ± 0.00007772 -0.11447877 ± 0.00012443

121 -0.11185582 ±0.00008862 -0.11205463 ± 0.00016620

141 -0.11031517 ±0.00009888 -0.11020254 ± 0.00021730

161 -0.10927209 ±0.00010877 -0.10903962 ± 0.00029130

181 -0.10822164 ± 0.00011865 -0.10814159 ± 0.00037220

201 -0.10732061 ± 0.00012816 -0.10776328 ± 0.00048264

Таблица 2. Параметр а асимптотики освещенности для а3 = 0.9, ас = 0.1

г Без учета поляризации С учетом поляризации

21 -1.01866919 ±0.00059075 -1 02013657 ± 0.00060360

41 -1.33101302 ±0.00170785 -1 32815022 ± 0.00185356

61 -1.39779342 ±0.00332768 -1 39760411 ± 0.00398409

81 -1.43054017 ±0.00537739 -1 41785750 ± 0.00733217

101 -1.44580067 ±0.00784958 -1 46235579 ± 0.01256675

121 -1.43455391 ± 0.01072290 -1 45861005 ± 0.02011034

141 -1.45443942 ± 0.01394202 -1 43855778 ± 0.03063944

161 -1.49280680 ± 0.01751208 -1 45537928 ± 0.04689971

181 -1.48811749 ± 0.02147555 -1 47362783 ± 0.06736764

201 -1.47144283 ± 0.02575923 -1 56041898 ± 0.09701034

Таблица 3. Параметры Л и а асимптотики освещенности, вычисленные методом параметрических временных производных для а3 = 0.95, ас = 0.05

Без учета поляризации С учетом поляризации

г Л а Л а

101 -0.0644246 -1.4568894 -0.0645330 -1.4678340

121 -0.0621121 -1.4655645 -0.0620577 -1.4589837

141 -0.0604434 -1.4725182 -0.0605189 -1.4831649

161 -0.0593987 -1.5131866 -0.0596728 -1.5573155

181 -0.0583921 -1.5189743 -0.0587141 -1.5772449

Л а

рических временных производных для а3 = 1, ас = 0.

Без учета поляризации С учетом поляризации

г Л а Л а

101 -0.0142817 -1.4424541 -0.0142357 -1.4378097

121 -0.0118947 -1.4392532 -0.0120252 -1.4550459

141 -0.0105072 -1.4815099 -0.0102723 -1.4484004

161 -0.0091509 -1.4732894 -0.0093292 -1.5019978

181 -0.0080240 -1.4523529 -0.0081027 -1.4665904

вании с поляризацией. Здесь значения параметров Л и а асимптотики освещенности в скалярном и векторном вариантах также совпадают с точностью до статистической погрешности.

Отметим, что полученные в расчетах значения а близки к значению -3/2, полученному в работе [14] аналитическими выкладками при соответствующих асимптотических предположениях.

4.2. Исследование дисперсий оценок

Как указывалось в разд. 2, полученные, например, из неравенств (2.4)-(2.6) теоретические выводы о конечности дисперсий оценок функционалов и их производных являются основанием для применимости метода Монте-Карло. Однако, являясь конечными, значения этих дисперсий в задачах с поляризацией могут быть очень велики, что приводит к большим затратам вычислительных ресурсов для решения конкретных задач. Для изучения этой проблемы проведены специальные численные эксперименты по исследованию поведения дисперсий векторных оценок при критических значениях коэффициента поглощения в среде.

В табл. 5 и 6 приведены вычисленные значения средних квадратов оценок функционалов, используемых для нахождения параметров временной асимптотики в момент времени Ь = 96 и при коэффициенте поглощения среды ас = 0.1. Индексы 5 и V у функционалов в таблицах относятся к скалярному и векторному вариантам соответственно. В табл. 5 приведены результаты при конечной дисперсии векторных функционалов, в табл. 6 — при расходящейся. Дисперсии скалярных функционалов в обеих таблицах — конечны.

Таблица 5. Средние квадраты оценок функционалов для ас = 0.1, Ь = 96. Аэрозольное рассеяние

Число траекторий ЕЩ2 ад о2 ЕЩ2 Е(.1и ')2

104 0.2079е—11 0.3871е—11 0.2026е—11 0.1989е—10

105 0.2077е—11 0.4108е—11 0.2325е—11 0.6383е—11

106 0.2074е—11 0.3931е—11 0.2215е—11 0.5691е—11

ю7 0.2025е—11 0.3964е—11 0.2062е—11 0.5594е—11

ю8 0.2023е—11 0.3837е—11 0.2050е—11 0.5303е—11

ю9 0.2021е—11 0.3793е—11 0.2047е—11 0.5285е—11

ю10 0.2019е—11 0.3787е—11 0.2042е—11 0.5278е—11

Таблица 6. Средние квадраты оценок функционалов для ис = 0.1, Ь = 96. Молекулярное рассеяние

Число траекторий ЕЩ2 Е(,]3'У ЕЩ2 Е(.1и ')2

104 0.4970е—12 0.7192е—12 0.5755е—12 0.5202е—12

105 0.6205е—12 0.9302е—12 0.2773е—11 0.2287е—11

106 0.6075е—12 0.1134е—12 0.2960е—10 0.1245е—11

10' 0.5939е—12 0.1079е—12 0.2659е—09 0.8205е—10

108 0.5935е—12 0.1049е—12 0.1922е—08 0.8695е—08

ю9 0.5938е—12 0.1052е—12 0.9523е—08 0.4131е—10

Таблица 7. Средние значения и дисперсии векторной оценки плотности столкновений. Молекулярное рассеяние. Коэффициент поглощения 0.05, точное значение функционала 20

Число Е3 Ш

траекторий

104 17.9073 7 728

105 19.108 17 129

106 19.078 64 638

107 19.55 565 938

108 19.758 1 202 333

Расчеты показали, что при значениях ас, больших критического (см. разд. 2), статистические оценки сходятся к их математическим ожиданиям, а статистические дисперсии остаются постоянными при увеличении числа траекторий. При ас, близких к критическому и меньше, наблюдается рост дисперсий с увеличением количества траекторий, хотя оценка вектора Стокса может оставаться достаточно точной. Так, в табл. 7 приведены результаты расчетов среднего значения и дисперсии векторной оценки плотности столкновений в пространственно-однородной задаче, в которой известно точное решение, но дисперсия оценки метода Монте-Карло является бесконечной. Из таблицы видно, что статистическая дисперсия сильно растет, а оценка — быстро сходится. Однако в более сложной задаче об оценке параметров асимптотики при молекулярной матрице рассеяния (см. табл. 6) в расчетах наблюдался не только рост дисперсии, но и очень слабая сходимость оценки к конечному математическому ожиданию, не позволяющая получить результат за приемлемое время.

Список литературы

[1] Дэвисон Б. Теория переноса нейтронов. М.: Атомиздат, 1960.

[2] Михайлов Г.А., Трачева Н.В., Ухинов С.А. Исследование временной асимптотики интенсивности поляризованного излучения методом Монте-Карло // Журн. вычисл. математики и мат. физики. 2007. Т. 47, № 7. С. 1264-1275.

[3] Михайлов Г.А. Весовые алгоритмы статистического моделирования. Новосибирск: Изд-во ИВМиМГ СО РАН, 2003.

[4] Метод Монте-Карло в атмосферной оптике / Г.И. Марчук, Г.А. Михайлов, М.А. Наза-ралиев и др. Новосибирск: Наука, 1976.

[5] Михайлов Г.А. Весовые методы Монте-Карло. Новосибирск: Изд-во СО РАН, 2000.

[6] Михайлов Г.А., Войтишек А.В. Численное статистическое моделирование. М.: Учебно-изд. центр "Академия", 2006.

[7] Ухииов С.А., Юрков Д.И. Оценки методов Монте-Карло для параметрических производных поляризованного излучения // Сиб. журн. вычисл. математики. 2002. Т. 5, № 1. С. 40-56.

[8] Михайлов Г.А., Ухинов С.А., Чимаева А.С. Дисперсия стандартной векторной оценки метода Монте-Карло в теории перноса поляризованного излучения // Журн. вычисл. математики и мат. физики. 2006. Т. 46, № 11. С. 2199-2212.

[9] Лотов а Г.З., Михайлов Г.А. Новые методы Монте-Карло для решения нестационарных задач теории переноса излучения // Журн. вычисл. математики и мат. физики. 2002. Т. 42, № 4. С. 569-579.

[10] MVAPICH: MPI over InfiniBand, http://mvapich.cse.ohio-state.edu/

[11] Сибирский суперкомпьютерный центр, http://www2.sscc.ru

[12] Марченко М.А., Михайлов Г.А. Распределенные вычисления по методу Монте-Карло/ / Автоматика и телемеханика. 2007. № 5. С. 157-170.

[13] Intel® Math Kernel Library. Vector Statistical Library Notes, http: / / download.intel.com / software/products / mkl/docs / vslnotes.pdf

[14] Романова Л.М. Предельные случаи функции распределения по пробегам фотонов, выходящих из толстого светорассеивающего слоя // Изв. АН СССР. Сер. ФАО. 1965. Т. 1, № 6. С. 599-606.

Поступила в редакцию 20 февраля 2008 г.