автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Решение некоторых задач математической физики методами Монте-Карло

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

На правах рукописи

КУЗНЕЦОВ Андрей Николаевич

РЕШЕНИЕ НЕКОТОРЫХ ЗАДАЧ МАТЕМАТИЧЕСКОЙ ФИЗИКИ МЕТОДАМИ МОНТЕ-КАРЛО

05.13.18 — Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

диссертации на соискание учёной степени кандидата физико-математических наук

2 6 ДПР 2012

Санкт-Петербург — 2012

005017959

Работа выполнена на кафедре прикладной математики факультета прикладной математики и компьютерных технологий Вологодского государственного педагогического университета

Научный руководитель:

кандидат физико-математических наук, доцент, СИПИН Александр Степанович

Официальные оппоненты:

доктор физико-математических наук, профессор МЕЛАС Вячеслав Борисович (Санкт-Петербургский государственный университет)

кандидат физико-математических наук, ведущий инженер-программист ТИМОФЕЕВ Константин Алексеевич (Санкт-Петербург, ЗАО «Транзас технологии»)

Ведущая организация: ФГБОУ ВПО «Московский государственный

университет имени М.В.Ломоносова» (факультет ВМК МГУ)

Защита состоится «Л» __ 2012 г. в А? часов на заседании совета Д 212.232.51 по защите докторских и кандидатских диссертаций при Санкт-Петербургском государственном университете по адресу: 198504, Санкт-Петербург, Петергоф, Университетский пр., 28, математико-механический факультет, ауд. 405.

С диссертацией можно ознакомиться в Научной библиотеке им. М. Горького Санкт-Петербургского государственного университета по адресу: 199034, Санкт-Петербург, Университетская наб., 7/9.

