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

С'. М. Ермаков, А. И. Рукавишникова, К. А. Тимофеев

О НЕКОТОРЫХ СТОХАСТИЧЕСКИХ И КВАЗИСТОХАСТИЧЕСКИХ МЕТОДАХ РЕШЕНИЯ УРАВНЕНИЙ*

В работе рассматриваются методы Монте-Карло (МК) и квази Монте-Карло (КМК) для решения некоторых типов уравнений в случае, когда решение представлено в виде обобщенного (в смысле главного значения) интеграла по траекториям. Для линейных уравнений указаны методы оптимизации алгоритма. Приводится пример рандомизации алгоритма вычисления решения сеточного аналога уравнений Навье-Стокса. Обсуждаются проблемы распараллеливания алгоритмов.

1. Введение

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

2. Линейный случай

Рассматривается интегральное уравнение второго рода

где x £ D С Rs, f и k —заданные функции на носителях меры « и « <g> « соответственно, а равенство по (mod «) обозначает его выполнение на носителе «. Если оператор с ядром k действует в нормированном пространстве F, F £ L(«), где Ь(р) —пространство («-интегрируемых функций, и норма оператора с ядром k(x, у) меньше единицы,

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

© С. М. Ермаков, А. И. Рукавишникова, К. А. Тимофеев, 2008

(mod «),

(1)

то, как известно [1], для любой функции Н € Г * возможно представление величины (Н, <^>) = / Н(х)^(х)р(йх) в виде интеграла по траекториям некоторой марковской цепи от специально построенной функции Ф(ш) траекторий ш. При этом предполагается интегрируемость Ф, что, в свою очередь, предполагает сходимость по норме Г мажорантного итерационного процесса

Трт+1 = J\к{х,у)\1рт{х)ц{<],х) + \/(х)\ тос1^, то = 0,1, Тр0 = /.

Интегралы такого рода могут вычисляться, в частности, с помощью МК метода, что предполагает моделирование траектории соответствующей марковской цепи и вычисление функции Ф на этой траектории. Погрешность при этом убывает как 0(Ы-1/2), где N — число промоделированных траекторий.

Применение детерминированного метода интегрирования КМК, сходного по своей структуре с методом МК, сталкивается со следующей трудностью. При вычислении интеграла по ¿-мерному гиперкубу от функции ограниченной вариации погрешность имеет порядок 0((1плN)/N). При каждом фиксированном й порядок убывания погрешности здесь асимптотически лучше, чем для метода МК, но интеграл по траекториям имеет очень большую кратность (теоретически бесконечную) и применение методов КМК становится проблематичным.

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

Пути преодоления этих трудностей указаны в ряде публикаций авторов и их коллег [2-4]. В данной работе получены новые результаты по обоснованию и улучшению предложенных ранее алгоритмов. В нелинейном случае построены и отработаны новые алгоритмы.

