О двух численных схемах метода Монте-Карло для решения уравнения Больцмана Текст научной статьи по специальности «Математика»

О ДВУХ ЧИСЛЕННЫХ СХЕМАХ МЕТОДА МОНТЕ-КАРЛО ДЛЯ РЕШЕНИЯ УРАВНЕНИЯ БОЛЬЦМАНА*

Н. М. Москалева

С.-Петербургский государственный университет, канд. физ.-мат. наук, [email protected]

Известно, что для решения ряда нелинейных кинетических уравнений динамики разреженных газов могут быть использованы столкновительные процессы [1, 2] и ветвящиеся случайные процессы [1, 3]. Моделируя движения молекул в соответствии со «столкновительным» процессом («прямая» схема), можно оценить решение уравнения сразу в заданной области и представить его в виде гистограммы распределения молекул. Но для надежного расчета больцмановской функции распределения ](Ь,£,х) в области высоких скоростей эта схема, как и метод прямого статистического моделирования [4, 5], требует значительных вычислительных ресурсов. Численная схема, основанная на моделировании ветвящихся процессов, может быть интерпретирована как «сопряженная» схема метода Монте-Карло. «Сопряженная» схема позволяет оценивать решение в заданной точке и ее преимущества проявляются при вычислении ](Ь, £, х) для больших значений скорости £ (на «хвостах») и областях физического пространства с малой вероятностью попадания в них молекул. Большинство задач решается методом прямого статистического моделирования, но для «тонких» задач кинетической теории газов требуется специальная техника, одним из вариантов которой является предлагаемая «сопряженная» схема.

Основной результат статьи состоит в модификации «прямой» и «сопряженной» схем и их совместном использовании в разных подобластях течения газа. Рассматривается комбинированная «прямая» схема, в которой движение частицы «наверх» (от начального момента до момента столкновения т) моделируется в соответствии с «прямой» схемой, а «обращенное» движение партнера по столкновению «вниз» (от момента т до начального момента) моделируется с помощью «сопряженной» схемы.

В [3] обсуждались трудности, возникающие при использовании численных схем, в основе которых лежит представление итерационного решения уравнения в виде математического ожидания несмещенной оценки на траекториях соответствующего случайного процесса. «Мажорантное» условие [1], обеспечивающее существование стандартной оценки, оказывается достаточно жестким для уравнений со знакопеременными ядрами. В частности, это условие ограничивает возможности применения метода для решения уравнения Больцмана. В связи с этим построение новых оценок с более слабыми «мажорантными» условиями остается актуальной задачей. Ослабление мажорантных условий может быть достигнуто путем «суммирования» определенным образом сгруппированных знакопеременных членов ряда, которым представлено итерационное решение уравнения. В методе Монте-Карло это приводит к необходимости построения специального случайного процесса, реализации которого определяют требуемый порядок объединения и суммирования членов ряда Неймана. Этот прием использован при конструировании оценки в [3]. Там же отмечено, что алгоритм не слишком прост и

* Работа выполнена при финансовой поддержке РФФИ (грант №08-01-00194-а).

© Н. М. Москалева, 2010

поэтому его встраивание в комбинированную «прямую» схему существенно усложнило бы вычислительную процедуру. В данной статье предлагается более простой способ построения оценки п для «сопряженной» схемы и оценки £ для модификации «прямой» схемы. Для обеих оценок способ объединения и суммирования слагаемых ряда Неймана одинаков, но отличен от [3]. Простота моделирования процессов и совместное использование легко реализуемых оценок п и £ позволили существенно повысить эффективность алгоритма.

Постановка задачи

Рассмотрим задачу Коши. Требуется найти функцию / : [0,Т] х К3 х К3 ^ [0, то), удовлетворяющую уравнению Больцмана В/ = Я3 /П^2 (/' /' 1 — //)Б(6, у)^\С6Се,

Ь € [0,Т], и начальным условиям / = /о,Ь = 0. Для степенных потенциалов Кт-Э с обрезанием по углу Б(в,ь) = Vй в(6), где V = |£ — ^х|, а = 1 — 4/в [6] .

Обозначим через р(Ь, £, х) отклонение больцмановской функции распределения / от глобального максвелловского распределения /т = (2п)3/2 ехр(—£2/2). Полагая в исходной задаче / = /т(1 + р), известным образом (напр., [6]) получаем уравнение

Вр + ри(/т ) = V (р,р), (1)

где у(/т) = /Й3101212 /т(£1)Б(6^)й£1в,6в,е, V(р,р) = ////т(£1)Б(6^)(р' + р\ + р'р' 1 — Р — Р1)С£1С6Се.

Проведя формальное интегрирование в (1), получим интегральное уравнение

рЬ с рп/2 р2п

р = р0в-1'(1т)Ь + )Ст / /т(£1 )Б(6^)Ф(р,р)^1 С6Се, (2)

■)о Зяз 3 о Зо