Автореферат разослан «{&» ^ууг/иа^к 2012 г.

Учёный секретарь

диссертационного совета „ Кривулин Н. К.

Общая характеристика работы

Актуальность темы

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

В работе проведён обзор существующих методов нахождения взаимных электростатических ёмкостей проводников и рассмотрена возможность их вычисления методами Монте-Карло.

Знание первого собственного числа для оператора Лапласа может потребоваться для оценки применимости различных алгоритмов решения краевых задач, например, при решении задачи Гельмгольца с помощью «блуждания по шарам и сферам». В работе рассмотрены алгоритмы нахождения оценок итераций оператора Грина для оператора Лапласа и. первого собственного числа для оператора Лапласа. Рассмотрена возможность использования оценок итераций оператора Грина для расширения области применимости методов Монте-Карло при решении уравнения Гельмгольца с помощью аналитического продолжения решения методом замены переменных. Цели работы

• Разработка алгоритмов вычисления взаимных электростатических емкостей методами Монте-Карло.

• Исследование оценок итераций оператора Грина с целью вычисления первого собственного числа оператора Лапласа и решения уравнения Гельмгольца.

• Реализация перечисленных алгоритмов для вычислений с использованием МР1-кластеров и графических процессоров общего назначения.

• Создание программного комплекса, упрощающего «обычную» и «параллельную» реализацию указанных программ.

Методика исследования

Методика исследования включает применение методов статистического моделирования к решению краевых задач математической физики. С помощью формул Грина краевые задачи сводятся к интегральному уравнению, которое затем решается методом Монте-Карло. При решении уравнения Гельм-гольца используется аналитическое продолжение решения методом замены переменных [1]. При нахождении первого собственного числа для краевых задач используются известные итерационные методы.Для реализации алгоритмов использованы язык программирования С++, реализация MPI MPICH-2, технология CUDA.

Основные результаты, выносимые на защиту

• Универсальные алгоритмы вычисления взаимных электростатических ёмкостей методами «блуждания по сферам» и «блуждания по сферам и полусферам».

• Алгоритмы для оценок итераций оператора Грина и вычисления первого собственного числа оператора Лапласа.

• Комплекс программ для реализации алгоритмов Монте-Карло. Реализация перечисленных алгоритмов.

Научная новизна

Разработан новый универсальный стохастический алгоритм вычисления взаимных электростатических емкостей с помощью «блуждания по сферам». Алгоритм не требует того, чтобы границы объектов были плоскими. В случае же, если все или часть границ являются плоскими, алгоритм «блуждания по сферам и полусферам» позволяет ускорить вычисления и повысить точность расчетов за счет построения несмещённой оценки. Представлены новые алгоритмы вычисления итераций оператора Грина. Практическая значимость

Разработанные статистические алгоритмы вычисления взаимных электростатических ёмкостей проводников могут быть использованы для оценки паразитных электростатических ёмкостей в больших и сверхбольших интегральных схемах и многослойных печатных платах. Проведены вычислительные эксперименты для алгоритмов, использующих методы Монте-Карло при решении задачи Дирихле для уравнения Гельмгольца. Рассмотрены различные методы вычисления итераций оператора Грина методами Монте-Карло. Проведены вычислительные эксперименты, позволяющие оценить возможность использования аналитического продолжения решения методом заме-

ны переменных. Представлены новые алгоритмы вычисления оценок первого собственного числа для оператора Лапласа с помощью оценок степеней оператора Грина. Разработан программный комплекс для реализации вычислительных задач, использующих методы Монте-Карло. Программный комплекс позволяет упростить реализацию и проведение вычислительных экспериментов как на выделенных компьютерах, так и на MPI-кластерс или графическом процессоре общего назначения. Апробация работы

Основные результаты обсуждались на семинарах кафедры прикладной математики в Вологодском государственном педагогическом университете, семинаре кафедры статистического моделирования математико-механичес-кого факультета СПбГУ и докладывались на

• пятой Всероссийской научной конференции с международным участием «Математическое моделирование и краевые задачи» (СамГТУ, Самара, 29-31 мая 2008 г.);

• II ежегодном смотре-сессии аспирантов и молодых учёных по отраслям наук. Секция «Математика. Информатика. Методика преподавания» (Вологда, 12-14 ноября 2008 г.);

• научной школе-конференции молодых исследователей «Математические идеи П. JI. Чебышева и их приложение к современным проблемам естествознания» (Обнинск, 14-18 мая 2011 г.).

Публикации

Материалы диссертации опубликованы в работах [А1-А4]. Из них работы [AI—A3] — в списке журналов, рекомендованных ВАК. Работы [AI-A3] выполнены в соавторстве. Соискателю в них принадлежат реализация задач, численные результаты, соавторам — общая постановка задач и верификация результатов.

Структура и объём диссертации

Диссертация состоит из введения, 4-х глав, заключения и списка литературы. Список литературы содержит 51 наименование. Общий объем работы составляет 147 страниц.

Содержание работы

Во введении обосновывается актуальность диссертационной работы, формулируются задачи исследований, и содержится краткий обзор результатов работы.

В первой главе рассматривается проблема нахождения взаимных электростатических ёмкостей.

В первом разделе приводится постановка задачи нахождения электростатической ёмкости уединённого проводника и взаимных электростатических ёмкостей нескольких проводников.

Рассмотрим основную задачу электростатики в обратной постановке: требуется найти функцию у?, удовлетворяющую уравнению Лапласа = О вне заданной системы проводников, обращающуюся в нуль на бесконечности, принимающую некоторые постоянные значения ^pi на поверхностях проводников Г* и удовлетворяющую интегральному соотношению на поверхностях проводников

где п — вектор внешней нормали к поверхности проводника. В [2] показано, что эта задача имеет единственное решение.

После нахождения констант ^ определение потенциала сводится к решению внешней задачи Дирихле для уравнения Лапласа. Потенциал существенно зависит от формы и расположения других проводников. Согласно [2], имеют место соотношения:

91 = Сц91 + С\2 (^2 - Ы + • • ■ + С1п{<Р,1 ~ VI); '

<72 = С21(¥>1 ~ ^г) + С22</>2 + . • - + С2п(</?„ - Уг);

>

Ча = Спу(ч> 1 - <л,) + С^уг ~ <Рп) + • • • + Спп<рп,, где и у?,- — соответственно заряд и потенциал г'-го проводника; С^ — взаимная ёмкость или ёмкостной коэффициент г-го проводника по отношению к :7-му (величина Си также называется собственным ёмкостным коэффициентом).

Взаимная ёмкость определяется как тот заряд, который должен быть сообщён г-му проводнику для того, чтобы все проводники кроме .7-го, имели нулевой потенциал, а проводник — потенциал 1. Соответственно, собственный ёмкостной коэффициент определяется как отношение заряда к потенциалу на этом проводнике при условии, что все остальные проводники заземлены. Из [2] известно, что матрица коэффициентов С,7 является симметричной, т.е. Сц = Ср.

После нахождения взаимных ёмкостей, данная система может быть решена любым удобным способом.

Do втором разделе приведено краткое описание существующих алгоритмов нахождения взаимных электростатических емкостей.

Аналитическое решение применимо лишь к небольшому количеству задач. Приведены известные решения для емкости заряженного шара, эллипсоида, системы из двух проводящих сфер и плоского конденсатора (системы из двух равных параллельных проводящих пластин). Данные решения использовались для проверки разработанных в диссертации алгоритмов.

Для реальных задач используются численные методы решения. В работе приведено описание основных методов. Их можно разделить на следующие группы:

• Алгоритмы «сопоставления е эталоном» (pattern matching), в которых емкость сложной системы проводников оценивается как сумма ёмкостей более простых систем.

• Алгоритмы, использующие метод граничных элементов. В этом случае поверхность проводников разбивается на «плитки», для которых строится система интегральных уравнений, связывающих потенциал в центре каждой плитки с зарядами на остальных плитках. Далее полученная система решается относительно зарядов проекционными методами Галёркина или методом моментов. Затем находятся взаимные ёмкости. Наиболее известны алгоритмы FastCAP [3] и FFTCap [4], используемые, в том числе, для решения промышленных задач. Открытые реализации этих алгоритмов, были использованы при оценке корректности разработанных алгоритмов.

• Алгоритмы, использующие метод Монте-Карло. Для решения данной задачи используются алгоритмы «блуждания по квадратам» [5] и кубам [6]. В работе [7| приведен способ нахождения емкости выпуклого проводника с помощью случайного «блуждания по границе».

Алгоритм «блуждания по кубам» имеет промышленную реализацию QuickCAP, но он ориентирован на задачи, в которых проводники имеют плоские границы, кроме того, использование функции Грина для куба приводит к необходимости сё аппроксимации рядом Фурье, что вносит дополнительную ошибку в вычисления. В данной работе представлены алгоритмы «блуждания по сферам» и «блуждания по сферам и полусферам», не имеющие таких недостатков.

Третий раздел содержит теоретические и практические результаты данной диссертации, касающиеся нахождения электростатической ёмкости усди-

немного проводника с помощью «блуждания по сферам».

Рассмотрим решение задачи (1) для уединённого проводника. В этом случае для ёмкости С можно построить интегральное представление

С =--- { /2Л4 7(Ж, у)- 4Д2

где Я — радиус сферы с центром в нуле полностью содержащей проводник; ж и у — случайные точки, равномерно распределённые на поверхности сферы с центром в нуле радиуса 2/2 и Л соответственно.

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

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

Взаимные ёмкости проводников находим исходя из физического смысла, решая внешнюю задачу Дирихле. Для фиксированного г, на одном наборе траекторий можно находить все С^: пусть блуждание окончилось на к-м проводнике, тогда оценка для Су^д.. равна нулю, а для С^ — полученному весу, т.к. вблизи границы проводника ^ = 1.

В этом случае интегральное представление для ёмкости С имеет вид:

= / ¿2 / й{у'п) ^у + х) ^ Г 5Г

где Г — поверхность, содержащая внутри себя проводник г и отделяющая его от других, а <7г — площадь поверхности Г.

Если г-й проводник можно отделить от других сферой радиуса с центром в точке О,-, формула приобретает вид:

5«,(0,) Я

Опишем процесс блуждания для вычисления взаимных ёмкостей С^ при фиксированном г и = 1,2,..., п:

1) Для каждого проводника к выбирается оболочка Г д., отделяющая его от других. Построенная оболочка не должна касаться поверхности проводника.

2) Находится радиус сферы Я с центром в нуле, содержащей все объекты и их оболочки. Инициализируем оценку £ := 1.

3) Выбор точек ж и у. Точка ж распределена равномерно на поверхности оболочки, отделяющей данный проводник от других. Точка у равномерно распределена на сфере радиуса г = dist ^ Г*) с центром в точке ж. Вес перехода из ж в у вычисляется по формуле £ (у, тг), где п --- вектор внутренней нормали к поверхности оболочки в точке ж; a^ — площадь поверхности Г;.

