close

Вход

Забыли?

вход по аккаунту

?

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

код для вставкиСкачать
Вычислительные технологии
Том 13, Специальный выпуск 4, 2008
Модиикация случайных блужданий для оценки
градиента решения эллиптических краевых задач
?
А. В. Бурмистров
Институт вычислительной математики
и математической геоизики СО АН, Новосибирск, оссия
Новосибирский государственный университет, оссия
e-mail: burmosmf.ss.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. Диеренциальная задача и определения
ассмотрим смешанную краевую задачу:
?u(r) ? c(r)u(r) + v(r), ?u(r) = ?g(r), r ? ? ? R3 ;
u(y) = ?1 (y), y ? ?1 ;
(?u(y), ?(y)) + ?(y)u(y) = ?2 (y), y ? ?2 ,
(1)
(2)
(3)
S
внутри области ? с границей ? = ?1 ?2 , которая предполагается кусочно-гладкой.
Будем считать, что коэициенты уравнения (1), граничные ункции и область ?
таковы, что существует единственное достаточно гладкое решение исходной задачи.
Область может быть неограниченной, в этом случае потребуем от решения u(r) следующего поведения на бесконечности: u(r) = O(r ?1 ) при |r| ? ?. Нас будет интересовать
оценка решения u(r) и градиента ?u(r) в некоторой точке r ? ?. Введем следующие
обозначения:
?? = ? ? ?2 стационарный диузионный оператор. Перепишем уравнение (1),
оставив только этот оператор в левой части:
?? u(r) ? ?u(r) ? ?2 u(r) = ? ?2 ? c(r) u(r) ? v(r), ?u(r) ? g(r);
(4)
R = R(r) = min{dist(r, ?), Rmax }; c0 R2 (r)
c0 константа такая, что p(r) = 1 ?
? (0, 1);
6
?
абота выполнена при инансовой поддержке оссийского онда ундаментальных исследований
(грант ќ 08-01-00334).
Институт вычислительных технологий Сибирского отделения оссийской академии наук, 2008.
20
Модиикация случайных блужданий для оценки градиента решения...
21
B(r0 , R) = {r ? R3 : |r0 ? r| ? R(r0 )} шар с центром в r0 , целиком лежащий в ?;
S(r0 , R) = ?B(r0 , R) = {r ? R3 : |r0 ? r| = R(r0 )} соответствующая сера;
?? = {r ? ? : ?r ? ? ?, |r ? r ? | ? ?} ?-окрестность границы ?;
? единичный вектор направления v(r), т. е. v(r) = |v(r)|?(r); a? = cos(r ?d
r ? , ?);
2
?c? = ? 2 /Rmax
первое собственное число оператора Лапласа в шаре радиуса Rmax ;
?
r проекция точки r ? ?1? на ?1 : |r ? r ? | = R(r). Условие (2) аппроксимируем в
?1? следующим образом: u(r) = ?1 (r ? ) + O(?) ? ?1 (r ? );
d? (r) расстояние от точки r до границы ?2 вдоль векторного поля ? ;
?(r) проекция точки r ? ?2? на ?2 вдоль поля ? . Условие (3) аппроксимируем в
?2? следующим образом [1?:
u(r) =
[1 + ?d? (r)]u(r ? ??)
??2 (?(r))
+
+ O(?2 ) ? p?(r)u(r ? ??) + ?e2 (?(r));
1 + ?(d? (r) + ?)
1 + ?(d? (r) + ?)
Pn (cos ?)r полиномы Лежандра [2?;
?
?n (x) =
I 1 (x), где In+ 1 (x) модиицированные ункции Бесселя [2?;
2
2x n+ 2
sinh(x)
1
sinh(x)
, ?1 (x) =
cosh(x) ?
?0 (x) =
;
x
x
x
Gr?0 (r, r ? ) нецентральная ункция рина для оператора ?? в шаре B(r0 , R), которая является решением следующей задачи Дирихле:
?? Gr?0 (r, r ? ) = ?(r ? r ? ), Gr?0 (r, r ? )r?S(r0 ,R) = 0,
где ?(·) дельта-ункция Дирака. Явные выражения для ункции Gr?0 и ее производных, используемые для решения поставленной задачи, выписаны далее.
2. Сведение к интегральному уравнению
Формально интегральное уравнение для ункции u(r ? ) из (4) и ее пространственной
?u ?
производной
(r ) можно получить, используя теорему о среднем. В нашем случае вос??
пользуемся нецентральной ункцией рина Gr?0 (r, r ? ) для оператора ?? в шаре B(r0 , R).
В результате по аналогии с [3? для r ? , r0 ? ? \ ?? будут иметь место представления:
?
u1 (r ) = ?
Z
S(r0 ,R)
?u1 ?
(r ) = ?
??
Z
S(r0 ,R)
Z
?Gr?0
?
(r, r )u1 (r)dSr +
Gr?0 (r, r ? )F (u, ?, r)dr;
?nr
B(r0 ,R)
?
Z
?Gr?0
? ?Gr0
?
(r, r )u1 (r)dSr +
(r, r ? )F (u, ?, r)dr.
?? ?nr
??
(5)
(6)
B(r0 ,R)
?u1
2
Здесь F (u1 , ?, r) ? |v(r)|
(r) + (? ? c(r))u1 (r) + g(r) ; nr внешняя нормаль к
??
сере S(r0 , R). Для r ? ? ?? полагаем u1 ? u. Далее будут сормулированы условия
существования и единственности системы уравнений (5) и (6) и совпадения их с решением исходной диеренциальной задачи и его пространственной производной. Теперь
22
А. В. Бурмистров
выпишем ункции, использованные в представлениях (5) и (6), а также найдем плотности вероятностей, пропорциональные им. Нецентральную ункцию рина Gr?0 (r, r ? )
будем искать в следующем виде:
1 sinh {?(R ? |r ? r ? |)}
?
?
? G(r, r ? ),
Gr0 (r, r ) =
(7)
4?
sinh {?R} |r ? r ? |
где ункция G(r, r ? ) является решением следующей краевой задачи:
1 sinh {?(R ? |r ? r ? |)}
?
?
, |r ? r0 | = R.
?? G(r, r ) = 0, |r ? r0 | < R, G(r, r ) =
4?
sinh {?R} |r ? r ? |
(8)
ешение краевой задачи (8) можно записать в серических координатах [2?:
?
G(r, r ) =
?
X
n=0
?n (r ? )Pn (?)
?n (?|r ? r0 |)
,
?n (?R)
где ? = cos ? косинус угла между векторами (r ? r0 ) и (r ? ? r0 ), а коэициенты ?n (r ? )
определяются по ормуле
n
o?
?
p
Z1
2 + |r ? ? r |2 ? 2R|r ? ? r |?)
sinh
?(R
?
R
0
0
2n + 1
1 ?
? Pn (?)d?.
p
?n (r ? ) =
2
4?
sinh {?R} R2 + |r ? ? r0 |2 ? 2R|r ? ? r0 |?
?1
В дальнейшем мы будем использовать в основном центральную ункцию рина для
оператора ?? , т. е. значения ункции u1 (5) и ее градиента (6) в точке r ? будем представлять в виде суммы интегралов по сере и шарам с центром в той же точке (r0 = r ?).
Отметим, однако, что для построения производных была использована именно нецентральная ункция рина (7), кроме этого, нецентральные ункции можно использовать при решении внешних краевых задач. При этом для ограниченности на бесконечности в выражениях для ункций ?n (x) модиицированные ункции Бесселя
In+ 1 (x) нужно заменить на ункции Kn+ 1 (x) [2?. После диеренцирования, обозна2
2
чая ? = |r ? r ? |, ??(?) = ?? cosh{?(R ? ?)} + sinh{?(R ? ?)}, а также учитывая, что
полиномы Лежандра ортогональны и P1 (cos ?) = cos ? [2?, получим
sinh {?(R ? ?)}
?Gr??
a?
??(?)
? ?1 (??)
?
?
?
Gr? (r, r ) =
,
(r, r ) =
?
,
4?? sinh {?R}
??
4? sinh{?R} ?2
R ?1 (?R)
?Gr??
1
1
? ?Gr??
a?
?
?
(r, r ) = ?
,
(r, r ? ) = ?
.
2
2
?nr
4?R ?0 (?R)
?? ?nr
4?R ?1 (?R)
Заметим, что при ? ? 0 оператор ?? стремится к оператору Лапласа ? и вышеприведенные ункции стремятся к соответствующим центральным шаровым ункциям для
оператора Лапласа [3?:
1
1
1
?Gr0?
a?
1
|r ? r ? |
?
0
?
Gr? (r, r ) =
?
,
(r, r ) =
?
,
4? |r ? r ? | R
??
4? |r ? r ? |2
R3
?Gr0?
1
? ?Gr0?
a? 3
?
(r, r ) = ?
,
=
?
.
?nr
4?R2
?? ?nr
4?R2 R
Модиикация случайных блужданий для оценки градиента решения...
23
Полученные ункции пропорциональны плотности распределения на сере FS и плотностям распределения в шаре F0 и F1 :
Gr? (r, r ?)
?
R2
=
C01 (?R)F0 (r, r ? ),
6
?Gr?
(r, r ?) = C00 (?R)FS (r, r ? ),
?nr?
?Gr?
3R
(r, r ? ) = a?
C11 (?R)F1 (r, r ? ),
?? 4
?
?Gr?
3
?
(r, r ?) = a? C10 (?R)FS (r, r ? ).
?? ?nr?
R
Здесь ункции Ckl (·), k, l ? {0, 1} монотонно убывают от 1 до 0:
x
,
3?1 (x)
4
t0 (x) ? 1
tanh(x/2)
C11 (x) =
t0 (x) ?
, t0 (x) =
.
3
x?1 (x)
x/2
C00 (x) = ?0 (x)?1 ,
C01 (x) =
6 ?0 (x) ? 1
,
x2 ?0 (x)
C10 (x) =
Функция FS является плотностью равномерного распределения на сере:
FS (?, ?)dS =
sin(?)d?d?
,
4?
здесь ? ? (0, ?), ? ? (0, 2?) и ? = |r ? r ? | ? (0, R) координаты локальной (с центром
в точке r ? ) серической системы координат. Плотности F0 и F1 в этой системе координат акторизуются: Fj (r, r ? )dr = FS (?, ?)d?d?Fj?(?)d?, j = 0, 1, а части Fj? (?) имеют
следующий вид:
?1
?1
6C01
(?R)
4C11
(?R)
??2 ?1 (??)
?
?
? sinh {?(R ? ?)} , F1 (?) = 2
??(?) ?
.
F0 (?) = 2
R sinh {?R}
R sinh{?R}
R ?1 (?R)
Моделирование случайных величин с такими плотностями можно осуществлять методом исключения [4?.
Для построения статистических алгоритмов решения системы уравнений (5) и (6)
объединим их далее в одно интегроалгебраическое уравнение с помощью специального
расширения множества азовых координат дискретной переменной j , принимающей
только два значения: j = 0 или j = 1. Обозначим точку нового азового пространства
через w = (r, j) ? R3 Ч {0, 1} и введем также следующее обозначение:
?
?u1 (r),
j = 0;
U(w) ? U(r, j) = R(r) ?u1
?
(r), j = 1.
3 ??
1 ? c0 R2 (r)
для рандомизации полученного уравнения:
6
?
?
Z
G(r ? , j ? ) ?
?
U(r ? , j ? ) = p(r ? ) ?
FS (r, r ? )U(r, 0)Qj0 (r, r ? )dSr +
?+
p(r ? )
S(r ? ,R)
Z
+(1 ? p(r ? ))
Fj (r, r ? ) pc1 Qcj? 1 (r, r ? )U(r, 1) + pc0 Qcj? 0 (r, r ? )U(r, 0) dr,
(9)
Используем величину p(r) =
B(r ? ,R)
24
А. В. Бурмистров
где pc0 = pc0 (r, r ? ) и pc1 = pc1 (r, r ?) тоже являются некоторыми вероятностями
pc0 + pc1 = 1, которые можно взять равными pc1 = pc0 = 0.5, однако более эективно
брать их пропорциональными ункциям |v| и |?2 ? c|:
pc1 =
|?2 ? c(r)|
3|v(r)|
3|v(r)|
, pc0 =
, где pcv =
+ |?2 ? c(r)|.
R(r)pcv
pcv
R(r)
Функция G(w) представима в виде интегралов [Ckl ? Ckl (?R(r ? )); k, l ? {0, 1}?:
Z
Z
R2
R2
?
?
?
G(r , 0) = C01
F0 (r, r )g(r)dr, G(r , 1) = C11
F1 (r, r ? )g(r)a? dr,
6
4
B(r ? ,R)
B(r ? ,R)
которые можно оценивать по одному случайному узлу [4?. Веса выражаются так:
3|v(r)|C01
9a? |v(r)|C11
a? C10
, Qc01 (r, r ? ) = c
, Qc11 (r, r ? ) =
,
?
p(r )
p1 c0 R(r)
2pc1 c0 R(r)
C00
(?2 ? c(r))C01
3a? (?2 ? c(r))C11
c
?
c
?
Q00 (r, r ? ) =
,
Q
(r,
r
)
=
,
Q
(r,
r
)
=
.
00
10
p(r ? )
pc0 c0
2pc0 c0
Q10 (r, r ? ) =
3. Алгоритм решения интегрального уравнения
В соответствии с (9) и с учетом введенных обозначений статистический алгоритм сопряженных блужданий по серам и в шарах для оценки величины U(w0 ) строится
следующим образом.
[w0 = (r0 , j0 ), Q0 = 1; траекторию моделируем до обрыва на шаге со случайным
номером N ?.
? Если rn 6? ?? , то ледующая точка rn+1 выбирается на сере S(rn ) с вероятностью
p(rn ) (случай 1) или в шаре B(rn ) с дополнительной вероятностью 1 ? p(rn ) (случай 2).
Для случая 1 точка rn+1 выбирается равномерно на сере S(rn ), т. е. по плотности
FS (rn , r ?). В счетчик добавляется величина Qn p?1 (rn )G(wn ). Величина jn+1 полагается
равной 0. Вес Qn умножается на величину Qjn 0 : Qn+1 = Qn Qjn 0 .
Для случая 2 точка rn+1 выбирается в шаре B(rn ) по плотности Fjn (rn , r ? ). Величина jn+1 полагается равной 1 с вероятностью pc1 (rn , rn+1 ) и равной 0 с вероятностью
pc0 (rn , rn+1 ). Вес Qn умножается на величину Qcjn jn+1 : Qn+1 = Qn Qcjn jn+1 .
Продолжаем алгоритм от места, обозначенного ?, для точки rn+1 .
Если rn ? ?1? , jn = 0, то в счетчик добавляется величина Qn ?1 (rn? ) и траектория
обрывается. Если rn ? ?2? , jn = 0, то в счетчик добавляется величина Qn ?e2 (?(rn )).
С вероятностью p?(rn ) следующая точка rn+1 = rn ? ?n ?n отражается в ? \ ?? . Величина
jn+1 полагается равной 0. Вес не изменяется. Продолжаем алгоритм от места, обозначенного ?, для точки rn+1 . С вероятностью 1 ? p?(rn ) траектория обрывается.
Если rn ? ?? , jn = 1, то траектория обрывается.
В результате для величины U(w0 ) получаем реализуемую оценку:
?? (w0 ) =
N
X
Qn G?(wn ),
n=0
причем G?(w) связана с G(w) следующим образом: G?(rn , jn ) = ?(jn+1 )G(rn , jn )/p(rn ). В
?? ункция G?(w) доопределяется следующим образом: G?(w) = 0 при r ? ?? , j = 1;
G?(w) = ?1 (r ? ) при r ? ?1? , j = 0; G?(w) = ?e2 (?(r)) при r ? ?2? , j = 0.
Модиикация случайных блужданий для оценки градиента решения...
25
Замечание. Для оценки производной от решения u по любому направлению µ (а не
только по направлению скорости ? ) на первом шаге алгоритма следует использовать
?u1
представление для
, аналогичное (6). После его рандомизации далее реализуется ал?µ
горитм моделирования, связанный с ортом ?(r). В частности, для оценки {?µ (r)}µ=x,y,z
вектора градиента gradu(r) верно представление
?µ (r) = 3
Gµ (r) + Qµ (w1 )?(w1)
, µ = x, y, z,
R(r)
где ?(w1 ) оценка для U(w1 ), а w1 ? (r1 , j1 ) состояние после первого перехода в цепи
Маркова; при этом коэициенты Gµ (r) и Qµ (w1 ) известны уже после первого шага.
Таким образом, все три компонента вектора градиента можно оценивать одновременно,
т. е. на одном наборе траекторий.
4. Связь с решением краевой задачи
Теоремы в этом разделе доказываются по аналогии с теоремами из [3?, где рассмотрены
только краевые условия (2). Доказательства теорем учитывают тот акт, что веса в
соответствующих рядах Неймана умножаются на множители Ckl (?R(r)) ? 1, а также
то, что среднее число отражений в алгоритме имеет порядок O(??1 ) [1?.
Теорема 1. Пусть для всех r ? ? выполнено условие
|?2 ? c(r)| +
c0 < 6c? /? 2 ? 0.6079c?, тогда
интегрального уравнения (9). Это
где
3|v(r)|
2c0
?
,
R(r)
3
существует единственное ограниченное решение
решение представимо соответствующим рядом
Неймана и совпадает с решением исходной краевой задачи
U(r, 0) = u(r), U(r, 1) =
Если
v ? 0,
то условие
(10)
(10)
(1)(3):
R(r) ?u
(r).
3 ??
|?2 ? c(r)| ? c0 .
u(r) ограничены в ?, то:
можно заменить на условие
Теорема 2. Если первые производные ункции
E??(r, 0) = u? (r)
существует и
E?? (r, 1) = f? (r)
существует и
|u(r) ? u? (r)| = O(?), ? > 0, r ? ?;
R(r) ?u
= O(?), ? > 0, r ? ?.
?
f
(r)
?
3 ??
c0 < 0.4881c? , тогда дисперсия V?? ограничена для любого ? > 0.
Использование в [3? ункции рина оператора Лапласа, для которого ? = 0 в (10),
накладывало неизическое требование малости на модуль коэициента c(r), необходимое только для сходимости соответствующих рядов Неймана. В данной статье путем
перехода к оператору ?? показана возможность решения уравнения (1) с б
ольшим по
модулю отрицательным коэициентом ?c(r), который близок к некоторой постоянной ??2 .
Автор выражает благодарность канд. из.-мат. наук Н.А. Симонову и чл.-корр. АН
.А. Михайлову за ценные замечания.
Теорема 3. Пусть
26
А. В. Бурмистров
Список литературы
[1?
[2?
[3?
[4?
Statistial modelling of solutions of stohasti dierential equations with
reetion of trajetories from the boundary // Russ. J. Numer. Anal. Math. Modelling. 2001.
Vol. 16, N 3. P. 261278.
Makarov R.N.
Владимиров В.С., Жарков В.В.
Уравнения математической изики. М., 2000.
Бурмистров А.В., Михайлов .А. Вычисление производных от решения стационарного
диузионного уравнения методом Монте-Карло // Журн. вычисл. математики и мат.
изики. 2003. Т. 43, ќ 10. С. 15171529.
Ермаков С.М., Михайлов .А.
Статистическое моделирование. М.: Наука, 1982.
Поступила в редакцию 20 евраля 2008 г.
Документ
Категория
Без категории
Просмотров
3
Размер файла
216 Кб
Теги
решение, оценки, случайных, эллиптическая, градиент, блужданий, модификация, задачи, краевых
1/--страниц
Пожаловаться на содержимое документа