В работах [2-4] в качестве альтернативы известным методам, основанным на моделировании цепей Маркова, предлагалось последовательно оценивать приближения ^п(х) = / к(х, у)(рп-1(у)^(с1у) + /(х). Речь идет о вычислении функции <£>„ во всех точках носителя меры ц при вычисленной заранее у>п_1. Как обычно, при решении задач методом МК рассматриваются «прямые» и «сопряженные» [1] оценки. Прямые предполагают оценивание <^п(х) в заранее выбранных точках. Пусть х' — точка, в которой мы оцениваем у>п. Одной из возможных оценок <^п(х) может быть

,()=ад^*) (Д

Рп-1(У)

где у есть случайная величина, распределенная с плотностьюрп-1(у), а фп-1(у) —оценка у>п-1, полученная на предыдущем шаге. Эта плотность должна удовлетворять условию согласования

Рп-1(у) > 0 для у, удовлетворяющих условию к(х', у)^п-1(у) = 0. (3)

Оценка ^1 (у) очевидным образом является несмещенной оценкой ^>п(х') при условии (3). Получив N независимых реализаций случайной величины у^ вектора у, мы можем

оценивать ^п(х') с заданной точностью с помощью величины

1 *

шю = й'Е$ш- (4)

¿=1

При таком способе вычислений предполагается, что <^п+1 будет оцениваться после вычисления (оценки) фп в некотором множестве точек х[, ...,х'я и указания методов интерполирования (восполнения)для вычисления фп+1 в любой точке носителя р.

Другую возможность предоставляют оценки фп в случайно выбранных точках. Пусть рп(х,у) есть совместное распределение х и у, так что //рп(х,у)р(йх)р(йу) = 1 и выполняется условие согласования

Рп(х,у), если к(х,у)фп-1(у) =0. (5)

Получив независимые реализации (х, у^) случайных векторов, распределенных с плотностью рп(х,у), мы можем получить случайный вектор

Ь(х, у) = /(*) + к{Х^Пу^Уг\ * = 1,2,... (6)

Для любой функции Н € Г * очевидным образом выполняется неравенство

Е(^2(хг,уг)Н(хг)) = (Н,/) ^ Л Н(х)к(х,у)фп-1 (у)^((Ъ)р((1у) = (Н,фп). (7)

Это, также, означает, что случайный вектор

N

1/^12 ЫхиУг) І=1

N

слабо сходится по

І= 1

вероятности к фп(х) в каждой точке X.

Теорема о существенной выборке ([1] стр. 138) позволяет решить вопрос об оптимальном в смысле величины дисперсии выборе рп при фиксированном фп-1 в каждом из двух рассмотренных случаев. Имеет место следующий результат.

Теорема 1. При фиксированном фп-1 для каждого п = 1, 2, . .. справедливы неравенства

Ю£і - (У 1к(х 1 ’У^п-іШКЛу^ - (/ к(Х,У)‘Рп-1(уЖЛу)^ ,

(8)

Ік(х,у)фп-1(у)Ік(<іх)іл(<іу)^ - (^11 к(х,у)фп-1(у)к(Лх)к(Лу)^ . (9)

Причем, знак равенства достигается в неравенстве (8) при

„ _ \Нх',у)Фп-і(у)\ ПГі,

” І\Кх',у)Фп-і{у)\^у)

и в неравенстве (9) при

р(Ху) = ________\Нх,у)Фп-і(у)\_______

" ’ И\к(х,у)фп-і(у)\^х)^у)

При этом рп(у) есть плотность по отношению к мере и, а рп(х,у) —по отношению к мере и ® и.

Доказательство. Имеем

вс. = 4Ыз-' «№п-М + /(Л = „ /к(*‘,у)Фп-М

) \ Рп-Ау)

\ Ри-Лу)

к2(х ,,у)$П-1(у)

-у р ^ ]~р{<1у) ~ {^1 к(х,,у)фп^1(у)ц(йу)^ >

>^! \к(х',у)фи-і(у)\КЛу)^ - (/ к(х' ,У)!?и-і{у)р{Лу)^

так как

к2(х',у)ф2п_1(у)

Рп-і(у)

и , [ ^(х^у^І^у)

= ] (рп_!(у))2 Р”-Чу)№у) -

>

\к(х ,у)срп-1(у) Л

-----------------Рп-і(уЖ<Іу) =

Ри-і{у) )

\к(х ',у)Фи-і(у)\р(<'1у)) ,

(неравенство Буняковского). Отсюда следует (8). Путем вполне аналогичных выкладок можно получить (9). Прямые выкладки также показывают, что знак равенства в (8) и (9) достигается для плотностей рп, задаваемых выражениями (10) и (11) соответственно.

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

Таким образом, имеем

фп+1 = К фи + / + вп+і, п = 0, 1,

(12)

Здесь К есть реализация случайного оператора, оценивающего методом Монте-Карло функцию Кфп при фиксированном п, а вп+і —систематическая погрешность, возникающая в результате восполнения на п +1 шаге. Предполагаем, что используются несмещенные оценки, так что Е(Кфп)\фп) = Кфп.

Если Ефп = фп + Фп, то в силу равенства фп+і = Кфп + / имеем «п+і = К«п + вп+і. Вводя обозначения V = тах ||«п||, в = тах ||вп||, получим неравенство V < ||КУ« + в,

п

V <

1 -ЦК II'

(13)

По предположению ||К|| < 1, а из теории непараметрического оценивания [5] следует, что при соответствующем методе восполнения в стремиться к нулю с увеличением числа независимых испытаний Мп при каждом п. Порядок убывания в зависит от свойств функции ф и меры и Как мы видим, в случае дискретной меры и оказывается в = 0. Далее, мы предполагаем ф достаточно гладкой, чтобы обеспечить порядок убывания, по крайней мере, в = 0(Кп 1/2).

в

Вернемся к изучению дисперсии случайной ошибки с ростом п. Введем обозначения:

Єп = фп — фп, Єп = Єп — Еєп = Єп — «ги д = К — К- (14)

Из (1) и (12) имеем

Єп+і = Кєп + «п+і + дєп + §фп- (15)

Обозначим через рп+і(х, у) = р(єп+і) ковариационную функцию ошибки Є"п+і. По определению

Рп+і(х,у) = Е(Єп+і(х),Єп+і(у)).

Используя равенство (15) и неоднократно упоминаемую независимость К и єп, получаем равенство

Рп+і = КрпКт + Е(5рп5т) + «п+і«п+і + Е(5фпф*п5т) =

= КрпКт + Е(дрпдт)+ 0(М-і). (16)

Здесь Т в виде верхнего индекса у оператора обозначает оператор, сопряженный к

исходному в смысле Лагранжа.

Можно показать (как это сделано, например, в [6]) что норма линейного оператора, действующего на рп в правой части равенства (16), не превосходит величины ||К||2 + 0(М-і). В каждом конкретном случае возможен более детальный анализ этого оператора. Отметим только, что слагаемое Е(5рп5т) зависит не только от числа повторений Ып на п-м шаге, но и от дисперсии оценок, вычисляющих фп+іпри фиксированном фп. Это могут быть, например, оценки вида (4), (6). Вопрос об уменьшении их дисперсий обсуждался нами выше.

Таким образом, при соответствующем выборе чисел Ып могут быть использованы построенные нами оценки для (Н,фп). Применение в этой схеме методов КМК позволяет иметь алгоритм с конструктивной размерностью ¿. Определенные неудобства такого подхода состоят в том, что в схеме метода КМК мы не можем, вообще говоря, построить доверительный интервал, оценивающий погрешность в процессе решения задачи.

3. Эволюционные уравнения с квадратичной нелинейностью

Рассмотрим уравнение

| = (17)

где ф = ф(г,Х),Х = (хі, ...,хв) Є О С Кя —искомая вектор-функция, Г — заданный оператор. Если Г — дифференциальный оператор, то, используя разностную аппроксимацию на сетке (_О предполагается ограниченной), приходим к системе обыкновенных дифференциальных уравнений.

(18)

где и (г) = (и і (г), ...,щ(г)) —вектор неизвестных функций и2 (г) = ф(Ь,Х^ ), X^ —точки сетки. Полученную систему можно, при соответствующих предположениях относительно Г, решать одним из известных методов (Эйлера, Рунге—Кутты). В простейшем случае имеем схему

и (г + Аг) = и (г) + Агф(и (г)). (19)

В задачах с квадратичной нелинейностью

т т т

Ни ) = \\іі + Т, аукик + ЕЕ Ь^кі икґщЦ=і

к=1

к=11=1

и при большом в вычисление вектора (20) может потребовать при каждом £ очень большого числа арифметических операций.

Возможен подход, при котором вместо Р(и) вычисляется его несмещенная (или «слабо» смещенная) оценка. При этом использование оценок, аналогичных £2, предполагает вычисление случайно выбранного подмножества компонент и. Если это подмножество не велико (существенно меньше в элементов), то возможно получить нужные оценки решения с меньшим числом арифметических операций.

Заметим, что условие несмещенности делает возможным простое распараллеливание алгоритма — аналогичные вычисления ведутся на разных компьютерах с независимыми датчиками случайных чисел, результаты осредняются. Наличие квадратичного члена требует, однако, модификации оценок предыдущего раздела. Здесь можно воспользоваться приемом, описанным в [4], и вычислить поправку, получаемую путем разложения квадратичной функции в ряд Тейлора в окрестности математического ожидания оцениваемых величин.

Приведенные общие соображения были проверены на модельном примере.

Рассматривается разностный аналог уравнений Навье—Стокса с введенной искусственной вязкостью, описывающий движение вязкой несжимаемой жидкости:

Здесь V = у(Ь,х) = (у-\_(1,х),Ю2(1,х),Юз(1,х)) —векторное поле скоростей в трехмерном пространстве (х = (х1,х2,хэ) € К3). Величины р и V являются скалярными параметрами, определяющими плотность и вязкость жидкости соответственно. Коэффициент £ — параметр искусственной вязкости.

Приведенное уравнение использовалось для расчета движения жидкости в кубической каверне с подвижной боковой стенкой: х € П, где П = [0,1]3.

В начальный момент времени (£ = 0) жидкость полагалась неподвижной:

Пространственная сетка содержала I х I х I точек — одинаковое число по каждой переменной. Использовалась стандартная аппроксимация производных с точностью до Н2, Н = 1/1. Полагалось £ = 0.5, р = 1, V в разных примерах меняется от 0.0001 до 0.25.

Вычислялись оценки решения в момент времени £ = 1 и полагалось АЬ = 1/Ы. Обычные детерминированные вычисления по явной схеме (19) налагают соответствующие условия на соотношения Н и А1 вида А1 < сН2, с — некоторая константа, обеспечивающая устойчивость относительно ошибок округления. Во всех случаях эти условия были выполнены. Кроме того, рандомизация требует выполнения условий стохастической устойчивости. Они уменьшают величину константы с.

(21)

«(0, х) = 0.

(22)

На границах задавалось условие прилипания, а именно:

(23)

В случае использования детерминистического метода Эйлера необходимо будет провести 3 • (I — 2)3М расчетов отдельных компонент функции /. Если использовать при расчетах рандомизацию, то для фиксированного числа усреднений т оценки на каждом шаге получим, что для расчета оценки решения к моменту времени 1 потребуется тМ расчетов компонент _Р. При этом мы предполагаем, что величина т должна быть достаточной для обеспечения стохастической устойчивости и, следовательно, распараллеливания алгоритма. Как показывают численные исследования, в данном примере наименьшее значение параметра т, при котором сохраняется стохастическая устойчивость, убывает с уменьшением ДЬ (вплоть до 1!). Это сближает описанную методику с методами, использующими аппарат стохастических дифференциальных уравнений. Таким образом, оказывается, что существует Ьо такая, что для всех 0 < ДЬ < ¿о рандомизированный метод Эйлера быстрее детерминистического в Ь0/ДЬ раз.

Для ряда параметров I и ДЬ приведем наименьшие значения т, при которых сохраняется стохастическая устойчивость, а также покажем, во сколько раз уменьшается время расчета компонент вектора ^ при использовании рандомизированного метода Эйлера, по сравнению детерминистическим (см. табл.).

1 Д* наименьшее т 3(1 — 2)6 /т - преимущество ранд, метода

10 0.001 461 3.33

10 0.0005 183 8.39

10 0.0001 19 80.8

10 0.00005 1 1536

20 0.0001 1731 10.12

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

Интерес представляет также поведение погрешности оценки, получаемой в момент времени Ь путем усреднения нескольких независимых оценок решения, полученных на разных компьютерах. Эксперименты показывают, что для достаточно большого диапазона числа усредняемых оценок N (при расчетах N бралось от 1 до 512) норма отклонения усредненной оценки решения от оценки, получаемой детерминистическим методом Эйлера, убывает как с/а/Ж, где с — некоторая константа.

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

Ниже, на рис. 1 в качестве иллюстрации приведен один из результатов расчета стохастическими методами при Ь = 1, I = 30, ДЬ = 0.0001, т = 25 000, V = 0.1 и V = 0.0001 соответственно.

На следующих графиках приведены значения среднеквадратических отклонений компонент оценок VI и V2 решения в плоскости Оху при г = 0.5, домноженные на 106, а также эти же значения среднеквадратических отклонений, деленные на модули математических ожиданий соответствующих компонент (относительная погрешность), рассчитанные для Ь = 0.12, I = 10, ДЬ = 0.1225, т = 1200, V = 0.25.

----------- О 1

0x1 X

Рис. 1. Линии тока в плоскости Оху при г = 0.5 для V = 0.1 и V = 0.001.

Рис. 2. Среднеквадратические отклонения оценок компонент V1 и г>2 в плоскости Оху при 2 = 0.5, домноженные на 106.

Рис. 3. Относительные погрешности оценок компонент VI и V2 в плоскости Оху при 2 = 0.5.

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

2. Ermakov S. M., Wagner W. Monte Carlo difference schemes for the wave equation // Monte Carlo Methods and Appl. 2002. Vol. 8, N1. P. 1-29.

3. Ermakov S. M., Rukavishnikova A.I. Quasi Monte-Carlo Algorithms for Solving Linear Algebraic Equation // MC Methods and Applications. 2006. Vol. 12, N 5. P. 363-384.

4. Ermakov S. M. MCQMC algorithms for solving some classes of equations // Monte-Carlo & Quasi Monte-Carlo Methods / Springer 2007. P. 23-33.

5. Devroye L. Nonparametric density estimation (five lectures) // Principles of Nonparametric Learning. Wien: Springer-Verlag, 2002. P. 211-270.

6. Адамов А. В., Ермаков С. М. Стохастическая устойчивость метода Монте-Карло (случай операторов) // Вестник СПбГУ. 2004. Сер. 1. Вып. 2. С. 3-8.

7. Рауч П. Вычислительная гидродинамика. М.: Мир, 1980. С. 612.

Статья поступила в редакцию 18 мая 2008 г.