4) Из точки у начинается процесс «блуждания по сферам»: у0 := у, п := 0.

5) Если после п шагов уп ф Кд, то уп+1 распределён на сфере Бц с плотностью

Р(Ж' У) =

£:=£Л/|уп| и осуществляется переход к шагу 7. В противном случае переходим к шагу 6.

6) Вектор уп+1 распределён равномерно на сфере 3Гп(уп) с центром в уп, где гп = (г)п, ^и Г^, а величина £ не изменяется.

п

7) Если точка уп+; не попала в г-окрестность и Г^, переход к шагу 5, иначе переходим к шагу 8.

8) В накопитель проводника, на котором остановилась траектория добавляется, вес, полученный на этой траектории, а в накопителях для остальных Сц лишь увеличивается счётчик траекторий. При этом счётчик не увеличивается в накопителях проводников, находящихся внутри г-го проводника, так как траектория лежащая снаружи не могла окончиться на этих проводниках.

9) В случае вложенных проводников, по окончании расчетов, ёмкость для внешнего проводника (г) вычисляется из полученных значений как Сц — ср. О;) где сумма берется по номерам проводников, содержащихся в г-м.

Для повышения точности оценки и ускорения работы, вместо «блуждания по сферам» можно использовать «блуждание по полусферам» [8]. В случае, если все границы проводников плоские, метод «блуждания по полусферам» дает несмещённую оценку. Моделирование осуществляется аналогичным образом.

Для модельных задач проведено сравнение оценок, полученных с помощью разработанного алгоритма, с аналитическими решениями и результата-

ми программ РавЮар и РРТСар.

Во второй главе рассматриваются алгоритмы вычисления оценок степеней оператора Грина для уравнения Лапласа и решение задачи Дирихле для уравнения Гельмгольца с их помощью. Представлены теоретические и практические результаты данной диссертации, касающиеся нахождения оценок степеней оператора Грина и применения метода аналитического продолжения решения к задаче Дирихле для уравнения Гельмгольца

Ди + Аи =-/;! ^

и|г = 0, /

где А — постоянное, возможно комплексное, значение.

В первом разделе рассмотрены алгоритмы решения задачи Дирихле для уравнения Гельмгольца методами «блуждания. по шарам и сферам» [9] и «блуждания по шарам» [10]. Эти алгоритмы не используют оценок степеней оператора Грина. Они реализованы для сравнения с разработанными нами методами.

Второй раздел содержит представление решения задачи Дирихле для уравнения Гельмгольца в виде суммы степеней оператора Грина. Приведено описание метода аналитического продолжения решения [1] посредством замены переменных.

Третий раздел посвящен вычислению оценок степеней оператора Грина для оператора Лапласа методами Монте-Карло. Рассмотрены как стандартные методы построения с помощью «блуждания по шарам и сферам», «блуждания по сферам» (через представление решения в виде суммы объёмного потенциала и гармонической функции), вычисления кратного интеграла, полученного из определения функции Грина, так и новый «векторный» алгоритм «блуждания по шарам», позволяющий находить одновременно несколько значений итераций оператора Грина. Наилучшие результаты показали алгоритмы «блуждания по шарам и сферам» и «блуждания по шарам», которые кратко описаны ниже.

Задача (2) эквивалентна интегральному уравнению и — А См + С/, ре-

00

шение которого при малых А представляется рядом и = А" С"+1 /•

п=0

Оценка степеней оператора Грина методом «блуждания по шарам и сферам». В трёхмерном случае степени оператора Грина можно вычислять на траекториях процесса «блуждания по шарам и сферам», описанном в [10], который естественно обобщается на случай комплексного параметра.

Пусть уравнение (2) однозначно разрешимо и0 ^ |Л| ^ с2, с = const, тогда интегральное представление для и имеет вид:

u(x) = 4dfe) / ^н5^1-«) /

Л'н(і)

Яс , ч c2sh[(/?-r)cl гдс q = у) = Mshific) - л с] - плот,,ость пирсхода из ж в у;

г = |ж - у\.

Тогда несмещённая оценка для п + 1-й степени оператора Грина имеет вид:

Л*.,,,),

с2

С(» + 1) — V^ V

« - -з^г0»«^'»«")-

В [10] доказано, что описанный процесс с вероятностью! сходится к границе области и N < оо с вероятностью 1.

Оценка степеней оператора Грина методом «блуждания по шарам». Для последовательности ип = G" / справедливо равенство Д а„ = -ií„-i, откуда следует, что

ип(х) = E(uh(xi) + s(x, xl)uri-i(xi)),

где cci точка, равномерно распределённая в шаре. Отсюда для Uk(xi) получены формулы:

N

щ(х0) = E^s(xi-!, х,)/(х,);

i=i

N N N

¿=i i=i N N N

U3{x0) = хг)и2{х,) = хг) ^ x^u^xj) =

i=i i=i j=¿+i N N N

i=i j=i+i ¿=j+i

которые очевидным образом позволяют оценивать конечное число степеней оператора Грина на одной траектории «блуждания по шарам».

Проведено сравнение полученных оценок с аналитическими значениями для модельной задачи. Рассмотрена возможность использования «аналитического продолжения решения» для расширения области применимости данных оценок при решении уравнения Гельмгольца.

В третьей главе представлены теоретические и практические результаты данной диссертации, касающиеся нахождения оценок первого собственного числа оператора Лапласа с помощью оценки степеней оператора Грина.

Рассмотрим задачу на собственные значения:

Аи + Хи = 0;1

I п (3)

и|г = О,J

где Л постоянное, возможно комплексное, значение. Известно [11], что множество собственных значений {А^} этой задачи не имеет конечных предельных точек, \k > 0, Ai простое (Ai < Аг < A3 < ...). Собственная функция ipi{x) неотрицательна, а щ можно выбрать вещественными и ортонормиро-ванными в С2{Т>).

Первый раздел содержит описание алгоритма вычисления первого собственного числа для оператора Лапласа с помошыо метода Келлога.

Второй раздал содержит описание алгоритмов оценки степеней оператора Грина для оператора Лапласа методами Монте-Карло, используемых для оценки первого собственного числа. Рассмотрена модельная задачи, и проведено сравнение полученных оценок с аналитическими значениями для модельной задачи.

Использование вероятностного представления для оценки итераций оператора Грина. С оператором Лапласа связан винеровский процесс, а именно является характеристическим оператором процесса. Пусты = т(х) момент первого выхода винеровского процесса £(i) из множества V, если 4(0) = ж, /(ж) = <ро = const > 0 и А < 0 и |А[ < Аь тогда для решения краевой задачи (2) справедливо вероятностное представление

0

Отсюда может быть получено представление

Е тп

и„. = G" ¥>0 = Vo^r

Несмещенной оценкой vn(x) = G"1 является случайная величина £„ = = N^/n\c2l\ где N ■ число переходов в шар на траектории «блуждания по

шарам и сферам», УУ'"' = Л^іУ - 1) ■ ■ • (М - п + 1). Оценка имеет конечную дисперсию.

Оценка итераций оператора Грипи на основе свойств броуновского даи-

Е т"

■жения. Рассмотрим представление С</>о = «А)^-]"' где г момент выхода броуновской траектории на границу области. Пустьжо = х,ху, ж2,..., ж,-,... процесс «блуждания по сферам», ть тг,..., г,,... времена выхода броуновской траектории на сферы, В, <т-алгебра, порождённая «блужданием по сферам» до момента времени і. п ЕтР

Для Ур = несмещенные оценки имеют вид

Ы*) = ¿ї Е {(ті + т2 + ... 4- Ті + ті)" І Вг} ,

где Ті - момент первого после Ті + Т2 + ■■■ + тг выхода броуновского движения на границу. На фиксированной траектории случайные величины 7~1, Т-2, . . . , Ті, тг независимы.

Пусть «о = 1, Уп(х) вектор ііц, Уі,..., ип, а V,, среднее значение V,, на сфере. Можно показать, что Уп(х) = АУп(х), где А нижнедиагональная матрица вида А = (I + -ш^Н + ш-2Я4Н2 + ... + ■шІІІїі"Н"), и:,а =

2"'т!

цг моменты случайной величины Ті, а матрица Н имеет отличные от нуля элементы = 1. Данное соотношение позволяет построить векторный

процесс, компоненты которого образуют мартингалы несмещённых оценок для ир(х), р = 1,2,..., п.

Пусть го, /'і, т'2, ■ • • -- радиусы сфер в процессе «блуждания по сферам» хо = х, х\, Х2,..., -- момент первого попадания блуждания в є-окрсстность границы, тогда

0\'г = Л(г0)Л(Г1) • • • А{гк.!)Уп{хк)

является несмещённой оценкой У{х) с конечной дисперсией, е-смещённая оценка У(х) с конечной дисперсией получается, если положить = 1,

а ук(хце) - 0.

С учётом структуры матрицы А, получено представление С, = А(го)^(п) х х •••!Д(г,,_1) в виде

С, = / + 7ІЯ + ... + 7п#", 1 = 1,2,....,

к

где ук определяются формулами = ги^Гц4'; 7= 7}и!к-іг^к~1\

1=0

В трёхмерном случае значение коэффициентов wt можно вычислить, используя рекуррентную формулу Wo = 1; Wn = —Wn i - -Ш„ 2 + ■■• + (-1)"+1

+ у—-—-w;n. Для пространства размерности более, чем 3, коэффициенты (2 п + 1)!

Wk можно определить, воспользовавшись результатами [12].

Четвёртая глава посвящена практической реализации описанных алгоритмов.

В первом, разделе описаны особенности реализации методов Монте-Карло при параллельных вычислениях на MPI-кластерах и графических процессорах общего назначения. Обращено внимание на особенности выполнения программного кода на графических процессорах и реализацию параллельного генератора псевдослучайных чисел.

Второй раздел содержит описание алгоритма, позволяющего моделировать плотность вероятности перехода р(х, у) = qön{y) + (1 - <])р(х, у), где

