CC BY
52
Математическое моделирование
УДК 519.245; 519.632.4
ОЦЕНКИ МЕТОДОМ МОНТЕ-КАРЛО ИТЕРАЦИЙ ОПЕРАТОРА ГРИНА И СОБСТВЕННЫХ ЧИСЕЛ ПЕРВОЙ КРАЕВОЙ ЗАДАЧИ ДЛЯ ОПЕРАТОРА ЛАПЛАСА
А. Н. Кузнецов, И. А. Рытенкова, А. С. Сипин
Вологодский государственный педагогический университет,
160035 Вологда, ул. С. Орлова, 6.
E-mails: pm_kan@uni-Vologda. ас .ru, profyservis@mail .ru, [email protected]
Рассмотрен алгоритм вычисления степеней оператора Грина методом Монте-Карло на траекториях марковской цепи. Полученные статистические оценки используются для. нахождения первого собственного числа первой краевой задачи. Проведено численное моделирование и сделан сравнительный анализ эффективности алгоритмов.
Ключевые слова: метод Монте-Карло, собственные числа первой краевой задачи, функция Грина, оператор Грина, распределённые вычисления.
1. Постановка задачи. Рассмотрим задачу на собственные значения
А<р + \<р = 0; 'І
Иг = °; \ С1)
<р Є С2(р)пСф), VcRn. )
Множество собственных значений {Л^} этой задачи не имеет конечных предельных точек, А& > 0, Аі — простое (Аі < Аг ^ Аз ^ ...). Собственная функция <рі(х) неотрицательна, можно выбрать вещественными и орто-нормированными в С2(Т)).
Пусть С{х, у) — функция Грина для задачи (1), тогда (1) эквивалентно задаче на собственные значения для интегрального оператора (С с ядром 0(х,у):
и{х) = А / С{х,у)и{у)(1у, иєС(Т>).
Jv
Интегральный оператор (С будем назвать оператором Грина, норму в С-2(Р) будем обозначать ||- |, в С(Т>) — ||-||с-
Для определения Аі можно использовать метод Келлога [1]. Пусть
||(/?(°)||=1; ір^ = &р <р(р)=<р^||; А(р) = ||(^р_1) Ц/Ц^^ ||-
Тогда последовательность {А(р)} сходится, монотонно убывая, к Аі, а последовательность сходится к^і в С,2(Т>) и в С{Т>), причём справедливы
Андрей Николаевич Кузнецов, старший преподаватель, каф. прикладной математики. Ирина Александровна Рытенкова, магистрант, каф. прикладной математики. Александр Степанович Сипин (к.ф.-м.н., доц.), доцент, каф. прикладной математики.
оценки
Аг/
\APVT^4
2 ’
\Р(р) ~ ^1\\с ^ 1
С1
Лі V yi^cf
Cl
р = 2,3,... р = 2,3,... р = 2,3,...
где Cl = (У°), <£i), L2 = max / (G(x,y))2dy.
xeD Jd
Построению стохастических алгоритмов для определения Л посвящено много работ. Наиболее общие результаты представлены в статье [3]. В данной работе анализируются проблемы реализации лишь одного стохастического алгоритма. В определённом смысле она является продолжением работы [4].
2. Вероятностный смысл метода Келлога для оператора Грина. С оператором Лапласа связан винеровский процесс, а именно, является характеристическим оператором процесса. Пусть т = т(х) —момент первого выхода винеровского процесса £(t) из множества V, если {(0) = х, тогда справедлива
Лемма. Пусть ср^Цх) = <ро = const > 0, тогда ср^Цх) = Еж тр.
Доказательство. Пусть Л < 0 и |Л| < Ai, тогда для решения краевой задачи
Av + Xv =—(ро, v|r = 0 (2)
справедливо вероятностное представление
00 7 \\р тр+1
(p + 1)!
В последнем ряде можно поменять местами знаки математического ожидания и суммы, т. к. Еж е!Л1т/2 < +оо. Тогда
¥>о^/AV ЕХТР+1
р=о
(р + 1)!'
С другой стороны, задача (2) эквивалентна интегральному уравнению у(х) = \Gv-\-Gtpo, которое имеет решение, представимое при малых Л рядом
СЮ
у(х) = ^ АРСР+1
р=о
Сравнивая коэффициенты при Лр, имеем (р^ = &р (ро = <Ро^^г- П Следствие.
1/2
Ai = lim 2р
р^-оо
Е ХГ
,-1\2
dx
/ (Ехтр)2 dx
IV
<р \{х) = lim
р—)>С©
L Jv
(Е xTp)2dx
1/2 •
Лемма. Пусть \р){х) = 2рЩгр-^-, тогда Ит \р)(х) = Х\. Доказательство.
^Р~1\х) ч>(р- 1)(ж) Н^-Щ! У(р_ 1)(ж)
(р) Х <р^{х) <Р(р)(х) ||у?(р)|| <Р(р)(х) (р)’
Ит (р(р)(х) = (р!(х) > 0, Ит Л(р) = Аь поэтому Пт Л(р)(ж) = = Ль □
р—ос р^ос р—)> ОС 1г1\х)
3. Построение статистических оценок итераций оператора Грина. Один из вариантов метода Монте—Карло для вычисления А1 может быть основан на определении функции
Уп(х) = № 1)(ж) = Ех Тп,
т. е. на оценке моментов случайной величины т. Займёмся получением несмещённых оценок величин ьп.
3.1. Преобразование сдвига. Для построения несмещённых оценок для уп(х) используем уравнение
Аи + \и = О, «|г = 1. (3)
При |Л| < А1 решение этого уравнения имеет вид
сю
и(х, А) = Ех еЛт/2 = ^ \пуп(х), Уо(х) = 1.
п=О
Пусть с > О, РсК3. Свяжем с задачей (3) интегральное уравнение. Для шара Кц(х) радиуса К с центром в точке х (Кц(х) С V) имеем [2]:
и(х) = д[ и(у)с1ш + (1 - д) [ р(х,у)(1 + ^)и(у)с1у, (4)
Зяп )кк V С1)
где ч = ЖШ’ ш ~ равномерное распределение на сфере, р(х,у) = 4ф_у\ х
вЬ((й— \х—у\)с)
X (аЬсД-сД) — переходная плотность.
Соотношение (4) позволяет строить несмещённые оценки для и(х, А) на траекториях {Жг}°^0 случайного процесса, для которого Хо = х, х^\ распределено с вероятностью q^ равномерно на сфере и с вероятностью (1 — <^)
в шаре с плотностью р(хг,у). Нетрудно доказать, что с вероятностью 1 существует Хоо = Ит Хг И Хоо € Г.
г—)>с©
С процессом {жг}°^0 естественно связать мартингал несмещённых оценок для и(х, А): (п = (1 + Х/с2)Мп и(хп), п = 1, 2,..., где — число переходов в шар на траектории Х0,Х\,... ,хп. Известно [2], что п. н. существует N = = Пт Л^га < +оо, причём конечны моменты Ех -/V8, •§ = 1, 2,....
П—)>С©
Лемма. Пусть |Л| <р<\\/3ир< с2, тогда мартингал £п квадратично интегрируемый.
Доказательство.
^1^(1 + ^) " Нхп)\2 < (1 + “f + ^4) " \u(xn)\2 < (i + ^г) " Кж.
с с л
Зр\
^|2;
и
(х)\ = (еж еЛг/2) < (еж егр/2) < Еж ерг < Еж e3pr/2 = и(х,3р). Тогда |{2(ж,А)| ^ (l + 3p/c2)Nn u(xn, Зр), значит, |ЕжС2(ж,А)| ^ и(х,3р) <оо.П
Следствие. Случайная величина £ = (1 + А/с2)М является несмещённой оценкой для и(х, А).
Доказательство. Действительно, квадратично интегрируемый мартингал равномерно интегрируем, поэтому
и(х, А) = Ех £n{x, А) = lim Ех £n{x, А) = Ех lim £п(х, А) =
п—У СО п—У со
= Еж с • lim и{хп, А) = Еж С- □
71—У ОС
Представим £ в другой форме:
га=О
Обозначим JVN = JV(JV - 1) • • • (N - n + 1), тогда
дгМ
га=О
П!С
,2га’
(5)
га «і"1
причем ряд (5J мажорируется рядом Р і 2п; поэтому
га=О
Лг[«
га=0
следовательно, справедлива
Теорема. Случайная величина является несмещённой оценкой уп(х) с конечной дисперсией.
Доказательство. Осталось проверить конечность дисперсии. Оценим дисперсию величины Сп = ^2! • Можно считать, что М^п, иначе (п = 0.
С2 =
S п
2N • 2(N - 1) • • • 2(N - п + 1) 2п ■ 2 (п — 1) • • • 2
0Ап
ГОЛТ\\2п] 2JV 2N—2 2N -2га+2
) 2N—1 2N—3 2N-2n+l
(2п)\с4п '
2 га 2га—2 2
2га—1 2га—З 1
Функция 2^1 убывает, поэтому 2(ДГ—к)-1 ^ 2(п—к) — 11 следовательно,
2 _(2Л01^ =
и (2п)!с2'2"
Докажем, что Ехг)п < +оо. Действительно,
І Еж Є21 <
Е-(1 + ^)2М <Ф>,Зр) и еж£2 = Х>*Еж(^^),
к=О
отсюда
(2ДГ)И Ех к\с2 к < +°°- 1“1
3.2. Реализация оценки сдвига. Оценка £га не является реализуемой в силу невозможности моделирования траектории бесконечной длины. Для получения реализуемой оценки остановим мартингал ^ в момент первого попадания процесса в е-окрестность границы. Обозначим — число переходов в шар до попадания в е-окрестность границы, х£ — последнюю точку процесса, тогда
£е = (1 + Х/с2)Ме и(хе) является несмещённой оценкой и(х, А):
оо [к] оо оо га [к]
«Л> = Е Е А"М*> = Е л" Е
к=0 т=0 га=0 к=О
Отсюда
гг
к\с2к
к=О
Смещённая оценка ьп(х) получится, если в (6) отбросить слагаемые, для которых к < п. Она имеет вид
Жп'
Сп,є — | 2га' (7)
ПІС
3.3. Оценка смещения и трудоёмкость алгоритма. Оценим величину смещения уп(х) — Еж{га)е. Пусть д(е) — модуль непрерывности функции У\(х) = = (С 1)(ж), т. е. |с 1(ж)| < д(е), при Г) < е, значит, Ьк(х) = Gk~1 У\(х) ^
^ ||СЙ_1||Г5(£)- Тогда смещение ьп не превзойдёт
і с-
П— 1
• І^ІІс < 5(Ф||с||с
к=О
Как правило, д(є) = д - є, и для достижения смещения 5 необходимо взять
е^/МсНГ1)-
Для выпуклой области можно показать, что среднее число шагов процесса до попадания в е-окрестность границы имеет порядок |1пе|.
4. Модельный пример.
4.1. Оценки итераций с помощью рядов. Рассмотрим задачу (1) в кубе V = [О, I]3. Собственные числа такой задачи известны:
Ак,1,т = К2(к2 + 12 + т2), к,1,т = 1,2,3,....
В частности, А1 = Зтт2, Х2 = 67Г2, причём А2 имеет кратность 3.
Собственные функции имеют ВИД (рк,1,т(х,у, г) = Xk{x)Yl(y)Zrn(z)1 где Хк(%) = у/Ъвттгкх, и Zm строятся аналогично.
Используя разложение 1 по собственным функциям, можно получить тождество
V 83 1
у„_о (2к + 1>2<2' + 1)2<2т + 1>2 Х
г. (8)
{тг2 [(2к + I)2 + (21 + I)2 + (2т + I)2]Г'
Для определения суммы ряда (8) можно вычислять его частичные суммы
П
Бп= Е ••••
к,1,т=0
Воспользовавшись неравенством а + Ь + с ^ 3 \J~abc (а ^ О, Ь ^ 0, с ^ 0), получим
1 1
[(2к + I)2 + (21 + I)2 + (2т + 1)2]р [3^(ЖTIjЧ2ІTIj2(2mTIj2F
что влечёт оценку остатка Ям = (1, Ср 1) — вм- А именно, Ям ^ Ям, где Ям построено для ряда
^ 83 1
Зртг6+2р (2к + 1)2+2р/3(2/ + 1)2+2р/3(2т + 1)2+2р/3 ’
к, /, тп—0
который распадается в произведение трех рядов.
о° 3
Пусть Ар = £ {2к+1)^Р/3 ’ ВР = зрД+^р> тогда
Ям = ЗВр
^*£-1 р»+1)2+«31
Нетрудно получить оценку для гм и Ар:
/ 1 1 / 1 ™ ^ 2 (1 + 2р/3)(2Ы + 1)1+2р/з ’ р ^ 1 + Г° “ 1 + 2(1 + 2р/3) ’
Ям ^ Длг ^ ЗВР ( 1 +
6 + 4р/ 6 + 4р (27У + 1)!+2р/3 ’
Таблица 1
Оценка итераций оператора Грина
р Е N (1, Єр1)
1 1(Г4 23 0,02017
2 1(Г7 25 6,243 6- 10~4
4 ю-12 18 6,942 341 • 10-7
6 1(Г17 16 7,905 050 4- Ю-10
8 10-22 15 9,015 860 062 6- 10-13
10 10-22 3 1,028 397 598- Ю^15
Таблица 2
Оценки А(р)
п А(о) А(1) А(2) А(з) А(4)
Точные значения
40,021 29,989 29,635 29,611 29,609
Расчёт на МРРкластере
ю5 39,576 29,522 30,233 32,835 36,919
ю6 40,048 30,222 30,797 33,055 36,796
ю7 40,020 30,008 29,820 30,381 31,703
ю8 40,011 29,976 29,657 29,843 30,458
ю9 40,021 29,996 29,673 29,768 30,092
Расчёт с использованием видеокарты
105 40,500 30,157 29,014 28,494 29,761
106 40,115 29,966 29,726 30,452 32,402
107 40,057 30,010 29,779 30,296 31,580
108 40,025 29,996 29,663 29,778 30,213
109 40,021 29,989 29,635 29,638 29,763
По заданной точности Е вычислялись N и итерации оператора. Результаты сведены в табл. 1.
Полученные на основе этих значений Л(р) = ^/(^’^ттху приведены в
табл. 2. Истинное собственное число А1 = 37Г2 = 29,608 813.
4.2. Реализация оценки (7). Наибольшую трудность при реализации оценки (7) представляет моделирование плотности вероятности перехода р(х, у) = = д5к(у) + (1 — д)р(х,у), где — равномерное распределение на сфере с
/ \ С2 ёЫсК— СГ)
центром в точке х, а р(х,у) = ш 8ьсд_сд •
Моделирование р(х, у) осуществлялось методом отбора.
Пусть а = сК и р(г) — плотность на [0, а], г) — случайная величина с
плотностью р(г).
Лемма. Пусть Т 8^а'а ^ < р(г), тогда метод отбора моделирует смесь
а . г$Ъ(а-г) . .
лн<(в-г) + -к-^- (9)
Доказательство. Определим случайную величину ( условиями:
^ _ I г], если ар(г)) < І а, в противном
г/з\ї(а—7]) . эЬ а ’
случае.
Далее
Г“ г вЬГа — г) , а
= 1-1 -------1------йг = -—.
впа впа
Jo
Пусть ж < а, тогда РС(х) =р((<х) = р[г]< х, ар (г/) < Г?8Ь^~Г?)) =
ГХ
= Л, вЪа-р(г) Уо йЬа
В качестве р(г) можно выбрать р(г) = гсЬ(а — г)/(сЬа — 1), тогда, очевидно,
г8Ь(а — г) гсЬ(а — г)
вЬ а сЬ а — 1
Лемма. Пусть {і и {2 независимые случайные величины, имеющие экспоненциальное распределение с плотностью е~х, х > 0, г\ = {і — а[{і/а], -22 = £2 — ^[Сг/о]) тогда случайная величина
( г\ + 22, если 21 + 22 < О/
^ \2а — (21+22), если 2і + 22 ^ а;
имеет плотность р(г).
Доказательство. Ясно, что 0 ^ 21 < а. Далее
Р(гі <х) = ^ Р(/са ^ {і < ж + ка) = ^ е ка — є
—ка—х
1-е-"’ &=0 &=0
Р(2 < х) = Р(2і + 22 < ж) + Р(2а — (21 + 22) < Ж, 21 + 22 ^ а) =
=--------------? (1 - (ж + 1)е-ж + (ж - 1)еж-2“ + е-2“) ;
(1 — є~а) 7
рг(ж) = ------------2 ((ж + 1>~Х - е~Х + еЖ_2“ + (ж - 1)еЖ_2“) =
(1 — е~а)
1 / г,. ™ о„ч 2е_“жсЬ(а —ж) жсЬ(ж —а)
- (же“ж + хех~2а) = ---- ---- -----/ = —^---------------±. □
' 'Л о/о—а і /о—2сь /".І-і /'і 1
(1 _ е_а)2 1 — 2е_“ + е_2“ сЬ а — 1
Моделирование переходной плотности р(х, у) теперь не вызывает затруднений.
Лемма. Пусть случайная величина г] имеет плотность (9) (а = cR), ш —равномерно распределено на единичной сфере, тогда случайный вектор у = х + г]ш/с имеет плотность распределения р(х, у).
4.3. Практическая реализация. От значения параметра с зависят средняя длина траектории и количество переходов в шар (табл. 3). Экспериментальным путём было выбрано значение с = 15; использованное значение е = = 10“ . Описанный алгоритм был реализован на языке C++.
Как показывают результаты вычислений (табл. 2, 4), для получения приемлемой точности требуется моделирование порядка 109 траекторий, что требует значительных вычислительных ресурсов. При моделировании на одном ядре процессора Intel® Core™2 Quad CPU Ц8400@2,66ГГц необходимое процессорное время превысило бы 27 часов.
Для получения значений, указанных в первой части таблиц, использовались распределённые вычисления с использованием технологии MPI (http://www.mpi-forum.org) В реализации MPICH2 (http://www.mcs.anl.gov/research/projects/mpich2/). С применением 13 ПК на базе Intel® Core™2 Quad CPU Q8400S2,66 ГГц, были запущены 50 вычислительных и один распределяющий процесс, что позволило сократить время расчёта до 36 минут. Для обеспечения независимости псевдослучайных чисел использовался параллельный 128-битный генератор псевдослучайных чисел [5]. Значения во второй части таблицы получены при вычислении на видеокарте с графическим процессором NVIDIA® GeForce® 9600 GT с применением технологии CUDA. Время расчёта составило 10 минут. Для компиляции использовался пакет CUDA Toolkit V4.0 RC2 (http://developer.nvidia.com/cuda-toolkit-40), псевдослучайные числа генерировались с помощью входящей в него библиотеки CURAND. В зависимости от возможностей видеокарты вычисления с плавающей точкой могут вестись с одинарной или двойной точностью (использованная видеокарта обладает «вычислительными возможностями» 1.1 и поддерживает лишь одинарную точность). В случае, если работа с двойной точностью не поддерживается, при моделировании большого количества траекторий нужно использовать те или иные способы уменьшения вычислительной погрешности при суммировании. Эксперименты показали, что для данной задачи применимы следующие методы:
- суммирование с компенсацией по формуле Кахана [6];
- вычисление на графическом процессоре за один раз не очень большого количества траекторий (10 -105) с последующим суммированием и осреднением на ПК с двойной точностью (это также помогает производит вычисления в случае, если нет возможности отключить монитор от видеокарты, чтобы избежать двухсекундного ограничения на время выполнения ядра).
- использование достаточного количества блоков (thread block) в сетке (grid) с последующим суммированием и осреднением результатов, полученных различными потоками сетки, не на видеокарте, а на ПК с двойной точностью.
Таблица 3 Зависимость длины траектории I и среднего количества переходов в шар п от с
с 1 п
5 121 0
10 123 2
15 125 4
20 128 8
25 132 12
Таблица 4
Оценки др = (1, (Ср 1)
п 91 92 94. 96 98 9 ю
Точные значения
0,020 170 6,243 6- 10-4 6,942 3- 10-7 7,905 0- Ю-10 9,015 8- 10~13 1,028 3- ю-15
Расчёт на МР1-кластере
ю5 0,020 344 6,384 5 ю-4 7,325 7- 10-7 8,014 6- Ю-10 7,433 7- ю-13 5,453 8 • 10-1б
ю6 0,020 177 6,2351 ю-4 6,826 4- 10-7 7,1972- Ю-10 6,586 9 • 10-13 4,865 0 • 10-1б
ю7 0,020 168 6,243 7 10-4 6,933 9- 10-7 7,797 5- Ю-10 8,447 8- 10-13 8,405 2- 10-1б
ю8 0,020 172 6,246 6 ю-4 6,9519- 10-7 7,9041 • Ю-10 8,874 9- 10-13 9,566 5 • 10-1б
ю9 0,020 169 6,243 5 ю-4 6,938 9- 10-7 7,880 5- Ю-10 8,892 9- 10-13 9,820 4- Ю-16
Расчёт с использованием видеокарты
105 0,019 944 6,096 8 10-4 6,703 9 ю-7 7,963 4 ю-10 9,808 5 ю-13 1,1074 ю-15
106 0,020101 6,2144 ю-4 6,920 4 ю-7 7,8316 ю-10 8,445 3 ю-13 8,043 8 Ю-16
107 0,020 142 6,232 3 ю-4 6,920 0 ю-7 7,803 3 ю-10 8,5019 10-13 8,5251 10-1б
108 0,020 167 6,242 3 10-4 6,937 9 ю-7 7,885 1 ю-10 8,892 5 10-13 9,7417 10-1б
109 0,020 170 6,243 5 ю-4 6,942 3 ю-7 7,904 7 ю-10 8,998 7 10-13 1,0158 ю-15
СО
Оценки методом Монте—Карло итераций оператора Грина и собственных чисел .
Работа выполнена при поддержке РФФИ (проект № 11-01-00769-а).
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Владимиров В. С. Уравнения математической физики. М.: Наука, 1967. 436 с. [Vladimirov V.S. The equations of mathematical physics. Moscow: Nauka, 1967. 436 pp.]
2. Ермаков С. М., Некруткин В. В., Сипин А. С. Случайные процессы для решения классических уравнений математической физики. М.: Наука, 1984. 208 с. [Ermakov S. М., Nekrutkin V. V., Sipin A. S. Random processes for solving classical equations of mathematical physics. Moscow: Nauka, 1984. 206 pp.]
3. Михайлов Г. А., Макаров P. H. Параметрическое дифференцирование и оценки собственных чисел методом Монте-Карло// Сиб. матем. журн., 1998. Т. 39, №4. С. 931-941; англ. пер.'.Mikhaylov С. A., Makarov R. N. Parametric differentiation and estimation of eigenvalues by the Monte Carlo method // Siberian Mathematical Journal, 1998. Vol. 39, no. 4. Pp. 806-815.
4. Кузнецов A.H., Сипин А. С. Статистические оценки для степеней оператора Грина// Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2009. №2(19). С. 114-123. [Kuznetsov A. N., Sipin A. S. Statistical estimates for the degrees of Green operator // Vestn. Samar. Cos. Tekhn. Univ. Ser. Fiz.-Mat. Nauki, 2009. no. 2(19). Pp. 114-123].
5. Михайлов Г. А., Марченко М. А. Параллельная реализация статистического моделирования и генераторов случайных чисел: Препринт РАН; Сибирское отделение; Институт вычислительной математики и математической геофизики; № 1154. Новосибирск, 2001. 20 с. [Mikhaylov С. A., Marchenko М. A. Parallel realization of statistical simulation and random number generators: Preprint of RAS; Siberian Branch; Institute of Computational Mathematics and Mathematical Geophysics; No. 1154. Novosibirsk, 2001. 20 pp.]
6. Kahan W. Pracniques: further remarks on reducing truncation errors // С ACM, 1965. Vol. 8, no. 1. Pp. 40.
Поступила в редакцию ll/V/2011; в окончательном варианте — 15/XI/2011.
MSC: 65C05; 35J08
MONTE-CARLO ESTIMATIONS FOR POWERS OF GREEN OPERATOR AND THE FIRST EIGENVALUE FOR DIRICHLET BOUNDARY VALUE PROBLEM
A. N. Kuznetsov, I. A. Rytenkova, A. S. Sipin
Vologda State Pedagogical University,
6, S. Orlova St., Vologda, 160035, Russia.
E-mails: [email protected], [email protected], cacSuni-vologda.ac.ru
In this paper, we examine the algorithm for computing the powers of a Green operator and the first eigenvalue for the Dirichlet boundary value problem using Monte-Carlo method. The efficiency of numerical realization of these algorithms is also discussed.
Keywords: Monte-Carlo method, eigenvalues of the Dirichlet boundary value problem, Green function, Green operator, distributed computing .
Original article submitted ll/V/2011; revision submitted 15/XI/2011.
Andrey N. Kuznetsov, Senior Teacher, Dept, of Applied Mathematics. Irina A. R'ytenkova, Master Student, Dept, of Applied Mathematics. Aleksandr S. Sipin (Ph.D. (Phys. & Math.)), Associate Professor, Dept, of Applied Mathematics.