где Ф(Р,Р) = Ф1(Р,Р) + Ф2(P,P), Ф1(Р,Р) = Р(т,£,хт ~) + Р(т,£1,хт ~) + р(т,£,хт)х P(т, £'l, хт ), Ф2(P, Р) = — P(т, £1,хт ) — P(T, ^ хт )P(T, £1, хт ).

Чтобы подчеркнуть связь уравнения (2) со случайным процессом, запишем его в виде

P(t, ^ х) = Ро(^ х — £t)g(t, £)+

рЬ с гп/2 р2п

+ (1 — д(Ь,£)) ЛтРт (Ь,т,£) с£1Р1&&) зт(6) РзФ(р,рЖ

о я3 о о

Здесь д(Ь,£) = ехр(—^(/т)Ь), и(/т) = /яз /т (£1 )|£ — £11а^1 Л12 в(6)С6Се = 2пво^т(£), Рт(Ь,т,£) = у(/т)ехр(—у(/т)(Ь — т))/(1 — д(Ь,£)), Р^,^) = 1£ —

Ыа/т(£1)/"т(0, Р2(6) = в(6)/во, Рз = 1/2п, Рт,РьР2,Рз — плотности распределения т, £1,6, е соответственно, д(Ь, £) имеет смысл вероятности поглощения (вероятность выхода случайного процесса за пределы 1).

В дальнейшем при конструировании оценок нам понадобятся следующие функции: Фо(ро(г1),ро(г2)) = Ф] (р0(т1),р°(т2))/д(п)д(т2), 3 = 1, 2, п = (т,£,хь), I =1, 2, р0(т) = Pо(r)g(r),

1 + р0(п)

д(т1)

если 3 = 1, I = 1, 2;

Рз(п) = { -(1 + ‘/(п) -071

------7—N---------------------1 если 3=2, 1 = 1;

д(т1)

ро(п), если 3 = 2, I = 2.

1. «Сопряженная» схема. Оценка п- Следуя общей теории метода [1, 7], будем искать решение в виде математического ожидания оценки п на траекториях ш ветвящегося марковского процесса.

Траектория модельной молекулы является реализацией ветвящегося процесса с начальной плотностью Ро(£, х) и переходной плотностью:

где (1 — д(Ь,£)) —вероятность столкновения частиц на промежутке времени [0,Ь]. В нашем случае Ро(£,х) является ^-мерой.

Моделирование ш осуществляется следующим образом.

Траектория молекулы (частицы) начинается в момент времени Ь из точки (£,х) € К3 х К3. Далее прослеживается ее «обращенное» движение от момента Ь до начального момента. С вероятностью д(Ь, £) частица поглощается (г = 0) в точке х—£Ь, и траектория обрывается. С вероятностью 1 — д(Ь, £) частица переходит в точку хт = х — £(Ь — т), где происходит столкновение с другой частицей, имеющей скорость после столкновения £1. Момент времени т € [0,Ь], £1 и параметры столкновения 6 и е распределены в соответствии с плотностями т,Р1,Р2,Р3. В результате столкновения с вероятностью г] (3 = 1, 2) рождается пара частиц типа 3, образующих первое поколение (к = 1). При этом типу 3 = 1 соответствует пара частиц со скоростями до столкновения £' и £' 1, типу 3 = 2 соответствует пара частиц со скоростями £ и £1 («фиктивное» столкновение). Каждая из рожденных частиц ведет себя независимо от других и так же, как начальная частица, поглощаясь (г = 0) или порождая следующее поколение (г = 2). Процесс заканчивается при поглощении всех частиц.

Введем обозначения:

к — номер поколения, к = 1,...,К. Поколение к образуют пары частиц с номерами шь = 1,..., Ми: I — номер частицы в паре (I = 1, 2); если номер одной из частиц пары равен I, то номер другой частицы равен ц = 3 — I. Каждая из частиц I пары шь, находясь в состоянии т^шь), может породить число потомков г^шь) = 0, 2.

