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

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

А. С. Сипин

Вологодский государственный педагогический университет, канд. физ.-мат. наук, доцент, [email protected]

1. Введение. В работе предложены алгоритмы статистического моделирования для решения задачи Коши для параболического уравнения второго порядка с переменными коэффициентами, в которых несмещенные оценки решения строятся на траекториях случайных блужданий.

Пусть в области ^„+1 = Еп х (0, Т) задан невырождающийся параболический оператор

( д д \ д п д2 п д

1-Чм"£-гм_з<-£ +2><*.‘)35;+«.<*,«, (1Л)

4 7 ъ,3 = 1 ъ=1

коэффициенты которого принадлежат классу На,а/2(пПП+>1), а < 1. Матрица

коэффициентов при старших производных предполагается симметричной, а ее собственные числа лежат в фиксированном отрезке [V, р] и V > 0.

Рассмотрим задачу Коши

( д д \

Ь ( М, —, — I и = /, и\г=0 = ^{х). (1.2)

В [1] показано, что уравнение (1.2) имеет фундаментальное решение Z(х,у,Ь,т). Пусть функция ] удовлетворяет условию Гёльдера по всем своим аргументам, р непрерывна, а / и р растут при |х| ^ то не быстрее е“|х| . Тогда решение задачи (1.2) может быть записано в виде суммы потенциалов:

и(х,Ь) = [ йт [ Z(х,у,Ь,т)f (у,т)йу + [ Z(х,у,Ь, 0)р(у)йу. (1.3)

•-! 0 о Еп •/ Еп

Константа а зависит от Т и коэффициентов уравнения.

Обсудим кратко известные методы решения задачи Коши. Для уравнений с постоянными коэффициентами задача (1.2) сводится к аналогичной задаче для оператора Ь = д/дЬ — Д, фундаментальное решение которого известно. Решение и(х,Ь) находится по формуле (1.3), в которой интегралы удобно вычислять методом Монте-Карло, так как фундаментальное решение является плотностью распределения вероятностей.

Если коэффициенты уравнения (1.1) и его правая часть не зависят от времени, то существует однородный диффузионный процесс Х8, стартующий из точки х, на траекториях которого справедливо вероятностное представление функции и(х,Ь):

и(х, Ь) = Ех J ехр ^ — J а0(Хт)йт^ f (Х8)йв + Ех ехр ^ — J ад(Хв)й^ р(Хг), (1.4)

* Работа выполнена при финансовой поддержке РФФИ (грант №11-01-00769). © А.С.Сипин, 2011

в котором символ Ех означает математическое ожидание следующей за ним случайной величины. Представление (1.4) явно указывает вид несмещенных статистических оценок для и(х,Ь) на траекториях процесса Х8. Отметим, что реализация этих оценок связана с моделированием траекторий случайного процесса. Для этого требуется решать систему стохастических дифференциальных уравнений численно, что приводит к смещенным оценкам для и(х,Ь) и затрудняет вычисление погрешности.

Наконец, для случая дифференцируемых коэффициентов уравнения можно получить, используя формулы Грина, интегральное уравнение для и(х,Ь) и решить его методом Монте-Карло. Такие алгоритмы рассмотрены в работах [2] и [3] для уравнений, главной частью которых является оператор Лапласа. На уравнение с переменной матрицей старших коэффициентов они непосредственно не переносятся.

В данной работе на траекториях легко моделируемых случайных блужданий построены несмещенные оценки решения задачи Коши и функционалов от него. В основе построений лежит формула (1.3), в которой каждый интеграл оценивается по одному случайному узлу. Фундаментальное решение является функционалом от решения некоторого интегрального уравнения Вольтерра, поэтому оценка для него находится по схеме Неймана—Улама [4].

2. Фундаментальное решение уравнения теплопроводности. В этом разделе представлены некоторые результаты из [1], связанные с построением фундаментального решения 2(х,у,Ь,т) для уравнения теплопроводности.

Фиксируем точку (у,т). Пусть А(у,т) —матрица, составленная из старших коэффициентов ац (у,т) оператора Ь, А(гц) (у,т) — элементы обратной матрицы А-1(у,т). Рассмотрим функцию 2о, которая при Ь > т определяется равенством

2о(х — у,у,і,т)

1

[47г(£ — т)]"/2 (сіеі А(у, т))1/2

(2.1)

При Ь < т функция 2о(х — у, у, Ь, т) = 0. Функция 2о удовлетворяет неравенству

(2.2)

где

Фундаментальное решение 2 можно представить в виде

Е,

2о(х — г, г,і, X)Q(z,y, Х,т)3,г, (2.3)

где функция Q является решением уравнения Вольтерра

Е

К(х, г, і, X)Q(z, у, X, т)<!г + К(х, у, і, т) = 0, (2.4)

а функция К определена равенством

, , \—л г , , , ,,д22о(х — z, z,t, X)

К(х, г, t, Л) = 22 1аг,з(г^) - аг,з(х^)\--------дх дх---------Ь

1,3=1 г 3

Е, ч д2о (х — z, z,t, X) . , , ,

аг(х, £)-----------—-------------\- ао(х, £).^о(ж — Л). (2.5)

дхг

1=1

Существуют положительные постоянные С и с, такие что при 0 ^ т < Ь

|К(х, у, т) | ^ с(£ — г)_(-"+2_“-)/2 ехр ^ . (2-6)

Пусть К1(х,у,Ь,т) = К(х,у,Ь,т). Для повторных ядер

Кт(х,у,Ь,т) = [ <1\ [ К(х^,Ь,Х)Кт-1^,у,Х,т)dz,m = 2, 3,... (2.7)

•у т •/ Еп

при 0 ^ т < Ь по индукции доказывается оценка

\Кт(х,уМ1 < с» (^)”(т-1)/2 ^||(^ -.)-(—)/2 ехр (-С^) , (2.8)

где Г(а) обозначает гамма-функцию.

Из оценки (2.8) следует, что ряд Неймана для уравнения (2.4) сходится равномерно при Ь — т > 0,

Q(x,y,t,т) = ^2( — 1)тКт(х,у,г,т), (2.9)

т=1

и справедливо неравенство

\С}{х,у,г,т)\ ^ С1(г - т)-^2-^/2 ехр ) • (2.10)

3. Несмещенные оценки решения задачи Коши. Из представления (1.3) для решения задачи Коши и формулы (2.3) для фундаментального решения следует, что решение задачи Коши распадается в сумму четырех потенциалов:

и(х, Ь) = и>1(х,Ь) + и2(х, Ь) + из(х, Ь) + и^(х, Ь) = г*

/ dт 2о(х — у, у, Ь, т)/(у, т^у +/ 2о(х — у,у,Ь, 0)ф(у№+

•-! о о Еп •/ Еп

/ dт / / dX / 2о(х — z,z,t,X)Q(z,y,X,т)dz /(у,т)dy+

Jо JEn\Jт -1еп )

+ [ ( [ dX ( 2о(х — z,z,t,X)Q(z,y,X, 0)dЛ ^(y)dy. (3.1)

Зеп \-)о Зеп )

В силу неравенства (2.2) для несмещенного оценивания ui(x,t) и u2(x,t) можно использовать плотность Z\. Нормально распределенный случайный вектор Y, имеющий плотность Zi(x — y,t — т), можно моделировать по формуле

У (ж, t, т) = х + \J 4 fj,(t — (3-2)

где y — случайная величина, имеющая гамма-распределение с параметром n/2, а ш — изотропный случайный вектор в En.

Пусть в — равномерно распределенная на отрезке [0,1] случайная величина,

Y = Y (x,t,te),

а

w = #■

Zi

тогда несмещенными оценками ui(x,t) и u2(x,t) будут

a(x,t) = t • W(x — Y,Y,t,te) • f (Y,te) (3.3)

и

6 (x,t)= W(x — Y, Y, t, 0) • p(Y) (3.4)

соответственно. Случайные величины в этих формулах должны быть независимы в совокупности. Если функции f и p ограничены, то дисперсии построенных оценок конечны, поскольку функция W также ограничена в силу неравенства (2.2).

Пусть функции f и p ограничены. Тогда в выражениях для u3 и u4 в формуле (3.1) можно менять порядок интегрирования, поскольку повторные интегралы сходятся абсолютно в силу неравенств (2.2) и (2.10).

Рассмотрим функции

v3(z,\)=[ dT I Q(z, y, X, т)f (у,т)dy, (3.5)

J 0 J En

V4(z,X)=[ Q(z, у, X, 0)p(y)dy. (3.6)

JEn

Тогда

u3(x,t)= dX I Zo(x — z, z,t,X)v3(z,X)dz, (3.7)

0 En

u4(x, t)= dX / Zo(x — z, z,t, X)v4(z, X)dz. (3.8)

0 En

В силу неравенства (2.10) функция V3(z, X) ограничена, а

v4(z, X) ^ const • Xa/<2-1.

Умножая уравнение (2.4) на f (y, т) и интегрируя, получаем для V3 интегральное уравнение

V3(x,t) + f dX f K(x, z,t,X)v3(z, X)dz + f с!т f K(x,y,t,т)f(y,т)dy = 0. (3.9)

0 En 0 En

Аналогично, подставляя в (2.4) т = 0, умножая его на р(у) и интегрируя по у, получаем для «4 уравнение

«4(х,€) + [ <1Л [ К(х,г,Ь,Л)«4(г,Л)йг + [ К(х,у,Ь, 0)р(у)йу = 0. (3.10)

•у 0 ^ Еп ^ Еп

Из неравенства (2.8) следует, что к уравнениям (3.9) и (3.10) применима схема Неймана—Улама [4]. Для ее реализации достаточно выбрать плотность вероятности перехода обрывающейся цепи Маркова, согласованную с ядром уравнения К(х, г, Ь, Л). В качестве такой плотности можно взять, например,

р{{х,Ь) —> (г, А)) = —- А)“/2 х^2(ж - г, £ - Л), (3-11)

где 0 < д < 1 является вероятностью обрыва цепи на текущем шаге, а

г2(х-г,Ь-Х)=(^-^ ехр (3.12)

при 0 ^ Л < Ь. При Л > Ь функция 22(х — г,Ь — Л) =0. Постоянная С в этих формулах берется из неравенства (2.6), которое также влечет согласованность плотности и ядра интегрального уравнения. Отметим, что вероятность обрыва цепи на каждом шаге постоянна, поэтому цепь обрывается с вероятностью 1 и среднее число шагов до обрыва равно д-1.

Для моделирования цепи {(хт,Ьт)}с^=1, стартующей из точки (хо,Ьо) = (х,Ь), будем использовать независимые в совокупности последовательность изотропных

г 1 оо

случайных векторов {^т|т=1, а также последовательность случайных величин {^т 1т=1, имеющих гамма-распределение с параметром п/2 . Предполагается также, что эти случайные элементы независимы с выбранными ранее элементами и, в, 7.

Пусть, кроме того, {вт)т=1 —последовательность случайных величин, независимых в совокупности с определенными ранее, имеющих на отрезке [0,1] бета-распределение с параметрами (1,а/2). Плотность распределения имеет вид

р(Ю = |(1-«Г/2-1 (3.13)

Тогда, при т = 1, 2,...

^т ^171—1^1717 хш хт— 1 \/~С ^(^т—1 (3*14)

Момент обрыва цепи N имеет геометрическое распределение с параметром д, то есть Р{N = т} = д(1—д)т при т = 0,1, 2,..., и независим от траектории. Для построения несмещенных оценок для решения уравнения

ю(х,€) + [ <1Л [ К(х,г,Ь,Л)ю(г,Л)йг + Г(х,Ь) = 0 (3.15)

Jo -'еп

определим последовательность весовых функций Ш(0) = 1,

цг(т) _ ^т+1цг(т-1) ^(Дт-Ь^ш^т-Ь^ш)

Р((хт—1 ,Ьт—1) ^ (хт,Ьт))

при т = 1, 2, . . .

(3.16)

Стандартными несмещенными оценками [4] для функции у(х,1) являются случайные величины

г](х,г) = ^2 ш('т)Р{хт,гт) и ((х,^ = ——_ т=0 д

В качестве функции Г(х,1) в уравнении (3.9) берется

Г(х,1)= [ &г \ К(х, у, I, т)/(у,т)<1у,

^0 л Еп

поэтому несмещенными оценками для «з(х,1) будут

N

П3 (х,Ь) = (1 — д)^2 Ш(т+1)1 (Хт+1,Ьт+1)

Сз(х,Ь)

т=0

(1 — д)Ш ^ +1)/(х„ +1,г„ +1)

д

Полагая У = У(х,г,гв), для функции из(х,г) получаем две несмещенные оценки:

Ь(х,г)= г • Ш(х — У, У,г, гв) • щ(У,гв) (3.17)

и

£з(х,г) = г • Ш(х — У,У,г,гв) • (3(У,гв). (3.18)

В уравнении (3.10) роль функции Г(х,г) выполняет функция

Р1(х,г)=[ К (х,у,г, 0)р(у)3,у, ■/е„

которая неограничена. Поэтому аналогичные несмещенные оценки для V4(х,г) имеют бесконечную дисперсию. Чтобы получить оценки с конечной дисперсией, нужно при оценивании Г(х,г) включить особенность ядра К(х,у,г, 0) в плотность. Для этого запишем V4(х,г) в виде разности V4(х,г) = «5(х,г) — Г1(х,Ь). Тогда для «%(х,г) справедливо уравнение (3.15), в котором функция

Г (х,Ь) = Г2(х,г) = — ( (1Л ( К (х,г,г,Л)[ К (г,у, Л, 0)р(у)с1ус1г (3.19)

^ 0 •/ Еп 'У Еп

уже является ограниченной. Таким образом,

П4(х,г) = — [ <1Л ( 2о(х — г, г,г, Л) ( К (г,у, Л, 0)р(у)ё1у3г + и5(х,г). (3.20)

•у 0 'У Еп Еп

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

у = У(х,г,г/3), уо = У, у1 = уо + у/с-1(г/з)71^,

= и^) = \¥{х-У,У,1,Щ-\¥1-р{у1). (3'21)

а Л2(у1 — уо,гв)

Случайная величина £4 (х,г) —несмещенная оценка интеграла.

и

Функцию и^(х,1) оцениваем на траекториях цепи Маркова. Для этого требуется получить несмещенную оценку Р2(хт,гт). Интегралы (3.19) и (3.20) имеют одинаковую структуру и одинаково оцениваются. Теперь случайная величина в имеет бета-распределение с параметрами (а/2, а/2). Построение оценки завершают формулы

УО = Хт -\- \/~С

VI = Уо + л/С^Ч

т /^7т+2^т+2?

ТТГ _ Д Ґ& 0-\ оІ-а/2/і _ Д',!-# К{хт,уо,Іт,Іт[3) К(уо,уі,Іт0,О) уут- гт& Р и Р) г2(хт-у0,іт-іті3)г2(у0-уі,іті3)’

фт = ^^т^(у1 )•

(3.22)

Здесь фт —несмещенная оценка Р2(хт ,іт), а

Г1

в (!•?)=/ 8“/2-1(і

есть бета-функция. По аналогии с оценками для из(х, г), получаем несмещенные оценки для п$(х, г):

N

&(х,г)= г ■ ш(х - У,У,г,гв) ■ ^ ш(т)фт, (3.23)

т=0

£5(х,г) = г-ш(х-у,у,г,гв) ■ту( )^>дг, (3.24)

ч

которые построены на траекториях цепи Маркова, стартующей из точки (У, Ю).

Почти очевидным следствием включения всех особенностей в плотности вероятностей перехода марковской цепи является конечность дисперсии построенных оценок. Действительно, дисперсии стандартных оценок в схеме Неймана—Улама заведомо конечны [4], если ряд Неймана для ядра К2 /р сходится в том функциональном пространстве, в котором решается интегральное уравнение (3.15). В нашем

(Т)

случае это пространство Ьж(В)п+1). Поскольку

\К\ 2 с /7г\"/2

р ^ а(1 — д) \С7/ 2’

для оператора Вольтерра с ядром К2 /р выполнено неравенство

к'-(х,у,г,т) 6од(,_т)-(-+2-„,,2ехр/ ск^£), ($Ж)

р((х,г) ^ (у, т)) V г - т

которое влечет сходимость ряда Неймана. Отметим, что дисперсии оценок равномерно ограничены при ограниченных / и р.

Таким образом, доказана

Теорема 1. Пусть коэффициенты параболического оператора (1.1) принадлежат классу На,а/2(пПТ+1 ),а < 1, матрица коэффициентов при старших производных симметричная, а ее собственные числа лежат в фиксированном отрезке [V, р]

и V > 0. Пусть функция /(х,г) из класса На,а/2 (0^++!) и непрерывная функция р(х)

ограничены. Тогда случайные величины £ = £1 +£2 +£з+£4 +£5 и £' = £1 +£2 +Сз+£4 +£5, определяемые формулами (3.3), (3.4), (3.17), (3.18), (3.21)-(3.24), являются несмещенными оценками решения задачи Коши (1.2). Дисперсии оценок равномерно ограничены по переменным х, Ь.

4. Определение постоянных с и С. Вычислив производные в представлении (2.5) ядра К, получим

К(х,у,г,т) = ———^—-[п - Тг(А(х,г)А-1(у,т))}го+

2(Ь - т)

+ 4^~~^2 (х -У)'А~1(У’Т)[А(У’Т) - а(х,*)1а~1(у,т)(х - у)г0-

- 7^7-—7а,(х,г)А~1(у,т)(х - у)го + а0(х,1)г0, (4.1)

2(Ь - т)

где функция ТГ(-) —вычисляет след матрицы, а(х,Ь) —столбец с компонентами а^(х, Ь), а операция а' означает транспонирование матрицы. В качестве С можно взять любую постоянную, удовлетворяющую неравенству 4рС < 1.

Пусть А,В — симметричные матрицы с постоянными элементами, удовлетворяющими условиям

|а^-(х,Ь) - а^- (у, т)| < а^- |х - у\а + Ъм \Ь - т \а/2,

а Б — ортогональная матрица, приводящая матрицу А-1(у,т) к диагональному виду. Тогда

п - Тг(А(х, Ь)А-1 (у, т)) = Тг(Б'(А(у, т) - А(х, Ь))ББ'А-1 (у, т)Б) =

п п п п

X-Б'(А(у,т) - А(х,Ь))Б- = ^2^2 "^2 Йк’3(ак,т(у,т) - ак,т(х,Ь))вт-,

3=1 3=1 к=1 т=1

где Б- — у-й столбец матрицы Б, а X- — у-е собственное число матрицы А-1 (у, т). Отсюда получаем неравенства

1 п п

|п - Тт(А(х,г)А-1 (у, т))| < \ак,т(у,т) ак,т(х,^')')\ <

к=1 т=1

Л п п

< - ^2 ^2 (ак,т\х - у\а + Ък,т\Ь - Т|“/2) = сЛ\х - у\а + с2\г - т\а/2. (4.2)

Несложные выкладки дают следующие оценки сверху для квадратичной формы в равенстве (4.1):

(х - у)'А-1(у,т)[А(у,т) - А(х,г)]А-1 (у,т)(х - у) <

п п

<^^2 \Хк\ (а-к,т\х - у\а + Ък,т\г - т\а/2) Ы <

к=1 т=1

< (||А|| \х - у\а + ||В|| \Ь - т\а/2) \г\2 <

<

(а \х-у\а+ В \г -г|“/2) (^-\х -у|^ =

= (сз\х - у\а + 'с-4\Ь - т\а/2)\х - у\2, (4.3)

где Хк —компоненты вектора г = А 1(у,т)(х - у).

Оценка билинейной формы в равенстве (4.1) имеет вид

|а'(х,г)А г(у, т)(ж - у)\ < ^\а(х,Щх - у\ < ^-\х - у\ = с5\х - у\, (4.4)

где а — вектор, компонеты которого оценивают компоненты вектора а(х,Ь), то есть \а^(х, Ь)\ < а при г = 1, 2,... ,п. Из неравенств (4.2)-(4.4) получаем оценку сверху для \К\:

\к(х,у,г,т)\ < _ (с!\х - У\а + с2\г - т\а/2^ г0+

1

2{Ь - т)

+ щ Ых ~ у\а + с4^ — т\а/2)\х — у|2^0 +

+ Х7Г~—\с5|ж - у\%о + а0г0, (4.5)

2(Ь - т)

где ао —верхняя граница для \ао(х,Ь)\. Запишем это неравенство в виде \К(х,у,Ь,т)\ < с(х,Ь,у,т)(Ь - т)а/2-12о,

где

, . 1 \х - у\а „\ \х - у\2 / \х - у\а „\

с(жЛу’т) = 2 + С2) + \Сз\Г^\^ + СА) +

+ г52|1 [1/2 (* ~ т^аУ2 + ^ ~ т^аП- (4'6)

г0 (х-у,у,г,т) ^ ехр(-С\^~—п/2ехр(-С———

Применяя к Zо неравенство (2.2), получаем

Ш) ' ■-

где

С1 = ^~С-

Пусть в > 0. Функция вв вхр(-Св) ограничена на интервале [0, +то) постоянной

М((3,С)=^^ ехр(—/3),

поэтому

\х - у\2 А 1 (а

с(ж,г,г/,т)ехр ^-С\ ^ ^ ^ [с.\М (7^1) + сг) +

+ ± (с3М (1 + 1,^) +^(1,(7!)) +с5^М ^С^Т1-^2 +а0Т1-а/2 = с.

Отсюда получаем окончательное неравенство

Щх,у,г,т)| < ^-г)-("+2-а)/2ехР ) .

(Т )

5. Оценка функционалов. Пусть Н(х,1) —интегрируемая на РП+1 по мере Лебега функция. Для оценки функционала

Ф(к) = & к(х,1)п(х,1)ё.х

■)о

обычно выбирают плотность распределения начальной точки п(х,Ь) и какую-либо несмещенную оценку £(х, I) решения. Тогда величина

/г(жо^о)£(жо^о)

п(хо ,и)

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

Ъ?(х, £) п(х, I)

(Т

интегрируема на РП+1 по мере Лебега, а в качестве £(хо,1о) выбрана одна из ранее построенных несмещенных оценок для п(х,Ь).

Для вычисления функционала Ф(к) можно использовать также «сопряженную» схему, предложенную в [5] для решения начально-краевой задачи. Для такой схемы, ядро Zо можно включить в плотность вероятности перехода, что может привести к уменьшению дисперсии. Отметим, что реализация «сопряженной» схемы не требует вычисления констант с и С.

Автор благодарен проф. С. М. Ермакову за полезные обсуждения проблемы и постоянное внимание к работе.

Литература

1. Ладыженская О. А., Солонников В. А., Уральцева Н. Н. Линейные и квазилинейные уравнения параболического типа. М.: Наука, 1967.

2. Сабельфельд К. К. Методы Монте-Карло в краевых задачах. Новосибирск, Наука,

1989.

3. Симонов Н. А. Стохастические итерационные методы решения уравнений параболического типа // Сиб. матем. журн. Т. 38, №5. 1997. С. 1146-1162.

4. Ермаков С. М., Михайлов Г. А. Статистическое моделирование. М.: Наука, 1982.

5. Сипин А. С. Блуждание по цилиндрам для параболических уравнений // Математические модели. Теория и приложения / под ред. проф. М. К. Чиркова. Вып. 11. СПб.: Изд-во ВВМ, 2010. С. 83-103.

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