6п равномерное распределение на сфере с центром в точке аз, а р{х, у) = £.2 Ы R _

=-----, без использования метода отбора, что позволяет эффек-

4тгг sh cR - cR

тивно производить моделирование на графических процессорах.

Третий раздел содержит краткое описание программного комплекса, разработанного для реализации описанных алгоритмов.

Комплекс состоит из двух частей: универсальной оболочки для запуска расчётов и модулей, выполняющих расчеты (URL: http://teachers. uni-vologda.ac.ru/stat_mod/dprg.htm). Разделение на две части позволило уменьшить дублирование кода и упростить отладку и исправление ошибок. Реализованы механизмы чтения параметров из конфигурационного файла, что позволяет производить моделирование с различными значениями параметров без перекомпиляции программ и обеспечивает возможность пакетного выполнения расчётов. Благодаря выводу результатов в формате XML, возможно использование XSLT-шаблонов для представления полученных результатов в требуемых форматах, например, в виде HTML для просмотра результатов в браузере, или в формате ЬЯ]еХ для помещения полученных результатов в статьи и печатные формы. Разработанный код может быть скомпилирован как под ОС Windows, так и под ОС Linux. «Универсальная оболочка» обеспечивает:

• загрузку модулей расчёта;

• чтение параметров из XML-файла;

• вывод результатов в XML-файл;

• учёт времени расчёта (процессорное/астрономическое);

• вывод дополнительной информации (оценка времени работы, вспомогательная и отладочная информация от модуля расчёта, информация о платформе и аппаратном обеспечении).

В вычислительном модуле производится:

• проверка и обработка переданных параметров;

• реализация вычислительного цикла;

• вывод необходимой информации.

Для унификации и упрощения реализации вычислительных модулей со-здань! компоненты, реализующие

• генераторы случайных чисел;

• базовый класс для модуля расчёта, содержащий

- функции-обёртки для запуска на МР1-кластере;

— функции-обёртки для запуска расчёта на графических процессорах;

- код взаимодействия с «оболочкой»;

— код для вывода промежуточных результатов.

В заключении перечислены основные результаты диссертации, показана их теоретическая и практическая ценность.

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

1. Кублапоаская, В. II. Применение метода аналитического продолжения посредством замены переменных в численном анализе / В. Н. Кублановская // Работы по приближенному анализу. - Т. 53 из Тр. МИ АН СССР. - М.-Л.: 1959. - С. 145-185.

2. Тихонов, А. Н. Уравнения математической физики / А. П. Тихонов, А. А. Самарский. — М.: Наука, 1977.

3. Nabors, К. Fast capacitance extraction of general three-dimensional structures / K. Nabors, S. Kim, J. White // IEEE Transactions on Microwave Theory and Techniques. — 1992. — Vol. 40, no. 7.- Pp. 1496-1506.

4. Phillips, J. R. A precorrected-FFT method for electrostatic analysis of complicated 3-D structures / J. R. Phillips, J. K. White // IEEE Trans. On Computer-Aided Design Of Integrated Circuits And Systems. — 1997. — Vol. 16, no. 10. — Pp. 1059-1072.

5. Iverson, R. B. A stochastic algorithm for high speed capacitance extraction in integrated circuits / R. B. Iverson, Y. L. Le Coz // Solid-State Electronics. - 1992. - Vol. 35, no. 7. -Pp. 1005-1012.

6. Iverson, R. B. A floating random-walk algorithm for extracting electrical capacitance / R. B. Iverson, Y. L. Le Coz // Mathematics and Computers in Simulation. — 2001. — Vol. 55, no. 1-3. - Pp. 59-06.

7. Mascagni, M. The random walk on the boundary method for calculating capacitance / M. Mascagni, N. A. Simonov // The Journal of Computational. Physics. — 2004.— Vol. 195, no. 2. - Pp. 465-473.

8. Si-pin, A. S. Random walks in the domain and their applications to the boundary value problems / A. S. Sipin // International Conference «Tikhonov and Contemporary Mathematics», Section 3.- Moscow: 2006.-June 19-25.- Pp. 113-114.

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

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

11. Владимиров, В. С. Уравнения математической физики / В. С. Владимиров, — М.: Наука, 1967. - 436 с.

1

12. Сабельфельд, К. К. Методы Монте-Карло в краевых задачах / К. К. Сабельфельд; Под ред. Г. А. Михайлова. — Новосибирск: Наука. Сиб. отделение, 1989.

Публикации автора по теме диссертации Статьи в рецензируемых журналах и изданиях:

Al. Кузнецов А. Н., Сипин А. С. Универсальный алгоритм расчета взаимных электростатических емкостей системы проводников методом Монте-Карло. // Матсм. моделирование — 2009. Т. 21, № 3. — С. 41-52.

А2. Кузнецов А. Н., Сипин А. С. Статистические оценки для степеней оператора Грина. // Вести. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. -2009. - Т. 2(19) - С. 114 -123.

A3. Кузнецов А. Н., Рытеикова И. А., Сипин А. С. Оценки методом Монте-Карло итераций оператора Грина и собственных чисел первой краевой задачи для оператора Лапласа. // Всстн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. 2011. Т. 4(25) С. 82 92.

Другие публикации:

А4. Кузнецов А. Н. Расчет взаимных электростатических емкостей системы проводников методом блуждания по полусферам. // Труды пятой Всероссийской научной конференции с международным участием (29 31 мая 2008 г.). Часть 2. Моделирование и оптимизация динамических систем и систем с распределенными параметрами. - Матсм. моделирование и краев. задачи. - СамГТУ, Самара: 2008. - С. 58-60.

Подписано к печати 06.04.12. Формат 60x84 'А . Бумага офсетная. Гарнитура Тайме. Печать цифровая. Печ.л. 1,00. Тираж 100 экз. Заказ 5421.

Отпечатано в Отделе оперативной полиграфии химического факультета СПбГУ 1УИ504, Санкт-Петербург, Старый Петергоф, Университетский пр., 26 Тел.: (812) 428-4043, 428-6919

61 12-1/944

Вологодский государственный педагогический университет

на правах рукописи

Кузнецов Андрей Николаевич

Решение некоторых задач математической физики методами

Монте-Карло

Специальность 05.13.18 — Математическое моделирование, численные

методы и комплексы программ

ДИССЕРТАЦИЯ

на соискание учёной степени кандидата физико-математических наук

Научный руководитель: кандидат физико-математических наук доцент Сипин А. С.

Вологда 2012

Оглавление

Введение 5

1. Нахождение взаимных электростатических ёмкостей 12

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

1.1.1. Электростатическая ёмкость............................15

1.2. Краткое описание существующих алгоритмов вычисления ёмкостей............................................................17

1.2.1. Аналитические методы нахождения электростатических ёмкостей..............................................17

1.2.2. Алгоритмы «сопоставления с эталоном»................19

1.2.3. Алгортмы вычисления погонных ёмкостей ............20

1.2.4. Алгоритмы вычисления трёхмерных ёмкостей .... 31

1.3. Решение задачи в случае уединённого проводника............69

1.3.1. Нахождение ёмкости......................................69

1.3.2. Вычисление нормальной производной..................70

1.3.3. Алгоритм вычисления ёмкости..........................71

1.3.4. Моделирование переходной плотности..................72

1.3.5. Экспериментальные данные..............................73

1.4. Решение задачи в общем случае................................74

1.4.1. Нахождение взаимных ёмкостей........................74

1.4.2. Вычисление нормальной производной..................75

1.4.3. Оценка взаимных ёмкостей «блужданием по сферам» 76

1.4.4. Оценка взаимных ёмкостей «блужданием по полусферам» ........................................................78

1.4.5. Экспериментальные данные..............................81

2. Решение задачи Дирихле для уравнения Гельмгольца 92

2.1. Решение задачи без использования оценок степеней оператора Грина..................................................93

2.1.1. Алгоритм «блуждания по сферам»......................93

2.1.2. Алгоритм «блуждания по шарам»......................94

2.2. Решение задачи через оценки степеней оператора Грина . . 96 2.2.1. Аналитическое продолжение решения..................96

2.3. Вычисление оценок степеней оператора Грина методами Монте-Карло......................................................99

2.3.1. Метод «блуждания по шарам»..........................99

2.3.2. Метод «блуждания по шарам и сферам» .......100

2.3.3. «Блуждание по сферам» с выделением объёмного потенциала ..........................101

2.3.4. Нахождение степеней оператора Грина из определения функции Грина......................................102

2.3.5. Экспериментальные данные...............102

3. Оценка первого собственного числа для оператора Лапласа 109

3.1. Оценка первого собственного числа первой краевой задачи

для оператора Лапласа ..........................................110

3.1.1. Метод Келлога............................................110

3.1.2. Вероятностный смысл метода Келлога для оператора Грина......................................................111

3.2. Способы оценки итераций оператора Грина для функции

Мх) = 1 ..........................................................112

