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

УДК 519.71

Вестник СПбГУ. Сер. 1. 2012. Вып. 1

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

А. С. Сипин

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

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

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

( д д\ д п д2 п д ь = ь Н м Ж = Ъь - £ + + ао(-хЛ (11)

\ / » з 1 »

коэффициенты которого принадлежат классу На>а/2{В^1), а < 1. Матрица коэффициентов при старших производных предполагается симметричной, а ее собственные числа лежат в фиксированном отрезке [г/, /х] и V > 0. Рассмотрим задачу Коши

( д д \

М Х,1,дх,д11и = ^ и^=0 = ¥(х). (1.2)

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

и(х,г) = (1т Z(x,y,t,т)f(y,т)dy + / г(х,у,г,0)1р(у)(1у. (1.3)

■) о -!еп -!еп

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

Пусть /г.(ж, ¿)—интегрируемая на по мере Лебега функция. Для оценки

функционала

Ф(к) = / Л Ъ,{х,1)и{х,1)<1х (1.4)

ио J Еп

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

7г(жо^о)

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

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

7г(ж, ¿)

интегрируема на -0„+1 по мере Лебега, а в качестве £(жо,£о) выбрана оценка для и{х,Ь) с конечной дисперсией [2]. Такие оценки функционалов рассмотрены в работе [3]. Там же имеется подробный обзор статистических методов решения задачи Коши.

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

Фундаментальное решение, в свою очередь, является функционалом от решения <3 некоторого интегрального уравнения Вольтерра [1], к которому применима схема Неймана—Улама [2]. Таким образом, речь идет о построении несмещенных оценок функционалов от <5. В работе [3] для этого использована «прямая» схема построения несмещенных оценок от решений интегральных уравнений, а в данной работе — «сопряженная» схема.

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

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

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

1

г0(х - у,у,г,т)

[4тг(* - т)]п/2 (с1е^(у,т))1/2 'Мг-т

ехР ( Е ] . (2.1)

Ь3=1

При £ < т функция Zo(x — у, г) = 0. Функция Zo удовлетворяет неравенству

г0(х -у,у,г,т) < [--) г^х-у^-т), (2.2)

где

/ - у I

гг(х -у,г -т) = ----—т-ехр -

[4тф(* - т)]п/2 V М* - т) Фундаментальное решение Z можно представить в виде

г(х,у,г,т) = г0(х - у,у,г,т) + ¿X ! г0{х - г, г, г, Х)(^(г,у, Х,т)йг, (2.3)

■'т .] Еп

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

<Э(х,у,г,т)+ [ ¿X [ К(х,г,г,\)С}(г,у,\,т)<1г +К(х,у,г,т) = 0, (2.4)

■1т ■) Еп

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

( \ \ \—л г / л \ / м^ (ж — Л) К{х,г,г,Х) = - ---Ь

Е, \ (ж --А) . \ Г7 I Л \ /Г» ^ \

аДх, ----Ь а0(х, £)А)(х - Л). (2.5)

г=1

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

\К(х,у,1,т)\ < ф -тГ(п+2-а)/2ехр (-сЦ^.) . (2.6)

t-т г

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

Km(x,y,t,T)= dX К(х, z,t, \)Km-i(z, у, А, г)dz, m = 2, 3,... (2.7)

Jr Jfî,!

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

!«.<.,< С" (§)""-"" - >/'ехр (-ci^f ) , (2.8)

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

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

оо

Q(x,y,t,T) = Y,(-l)mKm(x,y,t,T), (2.9)

m= 1

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

|д(ж, у, t, r)| < Cl(t - r)-("+2-«)/2 exp f-C^pj . (2.10)

3. Оценка функционалов. Для вычисления функционала Ф(/г) (1.4) будем использовать «сопряженную» схему Неймана—Улама, которая позволяет включить в плотность вероятности перехода ядро Zo, что может привести к уменьшению дисперсии. Для этого, используя (1.4), (1.3) и (2.3), запишем функционал в виде

Ф {h) = {h) + Ф 2{h) + Ф3(Ь) + Ф 4(h) =

= dt dxh(x,t) dr Zq(x — y, y, t, т)f(y, т)dy+

J 0 J En J 0 J En

+ сМ ¿х!г(х^) / - у,у,Ь,0)ч>(у)<1у+

J0 -1 Еп -1 Еп

+ / <М (1х1г(х^) ¿т / (¿Л / Zo(x — г, X)Q(z,y, X,т)dz \ f(y,т)dy+

■1 о -1 Еп Jo -1 Еп \7т Еп )

(1х1г(х^) / I ¿X Zo(x — z,z,t,X)Q(z,y,X,0)dz \ íp(y)dy. (3.1)

3 Еп \ио иЕп /

[ Л [

/о и Е„

Заметим, что функция Zо(х — у,у,Ь,т) по переменной х является плотностью распределения нормального случайного вектора X, у которого матрица ковариаций 2(4 — т)А(у,т), а среднее равно у. Такой вектор моделируется по формуле

Х = у+ у/2(1-т)у/А(У,т)Е,

где \[А — треугольная матрица квадратного корня из матрицы А, 2 — нормальный случайный вектор с единичной матрицей ковариаций и нулевым средним.

Предполагая ср(у) интегрируемой в Еп по мере Лебега, а к(х, 4) ограниченной, получаем несмещенную оценку для $2(/г):

Ф2 = Т^Щ-к (у + ^2Тву/А(У, Тв)Е, Тв) , (3.2)

7Г2(У) V /

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

вектор У имеет плотность распределения 7Г2(у), согласованную с функцией у>(у).

(т)

Предполагая /(у,т) интегрируемой на -0„+1 по мере Лебега, меняя порядок интегрирования по переменным 4 и г, аналогично получаем несмещенную оценку для

Ф1 (ну.

Ф1 = (Т- (у + (Т - в)в^А(У, (Т - 8)0)5, © + (Т - 8)0) , (3.3)

где случайный вектор (У,©) имеет плотность распределения тт^ (у, т), согласованную с функцией /(у,г).

Для оценивания Фз(/г) сделаем замену переменной и представим его в виде

Фз(/г) = / <1т <1у/(у,т) / ¿X ¿г сМх ■) о ¿т ./а

х / Zo(x — z,z,t,X)Q(z,y,X,т)h(x,t)dx. (3.4) Используя (2.9), получаем следующие несмещенные оценки для С^(г,у,Х,т)\

N , чт-1

= -£(- — ) Кт(г, у, А, г), (3.5)

т=1^

V/ = - (-^Г (3.6)

V 1-9/ ч

Здесь 0 < </ < 1, а N имеет геометрическое распределение, то есть при то = 1,2,... Р(М = то) = </(1 — д)т-1. Подставляя оценки (3.5) и (3.6) в (3.4), получаем несмещенные оценки для Фз(/г), которые подлежат дальнейшей рандомизации. Для этого нужно выполнить несмещенное оценивание функционалов

= [ ¿т [ ¿у/{у,т) [ ¿X [ ¿г [ (Их

.] Еп -1т .] Еп

х / Zo(x — г, Х)Кт(г,у, X,т)h(x,t)dx. (3.7)

Еп

Воспользуемся процедурой оценивания, предложенной в [6] для решения начально-

(т)

краевой задачи. Пусть А) —ограниченная функция в Рассмотрим интеграл

У(у,т)=[ д,Х ( Х)К(г, у, А, т)д,г. (3.8)

■)т ■) Е„

г-Т

¿X

1е.

Пусть

^0(у,г) = {Ше Еп\итА-\у,т)и = 1}

— эллипсоид с центром в нуле,

_ 27т"/2

~ Щ)

— площадь поверхности сферы радиуса 1 в Еп. Случайный вектор П распределен на 8®(у,т) с плотностью

р(у,т,ш) =-. -. (3.9)

Его можно моделировать по формуле

П = у/А(у,т)П ь (3.10)

где —изотропный единичный вектор в пространстве [7]. Тогда /■Т

/1 со

¿X у ¿гх

„ ,'п-Тг (А(у + Н1,Х)А-1(у,т)) ( г2 1 п_х ,

х Е ---1—-—-——--—--- ехр----- г у(у + г\1, А) —

^ (А-т)Г(§)(4(А-т))* Ч 4(А -г)] >)

- (Т (IX Г (1 Е2Г" т)А[-У + гП> ~

I Уо Г 4(А-т)2Г(|)(4(А-т))§ Х

х ехр (~4(а"Т7) ) Г"~1у<уУ + А')~

[Т „ Г , ^ {га<(у + гП,Х)А-1(у,т)П ( г2 N „ < , ^ Л - ¿X ¿гЕ [ — „ ' ехр----- )гп~1у(у + гП,Х +

У- Уо V (А-т)Г(|)(4(А-т))5 Р1 4(А-г); ^ ' V

Гт Г°° 2гп~1 ( г2 \

+ У>Уо ЛгЕа°{У + ^ ХНУ + ^ Л)Г(| )(4(Л - т))» еХР {-ЦХ^Г)) ' (311)

где знак Е означает математическое ожидание (по С1). Из условия принадлежности коэффициентов уравнения классам Гёльдера вытекает следующее представление для выражений в первых двух интегралах:

п-Тг (А(у + Н1, Х)А-\у, т)) = Тг ([.А(у, г) - А(у + гП, г)] А-1 (у, г)) + + Тг ([А(у + гП,т) - А(у + гП,Х)}А-\у,т)) =

= д1(у + гП, у, т)га +д2(у + гП, у, X, т)(А - т)а'\ (3.12)

[ПтА-\у, т)А(у + гП, Х)А-\у, т)П - 1] =

= ПтА-\у, т) [А(у + гП, т) - А(у, г)] А^1(у, т)П+ + ПтА-\у, т) [А(у + гП, А) - А(у + гП, г)] А'1 (у, т)П =

где <?1, д2, Ъ,\, к2 суть ограниченные функции. Подставляя эти выражения в (3.11) и выполняя замену переменной

/и(у + гП, у, т)га + к2(у + гП, у, А, т)(А - т)а/\ (3.13)

Г

5 4(А — г)'

получаем для У(у, т) представление

х Е (л(у + 2у/з{\-т)П, у, т)ь(у + 2у/з(Х-т)П, А))

х Е (д2{у + 2\/8{Х - т)П, у, А, т)ь(у + 2у/з{Х - т)П, А))

- / ¿А(А-г)2 1 / ¿в——1 ехр( —в)х

мт Jo 1(2)

¿Х(Х -г) 2 ^цпу8

х Е (ьл(у + 2у/з(Х-т)П, у, г)«(у + 2у/з(Х-т)П, А))

х Е (к2{у + 2у/з(Х-т)П, у, А, г)«(у + 2у/з(Х-т)П, А))

х Е (а'(у + 2у/з(Х-т)П, X)А-\у, т)Пу(у + 2у/з(Х-т)П, А)) +

х Е (а0(у + 2у/з{Х-т)П, X)«(у + 2у/з{Х-т)П, А)) . (3.14)

Пусть 7(то) —случайная величина, имеющая гамма распределение с параметром то. Тогда представление (3.14) можно записать в виде

У(у, г) = - £ ¿Х(Х - гГ^Е» ^ + 2^7 (А ~)П, У, Г ^ х

х « ^ + А) ^ -

~ ^ ¿Х(Х-т^-1^Ед2 (у + ^1^) (А -т)П,у,\,т^ х х у(у + 2^7 (А-т)П, А)) -

х « (у + 2^7 (А "о, А^ -

" £ ^А(А - ^ + 2^7 (А А^ г)Пх

х « (у + 2^7 (А А^ +

+ £ ЙА£а0 (у + 2^7 (А-т)П, + 2^7 (А - т)П, а) . (3.15)

Оценивая интегралы по переменной А, получаем несмещенную оценку г/(у, т) для У{у,т):

г,(у,т) = -(Т-т)* + (Т~т)Ж,у,т^ х

х ^у + (Г " г + (Т - т)^ -

_(Т_т)§132 ^ + 2^7 (Г " у, г + (Т - г)!?, т) х х V [у + 2^7 {Т-т)Ш, т + {Т- тЩ)^ -

^2а+1Г(п+2+а^ ( ^ /_ (п+ 2 + а

-% Ь+ 24/7 1 о ) (Т-т)№,у,т) х

ху[у + 2\Н[ п + 2 + а )(Т-т)т,т + (Т-тЩ ] -

-(Т-тр^ х к2 + ) (Т-т)№,у,т+(Т-т)#,т

XV I у + 2,/7 ( ) (Т-т)№,т+(Т-т)д I -

" (Т - (у + 2^7 ) (Т - т)ёП, т + (Т - т)ё ] А^1{у, т)Пх

ху[у + 2хН[ ^ ) (Т - т)6П, т+(Т-т)б\ +

+ (Т- т)а0 [у + 2^7 (Т-т)вП, т + (Т - т)в) х

х ь(у + 2^7 (Т - т)Ш, т + (Т - т)^ , (3.16)

где случайные величины (5, в распределены на отрезке [0,1]. Величины (5 имеют плотности а/2ва/2-1 и 1/(2л/в) соответственно, а в распределена равномерно. Выбирая с вероятностью 1/6 одно из слагаемых в (3.16) и умножая его на 6, получаем окончательную несмещенную оценку С(Утт) Для У(Утт)-

Оценки фт для функционалов Фт(/г) из (3.7) будем строить на траекториях

неоднородной цепи Маркова {(ук, Тк)}^о- Начальное состояние цепи выбирается в

(т)

с некоторой плотностью 7г(у, т), согласованной с функцией /(у, г). В оценку фт, при этом записывается дробь

Кур, то) 7г(з/о,т0)'

Далее, последовательно выполняются то переходов, на каждом из которых реализуется оценка С{ук,тк) по независимым в совокупности последовательностям случайных элементов {^к}^0, {^Л^о- При этом фт умножается на весовой

множитель в оценке С(Ук, т~к), & аргументы функции V определяют следующее состояние цепи. Например, если на шаге к = 1, 2,..., то в оценке (3.16) было выбрано первое слагаемое, то фт нужно умножить на

2ар( п+а ^ / / / о,\

-6(Г-гА;_1)а/2 91 {ук-1+2]11к-1 (—-—^ (Т - тк-1)'&к-1Пк-1,ук-1,тк-1

и перейти в точку (ук, тк) с координатами

ТЪ (У,

Ук = Ук-1 + 2«/7*_1 ( —-— ) (Т - Тк-^к-^к-!,

Тк = тк-1 + (Т - тк-1)$к-1-

На заключительном этапе, в соответствии с (3.7), выбираются время распределенное равномерно на отрезке [тт,Т], и координата хт, распределенная нормально с плотностью 2о(х — ут,ут^т,тт), а оценка фт умножается на (Т — 4т)/г(жт,¿т). Применяя формулы (3.4)—(3.7), получаем оценки функционала Фз(/г):

N , ^ \ т—1 т=1 ^ ^ '

Функционал $4(/г) оценивается также по формулам (3.17)—(3.18), в которых величины фт строятся на траекториях марковской цепи, стартующей из точки (г/о, 0). Первоначальное значение оценки фт равно

Аур)

тг(УО)'

где уо распределено в Еп с плотностью 7г(у), согласованной с функцией у (у).

Ограниченность дисперсий построенных оценок вытекает из ограниченности весового множителя в оценке С(у, т) и доказывается так же, как аналогичное утверждение в «прямой» схеме для первой краевой задачи [6].

Теорема 1. Пусть функция ограничена, а функция /2(у,т)/п(у,т) — ин-

(т)

тегрируема в -0„+1 = Еп х (0, Т). Тогда дисперсии оценок фз и ф'3 конечны. Доказательство. Оценим второй момент для ф'3:

00 ТГ I 2

ЭД)2 - Е я{1 ^ (3.19)

т=1

Заметим, что

ЕФ1 = Е^^У^П^СЧу^Т -

где С(Ук, Тк) —весовой множитель в оценке С(Ук, Тк). Рассмотрим интегральный оператор, определенный представлением (3.14), в котором все интегралы берутся со знаком «+», а функции </1, <72, /11, Н.2, ао и выражение а'(у + 2у/в(А — т)Г2, А)А~1(у, т)Г2 заменяются на их модули. Пусть К(х, у, г) — ядро этого оператора, а Кт(х, у, г) — его повторное ядро. Для ядра К(х, у, г) справедлива оценка (2.6) с некоторыми новыми постоянными с\ и Сь поэтому ряд из повторных ядер

оо

]Г ВтКт{х,у,±,т) (3.20)

т= 1

равномерно сходится для любой постоянной В. Весовой множитель С(Ук,Тк) в оценке С(Ук,Тк) ограничен в силу ограниченности функций <71, <72, /и, ао и выражения а'(у + 2у/в(А — т)Г2, А)А~1(у, т)С1. Поэтому, существует постоянная В, такая что \С(Ук,тк)\ < (1 -<7)В. Тогда,

Еф2т < [ <1т [ ¿у^ [ (IX [ сЬ [ (Их

Jo -1еп ъ{у,т) .1т .)вп Jx

2/

х/ г0(х--д)тВтКт(г,у,\т)Н2(х,г)ах. (3.21)

' Еп

Теперь, из равенства (3.19) следует неравенство

■т , /2(у>т) гт г

Е(ф'3)2 < [ <1т [ ¿у

■1О ■) Е„

О JEn ъ{у,Т, Jт л3п JX

¿А ¿г сИх

XV/ г0(х-г,г,г, Х)(^(г,у,Х,т)ТН2(х,г)(1х < оо, (3.22) 1 Зеп

где <3 — сумма ряда (3.20).

Для второго момента случайной величины фз

оо г? / 2 00 1/1 00

Е Ш- (3-23)

т= 1 Ч> т= 1 Ч> г=т+1

Уже доказано, что первое слагаемое в (3.23) конечно. Пусть Зт — ст-алгебра, порожденная цепью до момента времени то, а функция (у,т) определена равенством

(Утт) = / «¿А / (А с!г ! Zo{x — z,z,t,\)Q{z,y,\,т)\h{x,t)\dx^,

>1 т Jx J Еп J Еп

тогда

Е( Е ) = 1/г(У0'Т0?1ПГ=-01|СЫ,п)Ыут,тт).

\1=т+1 / ^УУО^ТО)

Теперь, аналогично (3.21) и (3.22), получаем неравенство

(•Т с -г2/„ лТ л сТ

т=1 г=т+1 ]Е

. (1у ^ ^ / (¿А / (£г / сЙх

1еп к{у,т) }т }Еп 7а

х А)(1 — </) / — А)<5(г;, у, А, г)|/1(ж, < оо

Еп

Отметим, что реализация «сопряженной» схемы не требует вычисления констант с и С.

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

Литература

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

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

3. Сипин А. С. Статистические алгоритмы решения задачи Коши для параболических уравнений второго порядка // Вестн. С.-Петерб. ун-та. Сер. 1. 2011. Вып. 3. С. 65-74.

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

1989.

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

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

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

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