Приведем алгоритм построения оценки п.

1. В исходной точке т = (Ь, £, х) траектории положить Ао ^ /т.

2. Далее возможны случаи:

Исходная частица поглотилась в точке т = (0,£,х — £ * Ь), т. е. К = 0.

Положить Ао ^ Ао * ро(т). Перейти на п. 3.

Число поколений К > 1.

Вычислить вклад в оценку единственной пары частиц Ш1 3-го типа 1-го поколения.

Положить А1 ^ Ао * Ш] (ш1), к ^ к +1.

Аналогично вычисляется вклад в оценку всех пар частиц к-го поколения (к = 2...К). Каждая пара шь = 1, ...,Мь вносит в оценку «вес» Ш](шь). мк

Ак ^ Ак-1 * П (шк).

т-к

3. Положить п ^ А к.

Алгоритм закончен.

Р(Ь, т, £, £1, 6, е) = (1 — д(Ь, £))Рт(Ь, т, £)Р1 (£, £1)Р2(6)Р3,

Описанный ветвящийся процесс связан с представлением итерационного решения р(т) в виде

р(т) = ^ Б~/(r),

■уЕГ

где (т) —сумма членов итерационного ряда, которым сопоставлены траектории ш, имеющие одинаковую плотность вероятности. Траектория ш имеет структуру дерева 7 из множества Г (понятия «дерева» и «поколение частиц» см. в [7]), из каждой вершины которого исходит либо 0 (поглощение), либо 2 ветви (столкновение частиц).

Можно записать (т) в виде интеграла по траектории ш:

(т) = J А1 (ш)С^,1 (ш),

где мера определена с помощью начальной и переходной вероятностей, описывающих ветвящийся процесс. Интегрирование ведется по точкам траектории ш, в которых произошли столкновения частиц. Функция А7(ш) конструктивно определена выше.

Можно доказать, что Еп = р(т), если выполняется мажорантное условие

J 1А7(ш)13,л1 (ш) < <х>.

Приведем пример, поясняющий правило вычисления А7(ш). Запишем первую итерацию уравнения Больцмана в виде

Р{1)(Ь,£,х) = р0(Ь,£,х — £Ь)+

+ 1Р (Ь,т, £, £1,6,е)д(т,£')д(т,£' 1)$1(ро(£' ,хт — т£' ),ро(£\,хт — т£ 1))СтС£1С6Се+

+ IР(Ь, т,£, £1,6,е)д(т,£)д(т, £1)Ф°(ро(£,хт — т£),ро(£1,хт — т£1))СтС£1 С6Се.

Рассмотрим второе слагаемое, представляющее собой сумму трех членов итерационного ряда, которым соответствуют траектории с плотностью вероятности (1 — д(Ь, £))Р(Ь, Т, £, £1, 6, е)д(т, £')д(т, £' 1). На этой траектории

А(ш) = ро(£\ хт — т£' )g-1(т, £/) + ро(£'l, хт — т£' l)g-1(т, £/1)+

+ ро(£/,хт — т£/)Ро(£/1 ,хт — т£/ 1) = Фо(ро,ро).

2. Комбинированная «прямая» схема. Уравнению Больцмана можно сопоставить процесс с переходной плотностью Р (Ь, £,х ^ т, £/, хт; т,£/ 1 ,хт), который описывает процесс не «размножения», а «слипания» частиц при движении «наверх». В. В. Некруткин в [1, 2] построил и подробно исследовал «столкновительные» процессы для решения нелинейных кинетических уравнений. Там же приведен алгоритм построения несмещенной оценки. В предлагаемой комбинированной «прямой» схеме построена оценка С, для которой члены ряда Неймана группируются и суммируются таким же образом, как для п. Соответствующий процесс моделируется следующим образом.

В начальный момент времени рождается частица, фазовая координата (х1, £1) которой имеет распределение Ро (х1 ,£1). Время свободного движения т имеет показательное

распределение с параметром у(/т)(£). Если т > Ь, то траектория обрывается, иначе в момент времени т в точке хт происходит столкновение с другой частицей. Скорость партнера по столкновению £д распределена в соответствии с плотностью Р1(£ ,£д). С вероятностью г] имеет место столкновение типа 3 (3 = 1, 2). Если 3 = 1, то разыгрываются параметры столкновения 6 и е согласно Р3, Р4 и вычисляется скорость исходной частицы после столкновения. Эта пара частиц образует первое поколение. Дальнейшее движение исходной частицы из точки хт происходит по выше описанному правилу, а «обращенная» траектория партнера по столкновению из точки (т, £д ,хт) до начального момента времени моделируется по «сопряженной» схеме.