3.2.1. Использование оценки моментов случайной величины г 112

3.2.2. Оценка итераций на основе свойств броуновского движения ......................................................116

3.2.3. Модельная задача........................................120

3.3. Экспериментальные данные...................122

4. Практическая реализация алгоритмов 126

4.1. Особенности реализации параллельных вычислений.....126

4.1.1. Реализация генератора случайных чисел ..............127

4.1.2. Особенности реализации с использованием технологии MPI.................................128

4.1.3. Особенности реализации вычислений на GPU..........129

4.1.4. Алгоритм суммирования с компенсацией ..............131

4.2. Моделирование переходной плотности..........................132

4.3. Краткое описание программного комплекса....................134

4.3.1. Взаимодействие вычислительного модуля с оболочкой 135

4.3.2. Базовый класс вычислительного модуля................135

Заключение 138

Литература 140

Введение

Актуальность темы

С расширением использования сверхбольших интегральных схем и развитием высокочастотной электротехники значительное влияние на работу электронных систем стали оказывать значения паразитных элементов сопротивления, индуктивности и ёмкости. Для вычисления паразитных элементов и определения их влияния на части системы на различных этапах проектирования и изготовления разрабатываются системы автоматизированного проектирования микроэлектромеханических систем (Microelectromechanical Computer-Aided Design system, MEMCAD) [1]. Основными элементами таких систем являются: программа для моделирования структуры объекта (structure Simulator), база данных со свойствами материалов, модуль механического анализа и модуль электростатического анализа (рис. 0.0.1, [1]). В зависимости от методов, используемых при электростатическом и механическом анализе, трёхмерная модель может Рис о.0.1. Архитектура являться приближением исходного объекта с помощью тех или иных примитивов

