ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 1, с. 10-21
УДК 519.676
ИССЛЕДОВАНИЕ И УЛУЧШЕНИЕ СМЕЩЕННЫХ ОЦЕНОК
МЕТОДА МОНТЕ-КАРЛО1
© 2015 г. Г. З. Лотова, Г. А. Михайлов
(630090 Новосибирск, пр-т Акад. Лаврентьева, 6, ИВМиМГСО РАН; 630090 Новосибирск, ул. Пирогова 2, НГУ) e-mail: lot@osmf.sscc.ru, gam@osmf.sscc.ru Поступила в редакцию 04.07.2014 г.
Численное статистическое моделирование свободного пробега частицы для столкновитель-ной модели процесса переноса с учетом ускорения внешним силовым полем реализуется шагами по времени; в работе построена новая конструктивная оценка соответствующей детерминированной относительной погрешности, которая позволяет выбрать подходящую величину шага. Стандартные статистические "локальные оценки" плотности потока частиц являются смещенными вследствие зануления вкладов от столкновений в "локальном шаре" малого радиуса для ограничения дисперсии; в работе представлены практически эффективные оценки соответствующей относительной погрешности. Дополнительно осуществлена равномерная оптимизация функциональной оценки плотности распределения частиц типа гистограммы в предположении "пуассоновости" соответствующего статистического ансамбля. Оказалось, что в оптимальных (по трудоемкости) вариантах рассмотренных алгоритмов детерминированная погрешность близка к статистической. Библ. 11. Табл. 3.
Ключевые слова: статистическое моделирование, свободный пробег, шаг по времени, трудоемкость, оценка по столкновениям, двойная локальная оценка, пуассоновский ансамбль, гистограмма.
DOI: 10.7868/S0044466915010159
1. ВВЕДЕНИЕ
Принято полагать, что метод Монте-Карло сравнительно эффективен для оценки ограниченного числа линейных функционалов от решения, которые представляются интегралами, вообще говоря, бесконечной кратности (влияние аргументов обычно экспоненциально убывает, что и дает возможность построения достаточно эффективной статистической оценки). На примере оценок одного функционала J продемонстрируем основной подход для минимизации трудоемкости, т.е. времени ЭВМ или среднего числа операций, которое необходимо для достижения заданной погрешности 5. Будем считать, что погрешность статистической оценки измеряется величиной среднеквадратичного отклонения, так как асимптотически при увеличении числа статистических испытаний, как правило, используемые на практике доверительные интервалы пропорциональны этой величине (см. [1]). Сначала рассматриваются несмещенные оценки, т.е. предполагается, что Е2, = J.
Итак, пусть дано семейство случайных величин {£,} с конечными дисперсиями Щ и таких, что Е2, = J. На основе N независимых статистических испытаний строим случайные оценки вида
N
:(0
Е
zn —1 N , E Zn —J, D Zn — dN .
1) Работа выполнена при финансовой поддержке Программы фундаментальных исследований Президиума РАН № 1 по стратегическим направлениям развития науки на 2014 г. "Фундаментальные проблемы математического моделирования"; грантов РФФИ (коды проектов 12-01-00727, 12-01-00034, 12-05-00169) и гранта НШ-5111.2014.1.
Обозначим через t среднюю трудоемкость реализации случайной величины т.е. средние затраты на проведение одного случайного испытания. Минимизация трудоемкости рассматриваемых оценок эквивалентна следующей задаче на минимум:
^ = Nt ^ min, ^ = S2, N
или tD£, ^ min (см. [1]). Таким образом, в классе несмещенных оценок с конечной дисперсией минимизация трудоемкости не зависит от требуемой погрешности 5, а также инвариантна относительно N.
Рассмотрим теперь однопараметрическое семейство смещенных оценок с трудоемкостью t(s), дисперсией d(s) и смещением порядка sa, т.е. |E£,e — J\2 ~ ее2а. Средний квадрат отклонения здесь определяется соотношением
E(£,e - J)2 = D£,e +(J - E^)2.
Следовательно,
E(Ze,n- J)2 - ^ + Cs2a. (1.1)
Таким образом, минимизация трудоемкости здесь эквивалентна следующей задаче на условный минимум:
Ss = Nt(e) ^ min, ^ + сs2a = 52
N,e
N
или
t(e) d(e) 2a o2
- - - - ^ min, ce <5 .
S2 2a
- ce E
В результате решения этой задачи получаются оптимальные значения N и e как функции 5, что дает некоторую связь между N и e. Таким образом, используемое иногда на практике уравнивание слагаемых в правой части выражения (1.1) не дает оптимального варианта оценки. Однако, как правило, способ уравнивания дает оптимальный порядок трудоемкости при 5 —► 0; это, в частности, показывают результаты настоящей работы.
Наибольшее практическое значение имеют аналитические и численные исследования смещения cea с помощью зависимых испытаний. В разд. 2 настоящей работы приведены результаты таких исследований для выбора шага по времени в статистическом моделировании переноса заряженных частиц с учетом их взаимодействия со средой и ускорения внешним силовым полем.
Разд. 3 посвящен изучению смещения стандартных численно-статистических локальных оценок интенсивности потока частиц с конечной дисперсией.
Дополнительно в разд. 4 осуществляется равномерная оптимизация функциональной оценки плотности распределения частиц типа гистограммы в предположении пуассоновости соответствующего статистического ансамбля.
2. О ВЫБОРЕ ШАГА ПО ВРЕМЕНИ И ВЕРОЯТНОСТИ СТОЛКНОВЕНИЯ ПРИ ПОСТРОЕНИИ ТРАЕКТОРИЙ ЧАСТИЦ С УЧЕТОМ УСКОРЕНИЯ ВНЕШНИМ СИЛОВЫМ ПОЛЕМ
1. В настоящем разделе работы решается задача выбора временного шага А? в алгоритме численного статистического моделирования переноса заряженных частиц под влиянием взаимодействий (столкновений) с частицами среды и внешнего силового поля. Критерием выбора является подходящая детерминированная погрешность (смещение) статистических оценок изучаемых интегральных характеристик. Полученные теоретические результаты применяются далее для моделирования процесса развития электронных лавин в газах под влиянием взаимодействия со средой и внешнего электрического поля.
Статистическое моделирование основано на предположении о том, что при А? —»- 0 изучаемые функционалы сходятся к соответствующим величинам для физической модели переноса. В рамках этой модели случайные свободные пробеги частиц между столкновениями реализуются по линиям, удовлетворяющим уравнению движения с учетом скорости и ускорения (см., на-
пример, [2]). Длина свободного пробега имеет неоднородное (вследствие ускорения) экспоненциальное распределение, которое определяет пуассоновский точечный поток столкновений. Интенсивность соответствующего точечного пуассоновского потока на оси времени равна 'к(Т) = ст(7)^7), V = |У(2)|. Здесь Т — энергия частицы, а(Т) — "сечение" взаимодействия (см. [2], [3]), V = (Ух, Vy, V) — векторная скорость частицы. Для последующих выкладок необходима оценка величины a = (А,(Т')Д(Т))т, где Т' — энергия электрона непосредственно после взаимодействия, (...)т обозначает осреднение по типу взаимодействия т. При этом, вследствие правила объединения пуассоновских потоков, величина Х(Т') получается суммированием вкладов от всех типов взаимодействия; в случае ветвления траектории суммируются соответствующие вклады.
Траектории квантов моделируются шагами At по времени I. В результате очередного /-го шага частица с энергией Т1 -1 перемещается из точки г; _ 1 = (х; _ 1, У1 -1, -1) в точку г;, причем ее координаты и скорость V изменяются следующим образом (см., например, [2]):
г,. = г,.-! + У(-Аг + у;(-Ар-2, V, = V, + V'Аг, (2.1)
где выражение для г; получается интегрированием уравнения движения с постоянным ускорением V', что уточняет стандартную схему Эйлера. В стандартном алгоритме в конце каждого временного шага реализуется столкновение с вероятностью (см. [2])
Р = 1 - ехр--а-Т._1)А/), (2.2)
где А1 = |г; — гг -1|. Затем разыгрывается тип т взаимодействия в соответствии с заданными сечениями, которые для простоты предполагаются независимыми от г.
Можно полагать, что детерминированная погрешность оценок определяется тем, что в интервале длины А1 реализуется v1(Аt) = v1(Аt, Т) столкновений, причем Р^(А0 = 0) = 1 — Р^(А0 = 1), а в физической модели реализуется v(Аt) столкновений с Р^(А?) > 1) > 0. Следовательно, погрешность оценки связана с "недоучетом" столкновений. Далее все оценки будут осуществляться с точностью до величин 0((А|)2) включительно. Соответственно (2.2) имеем
Р^1 -Аг) = 1) - аА/- -аА/)2/2.
Далее, в лемме 1 будет рассматриваться траектория частицы с энергией Т в конкретном интервале времени; для удобства выкладок нулевой момент полагаем перенесенным в начало этого интервала. Символ Е далее обозначает математическое ожидание, т.е. осреднение по соответствующему статистическому ансамблю, а УУ; — скалярное произведение векторов V и V'.
Лемма 1. При условии дифференцируемости функции а( v(t)) имеет место оценка = ЕЧА0-ЕУ-А0 = Г аV + ¿Ма^) -А + 0-Аг).
Ev-Аг) V йг V у 2
Доказательство. С точностью до 0((А1)2) включительно имеем
Ev-Аг) - Р-v-Аг) = 1) + 2Р^-Аг) = 2).
Используя представление
-аv)-г) = а V + -ау); г + о - г), 0 < г <А г,
получаем
Аг / г \ / Аг \
йг-
P-v-Аг) = 1)- + -ау)'г)ехр|-|-а V + -ау)'л)йл ехр -а |-а V + -а V)'л)йл I
0 0 г
Аг 2
- |-а V + -а V) 'г)-1 - аv г)-1 - aаv-А г - г)) йг -аvА г + --а V)' - -1 + а )-а V)2 )-А)-,
0
Аг с г \
а а у(Д г - г)
( J
О о
Р(у(Дг) = 2) - |(а V + (а V)'г) ехрI - |(а V + (а V)'л)йл
Аг
- |а (а V)2 (Д г - г) йг = а ) (Дг)2.
о
Используя выражения
(Д /)2 = ( уд г + у;(Дг)2 /2 )2 и д / - vД г( 1 + vv;v-2д г/2), далее получаем оценку
Е V - EVl - ^ Га (а V)2 + (а V)' - а УУ'),
2 V V ^
т.е. лемма 1 доказана.
Интересно, что модификация стандартного алгоритма с заменой (2.2) на
2
Р (VI (Дг) = 1) = аД/-р(-аД)-, |р|< 1, (2.3)
может улучшить оценки, так как справедлива следующая
Лемма 2. В условиях леммы 1 с учетом (2.3) имеет место оценка
Я(Дг; Т) = ((а + р - 1 )а V + - ^ Д + «(Дг).
Дополнительно отметим, что при использовании стандартной схемы Эйлера: г,- = г,- _ 1 + V,- _1Д?, в утверждениях лемм 1, 2 слагаемое УУ' V-2 отсутствует.
Интегральные характеристики процесса являются линейными функционалами от плотности потока частиц Ф(г, Т, 1) и, следовательно, плотности столкновений ф(г, Т, 1), так как выполняется соотношение ф(г, Т, 1) = а(г, Т)Ф(г, Т, 1).
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.