Траектории этого процесса имеют структуру деревьев 7, в каждом поколении которых рассматриваются только две частицы: исходная и партнер по столкновению Ць.

Здесь —номер исходной частицы в к-м поколении, = 1, 2. Значение указывает на принадлежность исходной частицы к левой (1ь = 1) или правой (1ь = 2) ветви и разыгрывается при каждом столкновении частиц, ць — номер партнера по столкновению (ць = 3 — 1ь). Координата частицы т(1ь) = (ть,£,хь), Б](т(1ь)) = Б](1ь) и

£ )) = в^.

Траектории частиц ць имеют структуру деревьев , описанных выше в «сопряженной» схеме.

Приведем описание алгоритма построения оценки £.

1. В начальный момент рождается частица с номером I = 1,

к <— 1.

Получаем реализацию случайной величины I с распределением

' 1 2 \ _ ть <— 0.

1/2 1/2

Фазовая координата частицы (хі,£і) распределена с плотностью Ро(£і,хі). Полагаем С <— Іт(£і)роШ,хі)/РоШ,хі), хк <— хі.

2 . к <— к + 1.

Получаем реализацию т из показательного распределения с параметром V(/т)(£і). Тк <- Тк-1 + т.

Если Тк > і, то переход на п. 7.

3. Моделирование столкновения с частицей ц в точке хк: хк <— хк-і + т£і,

Получаем реализацию случайной величины I с распределением

1/2 1/2 )’ Тк ^-°.

Скорость распределена согласно Р2(£і ,£д).

Получаем реализацию случайной величины 3 с распределением

1 2 1 Тк _ 0,.