геометрических, объектов со схожими свойствами и пр.), так и выражать-

База данных материалов

Электростатический анализ Механический анализ

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

Так как при проектировании схем может потребоваться моделирование сотен или даже тысяч различных их конфигураций, к программам, вычисляющим значения паразитных параметров, выдвигаются высокие требования по быстродействию и точности вычислений. Аналитическое решение данной проблемы возможно лишь в случае небольшого количества проводников достаточно простой формы (параллелепипеды, сферы, цилиндры), поэтому для оценки взаимных электростатических ёмкостей трёхмерных систем из большого количества объектов, требуется применение тех или иных способов численного решения этой задачи [2-13].

В работе рассмотрена возможность применения методов Монте-Карло для нахождения взаимных электростатических ёмкостей проводников. Показано, что методы «блуждания по сферам» и «блуждания по сферам и полусферам» дают хорошую оценку искомых величин и могут быть использованы для решения поставленной задачи.

Знание первого собственного числа для оператора Лапласа может потребоваться для оценки применимости различных алгоритмов решения краевых задач, например, при решении задачи Гельмгольца с помощью «блуждания по шарам и сферам». В работе рассматриваются алгоритмы нахождения оценок итераций оператора Грина для оператора Лапласа и первого собственного числа для оператора Лапласа. Исследуется возможность использования оценок итераций оператора Грина для расширения области применимости методов Монте-Карло при решении уравнения Гельмгольца.

Использование метода Монте-Карло для решения описанных задач позволяет заметно упростить реализацию алгоритма для параллельного вычисления, в том числе на графических видеопроцессорах (General Purpose Graphical Processing Unit, GPGPU), и не требует значительных затрат для

хранения промежуточных данных.

Цели работы

• Разработка алгоритмов вычисления взаимных электростатических ёмкостей методами Монте-Карло.

• Исследование оценок итераций оператора Грина с целью вычисления первого собственного числа оператора Лапласа и решения уравнения Гельмгольца.

• Реализация перечисленных алгоритмов для вычислений с использованием MPI-кластеров и графических процессоров общего назначения (видеокарт).

• Создание программного комплекса, упрощающего «обычную» и «параллельную» реализацию указанных программ.

Методика исследования

Методика исследования включает применение методов статистического моделирования к решению краевых задач математической физики. С помощью формул Грина краевые задачи сводятся к интегральному уравнению, которое затем решается методом Монте-Карло. При решении уравнения Гельмгольца используется аналитическое продолжение решения методом замены переменных [14]. Для нахождения первого собственного числа краевых задач используются известные итерационные методы [15]. Реализация алгоритмов выполнена на языке программирования С++ с использованием реализации MPI MPICH-2 [16] и технологии CUDA [17]

Практическая значимость

Разработанные статистические алгоритмы вычисления взаимных электростатических ёмкостей проводников могут быть использованы для оценки паразитных электростатических ёмкостей в больших и сверхбольших интегральных схемах и многослойных печатных платах. Проведены вычислительные эксперименты для алгоритмов, использующих методы Монте-Карло при решении задачи Дирихле для уравнения Гельмгольца. Рассмотрены различные методы вычисления итераций оператора Грина методами Монте-Карло. Проведены вычислительные эксперименты, позволяющие оценить возможность использования аналитического продолжения решения методом замены переменных. Представлены новые алгоритмы вычисления оценок первого собственного числа для оператора Лапласа с помощью оценок степеней оператора Грина. Разработан программный комплекс для реализации вычислительных задач, использующих методы Монте-Карло. Программный комплекс позволяет упростить реализацию и проведение вычислительных экспериментов как на выделенных компьютерах, так и на МР1-кластере или графическом процессоре общего назначения.

Апробация работы

Основные результаты обсуждались на семинарах кафедры прикладной математики в Вологодском государственном педагогическом университете, семинаре кафедры статистического моделирования математико-механического факультета СПбГУ и докладывались на

• пятой Всероссийской научной конференции с международным участием «Математическое моделирование и краевые задачи» (СамГТУ, Самара, 29-31 мая 2008 г.);

• II ежегодном смотре-сессии аспирантов и молодых учёных по отраслям наук. Секция «Математика. Информатика. Методика преподавания» (Вологда, 12-14 ноября 2008 г.); в научной школе-конференции молодых исследователей «Математические идеи П. Л. Чебышева и их приложение к современным проблемам естествознания» (Обнинск, 14-18 мая 2011 г.).

Основные результаты опубликованы в работах [18-21].

Структура и объём диссертации

Диссертация состоит из введения, 4-х глав, заключения и списка литературы.

В главе 1 рассматривается проблема нахождения взаимных электростатических ёмкостей.

Раздел 1.1 содержит постановку задачи нахождения электростатической ёмкости уединённого проводника и взаимных электростатических ёмкостей нескольких проводников.

В разделе 1.2 приведено краткое описание существующих алгоритмов нахождения взаимных электростатических ёмкостей.

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

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

аналитическими решениями и результатами программ РаяЮАР и РРТСАР.

