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

Вычислительные технологии

Том 13, Специальный выпуск 4, 2008

Модификация случайных блужданий для оценки градиента решения эллиптических краевых задач*

A.B. Бурмистров Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия Новосибирский государственный университет, Россия e-mail: burm@osmf . sscc. ru

An interior boundary-value problem for an elliptic equation is considered. The equivalent system of local integral equations is constructed with the help of the Green function for the Poisson—Boltzmann operator. The solution and its gradient are estimated using the "random walk on spheres and balls" algorithm.

1. Дифференциальная задача и определения

Рассмотрим смешанную краевую задачу:

Ди(г) - c(r)u(r) + (v(r), Vu(r)) = —g(r), r e П С R3; (1)

u(y) = Фг(y), y e Г1; (2)

(Vu(y),Y(y)) + а(У)и(У) = ф2(У), У e Г2j (3)

внутри области П с границей Г = Гг (J Г2, которая предполагается кусочно-гладкой. Будем считать, что коэффициенты уравнения (1), граничные функции и область П таковы, что существует единственное достаточно гладкое решение исходной задачи. Область может быть неограниченной, в этом случае потребуем от решения u(r) следующего поведения на бесконечности: u(r) = O(r-1) при |r| ^ то, Нас будет интересовать оценка решения u(r) и градиента Vu(r) в некоторой точке r e П. Введем следующие обозначения:

Дк = Д — к2 — стационарный диффузионный оператор. Перепишем уравнение (1), оставив только этот оператор в левой части:

Д^(г) = Д^г) — K2u(r) = — к2 — c(r))u(r) — (v(r), Vu(r)) — g(r); (4)

R = R(r) = minjdist(r, Г), Rmax};

/ C0R2(r)\

Co — константа такая, что p(r) =11--I G (0,1);

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 08-01-00334).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