*1 ) С — С/*і.

4. С вероятностью д(Тк, £д, хк) частица ц поглощается в точке Тд = (0, , хк — Тк * £д).

Полагаем С <— С * ,

шк =\ ф0(ті,тч), если к = 1;

і [ ^ (тд), если к = 1.

Переход на п. 6.

5. Если частица ц не поглотилась, то, применяя «сопряженную» схему, вычислем в точке гч = (тк, £ч,Хк) оценку вд интеграла Яд = Я(гд).

Полагаем С <— С * ^ К,

]¥к _ ) ГМ) * 5',, если к = I;

3 [ 5'д, если к ^ 1.

6. Если 3 = 1, вычислить скорость £1 после столкновения. Переход на п. 2.

7. Конец моделирования траектории.

3. Результаты моделирования. На примере решения нелинейной задачи исследовались работоспособность модифицированной «прямой» схемы и возможности новой «сопряженной» схемы для вычисления с высокой точностью функции распределения на «хвостах». Решалась нелинейная задача о релаксации максвелловского газа с начальными условиями специального вида. Результаты расчетов сравнивались с точным БКШ-решением [8, 9]. При моделировании предполагалась зависимость функции распределения от вектора скоростей, но в виду изотропности задачи в окончательных результатах (таблицы и графики) отражена зависимость от модуля скорости.

Целесообразно применять эти две схемы в разных областях течения газа. Расчеты течения в области и € [0; 3.5](и = |£|) наиболее просто и быстро реализуются с помощью комбинированной «прямой» схемы, «хвосты» (и € [3.5; 5]) функции распределения

/(“)

и

------"прям." сх ------ТОЧН. реш.

Рис. 1. Функция распределения /(и) (4 = 1).

Таблица 1.

и Выборочное Точное Выборочное Точное

среднее /(и), значение /Ь(и), среднее /(и), значение /Ъ(и),

* = 5 t= 5 1 = 10 1 = 10

0.1 5.747 * КГ- 5.752 * КГ2 6.237* КГ2 6.237* КГ2

0.9 4.080 * КГ2 4.078 * КГ2 4.210 * кг2 4.209 * КГ2

1.5 2.104 * КГ2 2.103 * КГ2 2.066 * КГ2 2.066 * КГ2

1.9 1.087 * КГ2 1.088 * КГ2 1.051 * 10-2 1.051 * 10-2

2.5 0.285 * КГ2 0.285 * КГ2 0.280 * КГ2 0.280 * КГ2

3.1 0.480 * КГ2 0.480 * КГ2 0.512 * КГ2 0.514 * КГ2

3.5 0.113 * КГ2 0.114 * КГ2 0.134 * КГ2 0.134 * КГ2

4.1 0.921 * КГ2 0.912 * КГ2 0.131 * КГ2 0.130 * КГ2

4.5 0.136 * КГ2 0.131 * КГ2 0.225 * 10-2 0.222 * 10-2

4.9 0.162 * КГ2 0.155 * КГ2 0.302 * КГ2 0.317* КГ2

/(«)

------- "прям." сх -----------"сопр." сх

Рис. 2. Относительная погрешность вычисления f (и) на «хвостах», 4 = 1.

Таблица 2.

4 Выбороч. знач. Мд Выбороч. знач. М| Точное знач. Мб Выбороч. знач. Мд Выбороч. знач. М| Точное знач. Ма

1 7.5 10 74.25 100.72 103.02 74.61 100.95 103.05 74.47 101.16 103.28 518.54 876.21 913.61 524.34 880.40 914.09 525.63 881.08 915.69

%

3,0 г

5 10 15

t

------"сопр." схема--------"прям." схема

Рис. 3. Относительная погрешность вычисления М%.

молекул оцениваются с помощью «сопряженной» схемы. Такое ограничение диапазона скоростей вызвано только тем, что для данной задачи в [3] при использовании менее эффективного алгоритма диапазон составил шесть тепловых скоростей.

На рис. 1 и табл. 1 представлены результаты расчетов /(и) при Ь = 1, 5,10. Зависимость относительной погрешности функции распределения от выбранной схемы расчетов при и € [3.5; 5] и Ь = 1 приведена на рис. 2. В табл. 2 даны значения моментов Мб(Ь) и М$(Ь) при Ь = 1, 7.5,10, вычисленные в области и € [3.5; 5] с использованием «прямой» схемы — МП, МП и «сопряженной» — МС, М§. Влияние на значения моментов Мб(Ь),М$(Ь) коррекции функции распределения на «хвостах» с помощью «сопряженной» схемы показано на рис. 3.

Таким образом, представленные результаты подтверждают эффективность совместного использования рассмотренных схем для получения достаточно точного решения в диапазоне до пяти тепловых скоростей при небольших затратах машинного времени.

Литература

1. Ермаков С. М., Некруткин В. В., Сипин А. С. Случайные процессы для решения классических уравнений математической физики. М., 1984. 205 с.

2. Некруткин В. В. О случайном процессе, решающем уравнения типа Больцмана // Методы Монте-Карло в вычислительной математике и математической физике. Новосибирск, 1976. С. 155-162.

3. Ермаков С. М., Москалева Н. М. О моделирующих алгоритмах высокой точности // ДАН. 1995. Т. 343. №4. С. 475-477.

4. Абрамов А. А., Гимельшейн С. Ф., Иванов М. С., Макашев Н. К. Высокоскоростные «хвосты» функций распределения молекул для одномерных неизотермических течений // МЖГ. 1999. №2. С. 159-169.

5. Иванов М. С., Рогазинский С. В. Экономичные схемы прямого статистического моделирования течений разреженного газа // Мат. моделирование. 1989. Т. 1. №7. С. 130-145.

6. Маслова Н. Б. Теоремы о разрешимости нелинейного уравнения Больцмана // Черчи-ньяни К. Теория и приложения уравнения Больцмана. М.: Мир, 1978. С. 495.

7. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1975. С. 471.

8. Бобылев А. В. О точных решениях уравнения Больцмана // ДАН. 1975. Т. 225. №6. С. 1296-1299.

9. Krook M, Wu T. T. Phys. Fluids. 1977. Vol. 20. N 10. P. 1589-1595.

Статья поступила в редакцию 15 июня 2010 г.