Глава 2 посвящена алгоритмам вычисления оценок степеней оператора Грина для уравнения Лапласа и решению задачи Дирихле для уравнения Гельмгольца. В ней представлены теоретические и практические результаты данной диссертации, касающиеся нахождения оценок степеней оператора Грина и применения метода аналитического продолжения решения к задаче Дирихле для уравнения Гельмгольца.

В разделе 2.1 рассмотрены алгоритмы решения задачи Дирихле для уравнения Гельмгольца методами «блуждания по шарам и сферам» [22] и «блуждания по шарам» [23]. Эти алгоритмы не используют оценок степеней оператора Грина. Они были реализованы для сравнения с разработанными нами методами.

Раздел 2.2 содержит представление решения задачи Дирихле для уравнения Гельмгольца в виде суммы степеней оператора Грина. Приведено описание метода аналитического продолжения решения [14] посредством замены переменных.

Раздел 2.3 посвящён вычислению оценок степеней оператора Грина для оператора Лапласа методами Монте-Карло. Рассмотрены как стандартные методы построения с помощью «блуждания по шарам и сферам», «блуждания по сферам» (через представление решения в виде суммы объёмного потенциала и гармонической функции), вычисления кратного интеграла, полученного из определения функции Грина, так и новый «векторный» алгоритм «блуждания по шарам», позволяющий находить одновременно несколько значений итераций оператора Грина.

В разделе 2.3.5 проведено сравнение полученных оценок с аналитическими значениями для модельной задачи. Рассмотрена возможность использования «аналитического продолжения решения» для расширения области применимости данных оценок при решении уравнения Гельмгольца.

В главе 3 представлены теоретические и практические результаты дан-

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

Раздел 3.1 содержит описание алгоритма вычисления первого собственного числа для оператора Лапласа с помошью метода Келлога.

Раздел 3.2 содержит описание алгоритмов оценки степеней оператора Грина для оператора Лапласа методами Монте-Карло, используемых для оценки первого собственного числа.

В разделе 3.3 рассмотрена модельная задачи и проведено сравнение полученных оценок с аналитическими значениями для модельной задачи.

Глава 4 посвящена практической реализации описанных алгоритмов.

В разделе 4.1 описаны особенности реализации методов Монте-Карло при параллельных вычислениях на МР1-кластерах и графических процессорах общего назначения. Обращено внимание на особенности выполнения программного кода на графических процессорах и реализацию параллельного генератора псевдослучайных чисел.

Раздел 4.2 содержит описание алгоритма, позволяющего моделировать плотность вероятности переходар(х, у) = + — я)р(х, у), где 5ц —

равномерное распределение на сфере с центром в точке х, а р(ж, у) =

с2 вЬ(сД — сг) ^ ^

=--—, без использования метода отбора, что позволяет эф-

47ГГ эЬ сК — сН

фективно производить моделирование на графических процессорах.

В разделе 4.3 содержится краткое описание программного комплекса, разработанного для реализации описанных алгоритмов.

В заключении перечислены основные результаты диссертации, показана их теоретическая и практическая ценность.

Глава 1.

Нахождение взаимных

электростатических

ёмкостей

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

Если не оговорено отдельно, формулы и значения электростатических ёмкостей будут приводится согласно системе СГСЭ («сантиметр-грамм-секунда» электростатическая). В этом случае £о = 1 и несколько упрощаются математические выкладки и вычисления. Также считается, что проводники разделены однородным изотропным диэлектриком с диэлектрической проницаемостью 1.

Согласно [24] основной задачей электростатики является отыскание поля, создаваемого системой зарядов на заданных проводниках.

Прямая постановка предполагает нахождение поля вне проводников и плотностей зарядов на проводниках при известных потенциалах проводников. Решение задачи сводится к нахождению функции </?, удовлетворяющей уравнению Лапласа Д</? = 0 всюду вне заданной системы проводников, об-

ращающейся в нуль на бесконечности и принимающей заданные значения срг на поверхностях проводников Г^:

То есть мы имеем первую краевую задачу для уравнения Лапласа, которая имеет единственное решение. Процесс решения внешней задачи Дирихле для уравнения Лапласа методом Монте-Карло подробно описан в [23].

Приведём краткое описание алгоритма «блуждания по сферам», описанного там. Пусть V — ограниченная область в Мт, т > 3, Г — её граница, Х>1 = Жт \ Т), ср — непрерывная на Г функция, а Кц — шар радиуса К с центром в нуле, такой, что V С К л, тогда для точки х ф справедлива формула, дающая решение внешней задачи Дирихле для шара

где ат — площадь поверхности сферы радиуса 1 в Мт;

— сфера радиуса Я с центром в нуле. Тогда марковская цепь в Т>\ может быть определена следующим образом: Хо :=х; если после п шагов хп ^ Кц, то вектор хп+\ распределён на сфере с плотностью

иначе хп+\ равномерно распределён на сфере Зг(хп) радиуса г с центром в хп, где г — сИ81(жп,Г). Построенная марковская цепь с вероятностью 1 сходится к границе области V, а дисперсия полученной оценки ограничена.

В случае обратной постановки считаются известными полные заряды проводников Цг. Искомыми величинами являются потенциалы проводников, распределение зарядов по их поверхностям и поле вне проводников.

А<р = 0;

(¿>1 = щ = ал^ .

1 г

(1.1.1)

То есть, требуется найти функцию ср, удовлетворяющую уравнению Лапласа А(р = 0 вне заданной системы проводников, обращающуюся в нуль на бесконечности, принимающую некоторы