В(го, R) = {г G R3 : |Го — г| < Д(го)} — шар с центром в го, целиком лежащий в S(r0, R) = dB(r0, R) = {r e R3 : |r0 — r| = R(r0)} — соответствующая сфера; г = {r e П : 3r' e Г, |r — r '| < e} — е-окрестпость границы Г; ш — единичный вектор направления v(r), т, е, v(r) = |v(r)|w(r) аш = cos(r — r ',ш); —c* = n2/Rmax — первое собственное число оператора Лапласа в шаре радиуса Rmax; r* — проекция точки r e Г1е на Г^ |r — r*| = R(r), Условие (2) аппроксимируем в Г1е следующим образом: u(r) = ^1(r*) + O(e) « ^1(r*);

dY(r) - расстояние от точки r до границы Г2 вдоль векторного поля 7; n(r) — проекция точки r e Г2е на Г2 вдоль поля 7, Условие (3) аппроксимируем в Г2е следующим образом [1]:

. . [1 + adY(r)]u(r — £7) e-02(n(r)) 2. . . . . ~ . . ..

м r = L 1 , '):) V, V " + 1 , ¡\ ; , + 0{е2) «prur - ет + Ь 7Г г ; 1 + a(dY (r) + e) 1 + a(dY (r) + e)

Pn(cos 0) — полиномы Лежандра [2];

гт

vn(x) = \ —In+i(x), где In+i(x) — модифицированные функции Бесселя [2];

у 2 2

sinh(x) 1 / sinh(x)

щ(х) =-, /^i(x) = — cosh(x) —

/V» /V» \ /V»

Ду •Х' \

Ор:(г, г') — нецентральная функция Грина для оператора Дк в шаре В(г0,Л), которая является решением следующей задачи Дирихле:

д^(г, г') = ¿(г - г'), (г, г')|рея(Р0>д) = 0,

где $(•) — дельта-функция Дирака, Явные выражения для функции Gр0 и ее производных, используемые для решения поставленной задачи, выписаны далее.

2. Сведение к интегральному уравнению

Формально интегральное уравнение для функции и (г ') из (4) и ее пространственной „ дм '

производной —— (г ) можно получить, используя теорему о среднем, В нашем случае вое-д^

пользуемся нецентральной функцией Грина &р0(г, г') для оператора Д к в шаре В(г0, Л), В результате по аналогии с [3] для г ', г0 € П \ Г£ будут иметь место представления:

щ

dn

-(г) = - I — I I (г, г lu^rjdSr + I -^(r,r)T(u,u),r)dr. (6)

м^г') = — J -j^ijr.r^u^^dSr + J Qp0{r,r')T{u,u,r)dr\ (5)

S(ro ,R) B(ro,R)

дш

S(ro ,R) B(ro ,R)

Здесь F(u1 ,w,r) =

du

+ (к2 - c(r))uAr) + <jf(r)

дш

внешняя нормаль к

сфере 5(г0,Л), Для г ' € Г£ полагавм и1 = м. Далее будут сформулированы условия существования и единственности системы уравнений (5) и (6) и совпадения их с решением исходной дифференциальной задачи и его пространственной производной. Теперь

выпишем функции, использованные в представлениях (5) и (6), а также найдем плотности вероятностей, пропорциональные им. Нецентральную функцию Грина GK0 (г, г') будем искать в следующем виде:

(г г') = — УтоК1-*1 ) 4?г

sinh {k(R — |r — r'|)} sinh {kR} |r — r'|

— G(r, r0,

(7)

где функция G(г, г') является решением следующей краевой задачи:

Д«£(г, r') = 0, |r — ro| < R, G(r, r')

1

4n

sinh {k(R — |r — r'|)} sinh {kR} |r — r'|

|r — r0|

R. (8)

Решение краевой задачи (8) можно записать в сферических координатах [2]:

Г( >\ V^ ( D f*\UMr ~ ðРG(r,r) = J2on(r)Pn(0 Un{KR)

n=0

где £ = cos 9 — косинус угла между векторами (r — r0) и (r' — r0), а коэффициенты ^n(r') определяются по формуле

£n(r')

2n + 1

1

4п

-i

sinh j k(R —

sinh {кЩ sjR2 + \r> - r0|2 - 2R\r' - r0|£

В дальнейшем мы будем использовать в основном центральную функцию Грина для оператора Дк, т.е. значения функции u1 (5) и те градиента (6) в точке r' будем представлять в виде суммы интегралов по сфере и шарам с центром в той же точке (r0 = r'). Отметим, однако, что для построения производных была использована именно нецентральная функция Грина (7), кроме этого, нецентральные функции можно использовать при решении внешних краевых задач. При этом для ограниченности на бесконечности в выражениях для функций vn(x) модифицированные функции Бесселя In+i(x) нужно заменить на функции Kn+i{x) [2], После дифференцирования, обозначая р = |r — r'|, />(р) = крcosh{K(R — р)} + sinh{K(R — р)}, а также учитывая, что полиномы Лежандра ортогональны и P^cos 9) = cos 9 [2], получим

sinh {k(R — р)} d£V Sr'r') = --\ r p J , (r> r)

dnr

(r, r') = —

4пр sinh {kR}

11

4n sinh{KR}

4nR2 v0(KR)' dw \ dn.

dw

д (<9£K' К A аш

(r,r ) = —

/>(р) к ^1(кр) ~ RV1(KR)_

k

4nR2 vi(KR)'

Заметим, что при к ^ 0 оператор А к стремится к оператору Лапласа А и вышеприведенные функции стремятся к соответствующим центральным шаровым функциям для оператора Лапласа [3]:

4п \ |r — r'|

дбг', п

\J I Or

4nR2

a7(r'r) д_

dw

1

dG0'

dnr

4п \ |r — r'|2 йм 3

|r — r'| R3

4nR2 R

a

м

1

Полученные функции пропорциональны плотности распределения на сфере FS и плотностям распределения в шаре F0 и F^

r2 dGк 3 R

G?(r,r') = —C01(KR)F0(ry), i^(r,r') = fl.-CiiHFiir,/), 6 дш 4

dGK d fdGк\ 3

^(гУ) = Coo(KR)Fs(r,r'), -— (r,r') =

our' дш \ on,rU R

Здесь функции Ckl(-), k,l £ {0,1} монотонно убывают от 1 до 0:

Coo(x) = vo(x) , Cio(x) =

x

Coi(x) — ———, Сц(х) — -x2 v0(x) 3

to(x) -

3v1(x)

tp(x) - 1 xv1(x) J '

to(x)

tanh(x/2)

x/2

Функция FS является плотностью равномерного распределения на сфере:

sin(0)d0d^

Fs (tp,e)dS =

4п

здесь 9 £ (0,п), р £ (0, 2п) и р = |г — г'| £ (0,Д) — координаты локальной (с центром в точке г') сферической системы координат. Плотности Fo и в этой системе координат факторизуютея: Fj• (г,г')в,г = FS(р,9)dрd9Fj(р)dр, ] = 0,1, а части Fp(р) имеют следующий вид:

Fq (р)

бСр-^/сД) Д2 sinh {kR}

рsinh {k(R — р)} , Fp(p)

Д2 sinh{^}

/>(р) —

кр2 ^1(кр) Д 1/1 (кД)

Моделирование случайных величин с такими плотностями можно осуществлять методом исключения [4].

Для построения статистических алгоритмов решения системы уравнений (5) и (6) объединим их далее в одно интегроалгебраическое уравнение с помощью специального расширения множества фазовых координат дискретной переменной принимающей только два значения: ] = 0 или ] = 1, Обозначим точку нового фазового пространства через w = (г,]) £ R3 х {0,1} и введем также следующее обозначение:

П(г),

U(w) = U(г, J) =t R(r) On,

3 дш

J = 0; (г), J = 1

1 — CoR2(r)

Используем величину руг) =-для рандомизации полученного уравнения:

6

U (г', J') = р(г')

Fs (г, r')U (г, 0)Qjo(r, r')dSr +

S(r',R)

G(r',f)

р(г')

+

+(1 — р(г')) J Fj(г, г') [plQ^^U(г, 1)+ pQQ^y^U(г, 0)] dг, (9)

B(r',R)

где р0 = р0 (г, г') и = р1(г, г') тоже являются некоторыми вероятностями Р0 + Р1 = 1 которые можно взять равными р1 = р0 = 0.5, однако более эффективно брать их пропорциональными функциям и |к2 — с|:

с 3Кг)| с |к2 — с(г)| 3|^(г)| |2 ...

^ = -£ПГ-> = - ' где Р™ = + к - с г •

Я(г)рст раи я(г)

Функция С^) представима в виде интегралов [Си = Сы(кД(г')); к,/ € {0,1}]: С(г/,0) = СО1^ I Р0{г,г')д{г)(1г, С(г,,1) = С11^ I Р^г'Ыф^г,

В(г' ,Я) В(г' ,Я)

которые можно оценивать "по одному случайному узлу" [4]. Веса выражаются так: /о (V - пс (г г') 3|^(г)|С<01 Ос Гг г') 9аХг)1С"

о (г г<\ 0е Гг гЧ - ~ Ф))Со1 0е , А _ 3^(^-с(г))Сп

3. Алгоритм решения интегрального уравнения

В соответствии с (9) и с учетом введенных обозначений статистический алгоритм сопряженных блужданий по сферам и в шарах для оценки величины и ^0) строится следующим образом,

= (г0,^0), 00 = 1; траекторию моделируем до обрыва на шаге со случайным номером Щ.

• Если гп € Г£, то следующая точка гга+1 выбирается на сфере Б(гп) с вероятностью р(гп) (случай 1) или в шаре В(гп) с дополнительной вероятностью 1 — р(гп) (случай 2), Для случая 1 точка гга+1 выбирается равномерно на сфере Б(гп), т.е. по плотности ^(гп, г'). В счетчик добавляется величина 0гар-1(гга)С^га). Величина полагается равной 0. Вес умножается на величину 0,п0: 0га+1 = 0га0,п0.

Для случая 2 точка гга+1 выбирается в шаре В(гп) то плотности (гп,г'). Величина полагается равной 1 с вероятностью р1(гп, гга+1) и равной 0 с вероятностью р0(гп, гп+1). Вес умножается на вел ичину 0га+1 = •

Продолжаем алгоритм от места, обозначенного •, для точки гга+1. Если гп € Г1е, ^ = 0, то в счетчик добавляется величина 0га^1(гП) и траектория обрывается. Если гп € Г2е, ^ = 0 т0 в счетчик добавляется величина 0га^/'2(п(гга)). С вероятностью р(гп) следующая точка гга+1 = гп — ега7га отражается в О \ Г£. Величина полагается равной 0. Вес не изменяется. Продолжаем алгоритм от места, обозначенного • , для точ ки гп+1. С вероятностью 1 — р(гп) траектория обрывается. Если гп € Г£, ^ = 1, то траектория обрывается, В результате для величины и^0) получаем реализуемую оценку:

N

СЮ = £ О

п=0

причем С^) связана с С^) следующим обр азом: <5(г„, = ^(^'га+1)С(гга, ^п)/р(гп). В Г функция С^) доопределяется следующим образом: (5^) = 0 при г € Г£, ] = 1; С^) = ^1(г*) при г € Г1е, = 0 С^) = -02(п(г)) при г € Г2е, = 0,

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

представление для ^ аналогичное (6). После его рандомизации далее реализуется алгоритм моделирования, связанный с ортом и (г). В частности, для оценки {С^(г)}^=х,у,г вектора градиента gradu(r) верно представление

«г) = 3„ = х,у,г,

где (^1) — оценка для и^1), а Wl = (г\,]\) — состояние после первого перехода в цепи Маркова; при этом коэффициенты См(г) и QMl(w1) известны уже после первого шага. Таким образом, все три компонента вектора градиента можно оценивать одновременно, т, е, на одном наборе траекторий.

4. Связь с решением краевой задачи

Теоремы в этом разделе доказываются по аналогии с теоремами из [3], где рассмотрены только краевые условия (2), Доказательства теорем учитывают тот факт, что веса в соответствующих рядах Неймана умножаются на множители Сы(кЯ(г)) < 1, а также то, что среднее число отражений в алгоритме имеет порядок 0(е-1) [1], Теорема 1. Пусть для всех г € П выполнено условие

где с0 < 6с*/п2 ~ 0.6079с*, тогда существует единственное ограниченное решение интегрального уравнения (9). Это решение представимо соответствующим рядом Неймана и совпадает с решением исходной краевой задачи (1)-(3):

и (г, 0) = и(г), =

Если, V = 0 то условие (10) можно заменить на, условие |к2 — с(г)| < с0. Теорема 2. Если, первые производные функции и(г) ограничены в П, то:

Е(е(г, 0) = и£(г) существует и |и(г) — и£(г)| = О(е), £ > 0, г € П;

Е(г) ди

Е(е(г, 1) = /£(г) существует и

3 ди

- /е(г)

= O(e), £ > 0, r е П.

Теорема 3. Пусть c0 < 0.4881c*, тогда, дисперсия, V(e ограничена для любого £> 0. Использование в [3] функции Грина оператора Лапласа, для которого к = 0 в (10),

c(r)

димое только для сходимости соответствующих рядов Неймана, В данной статье путем перехода к оператору Дк показана возможность решения уравнения (1) с большим по

-c(r)

ной -к2.

Автор выражает благодарность канд. физ.-мат, наук H.A. Симонову и чл.-корр. РАН Г.А. Михайлову за ценные замечания.

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

fl] Makarov R.N. Statistical modelling of solutions of stochastic differential equations with reflection of trajectories from the boundary // Russ. J. Numer. Anal. Math. Modelling. 2001. Vol. 16, N 3. P. 261-278.

[2] Владимиров B.C., Жарков В.В. Уравнения математической физики. \!.. 2000.

[3] Бурмистров А.В., Михайлов Г.А. Вычисление производных от решения стационарного диффузионного уравнения методом Монте-Карло // Журн. вычисл. математики и мат. физики. 2003. Т. 43, № 10. С. 1517-1529.

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

Поступила в редакцию 20 февраля 2008 г.