УДК 519.642.4
А. А. Никитин, М. В. Николаев2
ИССЛЕДОВАНИЕ ИНТЕГРАЛЬНОГО УРАВНЕНИЯ РАВНОВЕСИЯ С ЯДРАМИ-КУРТОЗИАНАМИ В ПРОСТРАНСТВАХ РАЗЛИЧНЫХ РАЗМЕРНОСТЕЙ*
В настоящей работе продолжается изучение интегральных уравнений, возникающих в модели стационарных биологических сообществ, с ядрами, имеющими переменные коэффициенты эксцесса, ядрами-куртозианами. Рассматривается зависимость первого и второго пространственных моментов от размерности окружающей среды. Разрабатывается алгоритм быстрого вычисления многомерной нелинейной свертки. Доказывается существование радиально-симметричного решения.
Ключевые слова: математическое моделирование, интегральные уравнения, численные методы, математическая биология.
1. Введение. Известно, что использование пространственных моделей, учитывающих поведение каждого индивида в отдельности (individual-based models), приводит к более реалистичным результатам, чем использование однородных моделей, работающих с "хорошо перемешанными" (well-mixed) популяциями [1]. Несмотря на это, пространственно-неоднородные модели долгое время были весьма сложны для анализа. Однако со временем стали появляться доступные вычислительные технологии, а также более простые модели, учитывающие пространственную структуру биологических сообществ. Об одной из таких моделей и будет идти речь в этой статье.
2. Пространственная модель стационарных сообществ. Модель, рассматриваемая в нашей работе, была предложена У. Дикманом и Р. Лоу в [1, 2] и является обобщением модели [3]. Рассмотрение этой модели привело к изучению системы интегро-дифференциальных уравнений, численное решение которой и является главным предметом изучения настоящей работы. Прежде чем перейти к описанию результатов проведенного исследования, кратко опишем вышеназванную модель. Более подробно о ней можно прочитать в работах [1, 2, 4].
Рассматривается популяция стационарных биологических сообществ, обитающих в некоторой области А С 1". Все индивиды являются материальными точками и ничем не отличаются друг от друга, кроме положения в пространстве. Время в модели непрерывно, и в каждый момент с любым из индивидов могут произойти два события: рождение и смерть. Не вдаваясь в подробности, будем считать, что за эти события отвечают заданные моделью функции и параметры: т(х) — ядро рассеивания; ш(х) — ядро конкуренции; b — темп рождаемости; d — темп смертности; d! — сила конкуренции.
В настоящей статье мы продолжаем начатое в [5] (где рассматривались линейные интегральные уравнения с ядрами, представляющими собой распределение Стьюдента) исследование интегральных уравнений, возникающих в вышеназванной модели в случае, когда ядра имеют переменный коэффициент эксцесса. Эта характеристика отражает меру "пиковости" распределения случайной величины. Здесь будут приведены результаты, касающиеся изучения уравнений с ядрами рассеивания и конкуренции, имеющими следующий вид:
1 Факультет ВМК МГУ, доц., к.ф.-м.н., НИУ ВШЭ, доц., e-mail: nikitinQcs.msu.ru
2 Факультет ВМК МГУ, студ., e-mail: nikolaev.mihailQinbox.ru
* Публикация подготовлена в результате проведения исследования (№ проекта 18-05-0011) в рамках Программы «Научный фонд НИУ ВШЭ» в 2018-2019 гг. и в рамках государственной поддержки ведущих университетов Российской Федерации "5-100".
где константы ст и сш находятся из условий: f т(х) dx = 1, f ш(х) dx = 1. Данные функции при-
Кп Кп
нято называть куртозианами (kurtosian) [6]. В случае, когда параметры s0 и si равны, функция называется мезокуртозианой (mesokurtic function), если s0 < si, то — плэтикуртозианой (platykurtic function), если s0 > si, то — лептокуртозианой (leptokurtic function). В случае мезокуртозианы, которой является, например, плотность распределения Гаусса с нулевым математическим ожиданием и единичной дисперсией, коэффициент эксцесса будет равен нулю. В случае плэтикуртозианы и лептокуртозианы коэффициент эксцесса будет отрицательным и положительным соответственно.
Рассматриваемая модель У. Дикмана и Р. Лоу является стохастической, и одноразовые симуляции описываемого моделью процесса не позволяют сделать какие-либо выводы о конкретном поведении индивидуумов. Поэтому было предложено перейти к анализу математических ожиданий величин, которые создатели модели назвали пространственными моментами.
Для того чтобы ввести данные величины, зафиксируем некоторый момент времени и рассмотрим функцию р(х), называемую пространственным паттерном, и равную р(х) = ^ 8(х — х'), где
х'ех
X — множество позиций всех индивидов, а 8(х) — многомерная дельта-функция Дирака. Далее вводятся плотность популяции N(p) (нормированное на площадь области число индивидов в этой области):
N(p) = щ J Р(х) dx,
А
плотность пар популяции С(£,р) (нормированное на площадь области число пар индивидов, расстояние между которыми равно £):
= щ Jр{х){р{х + С) - 5(C)) dx,
плотность троек популяции Т(£, ,р) (нормированное на площадь области число троек индивидов, в которых расстояние от первого до второго равно а от первого до третьего — £'):
•т. с,р) = щ Iр(х)(р(х + о - б(о)(р(х + с') - - б(£ - О) йх.
А
Вместо введенных пространственных моментов будем рассматривать их математические ожидания по паттерну р(х) и обозначать их Ж, С(£) и соответственно. Далее мы будем работать с интегро-дифференциальными уравнениями динамики пространственных моментов, вывод которых достаточно громоздок (см. [2]):
^ = J C{i)w{i)di,
dC(Q dt
= bmiON + / иоск + С) < - („ + ^(0)С(0 - / Ж«, О -if.
(i)
Заметим, что уравнение динамики первого пространственного момента использует второй момент, уравнение динамики второго момента — третий. Уравнение динамики третьего пространственного момента будет использовать четвертый момент и так далее. Поэтому количество неизвестных функций (пространственных моментов) всегда будет превышать число уравнений. Для того чтобы избавиться от подобных зависимостей, было предложено "замкнуть" уравнения динамики, исходя из биологических ограничений на исследуемые функции. В настоящей работе мы использовали так называемое "параметрическое семейство замыканий" третьего момента, предложенное в [7]:
C(QC(0 C(QC(e - О C(?)C(?-Q _ з N N N
(1 — а)-jy-• (2)
Цель настоящей работы состоит в нахождении зависимости стационарного первого момента от коэффициентов эксцесса распределения ядер, т. е. в зависимости от параметров куртозианов.
3. Уравнение равновесия. Рассмотрим далее стационарную точку системы (1), т. е. точку, в которой производные первого и второго моментов обращаются в нуль. Для упрощения дальнейших вычислений будет удобно ввести ряд обозначений. Будем выделять неизвестные величины жирным шрифтом, при этом опустим аргументы у всех функций. Скалярное произведение, / а(х)Ь(х) йх будем записывать в виде (а, Ь). Свертку двух функций j' а(х — у)Ь(у) йу обозначим Кп Кп
через [а * Ц. Кроме того, положим й'ш = ш ж Ьт = т.
В новых обозначениях после подстановки замыкания (2) система (1) примет вид
(С,¿5),
О =гШ + [С * т] - ёС - шС - [С{£, С) + С[ш * С] + [шС * С] - ё'Щ + —С(ш, С).
_ _
Разделим оба уравнения последней системы на К2 и обозначим С = —^, (ш, С) = У. Тогда из первого уравнения получим
Ь — й Ь — й
Подставив данное выражение во второе уравнение системы, будем иметь
О = + [С * т] - (Ю - шС - (уС + Ср) * С] + [а?С * С] - ё') + (1 - а)(Ъ - ё)С.
О (X /а X \ /
В силу особенностей численного метода удобно,_чтобы искомая функция стремилась к нулю на бесконечности. Сделаем замену переменных: С = + 1. Тогда после упрощения получим уравнение
£ + = + } ^ а^^ ((Ц + 2)[ш + Щ (3)
где У = (ш, С} + 1). Будем называть (3) уравнением равновесия.
Прежде чем приступать к дальнейшему анализу этого уравнения, заметим, что, если оно имеет какое-либо решение, то оно имеет и радиально-симметричное решение. В самом деле, допустим 0 — решение (3). Значит, эта функция обращает уравнение равновесия в верное тождество. Тождество сохранится, если взять от его левой и правой частей усреднение по поверхности п-мерной сферы. Иными словами, для любого г > 0 верно равенство
1 /ГЛ. , а/, , ¿'(Ъ-й)
5
- ((<Э + 2)[а? * д] + т * 0]) |
где У = {а?, Я + 1), а |5Г| = пъп12гп~1 / Г + .
Введем функцию V = ¥(х) = —-—- <1) С^йБ. Далее учтем, что константу и, вообще гово-
\ьы\->
ЬЫ
ря, любую радиально-симметричную функцию можно вынести за знак интегрирования по сфере. Кроме того, используя переход к сферическим координатам, а также теорему Фубини, нетрудно показать, что интегрирование по сфере можно внести внутрь свертки. В итоге будем иметь
Остается лишь сказать, что функция V является радиально-симметричной, и взять еще раз усреднение по сфере произвольного радиуса, чтобы получить верное равенство:
Итак, функция V является радиально-симметричным решением уравнения (3). Таким образом, мы имеем право изначально искать исключительно такие решения. Для этого применим метод последовательных приближений или метод Неймана, суть которого в данном случае заключается в следующем.
Введем оператор
+ + 2 )[«*/]+ [£/*/])
Будем считать, что он действует в пространстве функций нужной нам гладкости, причем это пространство полно (допустим, это множество непрерывных функций с метрикой p(f,g) = = sup|/(a;) — д(х)|). Зададим какое-либо начальное приближение Q0 (например, ш или т). Далее каждое Qn, п= 1, 2, 3,..., вычисляется по формуле Qn = ДС^^.
Принцип сжимающих отображений гласит, что, если для любых рассматриваемых нами функций / и д имеет место p(f,g) ^ cp(Af,Ag), где 0 < с < 1, т.е. если оператор Л является сжимающим, то последовательность сходится к решению уравнения *4Q = Q по метрике р. Построим дискретный аналог метода последовательных приближений.
В определении оператора Л присутствует интегрирование по всему пространству Шп, однако
функции ш,т — интегрируемы в несобственном смысле, а функция Q -> 0, т. е. мы можем
г—Ц-оо
выбрать область достаточно большого радиуса (область, граница которой достаточно отделена от нуля) и проводить интегрирование только по ней в силу того, что интеграл по внешности этой области будет мал. В качестве такой области, в зависимости от замены переменных, можно выбрать куб или шар с центром в нуле. Кроме того, проведем дискретизацию, а именно, выберем сетку с достаточным количеством узлов и будем рассматривать все функции только на ней.
При сделанных допущениях оператор Л переходит в дискретный оператор В, задаваемый формулой, аналогичной (4). Вопрос о том, является ли данный оператор сжимающим, — предмет глубокого анализа, выходящего за рамки данной статьи. Анализировать сходимость процесса итераций мы будем лишь численно.
4. Ускорение вычисления многомерной свертки. В данном разделе мы будем пользоваться идеями из работы [8]. Как видно из уравнения (3), наиболее вычислительно сложная часть — это вычисление свертки двух функций. "Наивный" вариант алгоритма подразумевает сложность 0(iV2), где N — количество узлов сетки. Однако приведенная ниже теорема о свертке позволяет упростить алгоритм.
Прежде всего введем понятие преобразования Фурье в n-мерном случае. В нашей статье мы будем пользоваться формулами:
= J/(ж1,ж2, • • • ,хп) ■ ехр ^-г ^ XjUj^j dx i dx2 • • • dxn,
Rn J"1
Т~1д(ш1,ш2, ...,шп)= J g(xi,x2, ...,xn)- ехр | i^Xjioj J dx i dx2--- dxn.
Rn \ J"1 /
Отметим, что все рассматриваемые в статье функции имеют фурье-образы, суперпозиции которых имеют фурье-оригиналы.
Теорема (о свертке). Пусть функции fug принадлежат пространству Li(Wl). Тогда верна формула
If *9}(х 1
, X2i • • •, хп
) = Т 1{37/(ш1,ш2, .. .,шп) • 37д(ш1,ш2,... ,ШП)}(Ж1 , х2, • • • , Хп ) .
Как видно, данная теорема позволяет свести операцию взятия свертки к трем преобразованиям Фурье (двум прямым и одному обратному). Чтобы понять важность этого факта введем прямое
и обратное дискретные преобразования Фурье (ДПФ). Они задаются формулами
= £ Л ехР (-' = ^ X Р3 ехР • (5)
к=0 4 ' 3=0 4 '
Здесь N — количество узлов дискретной равномерной сетки. Несмотря на то, что формулы (5) имеют квадратичную вычислительную сложность, существует алгоритм быстрого дискретного преобразования Фурье (БПФ), сложность которого 0(Ж1пЖ). В итоге получаем, что сложность свертки в случае Ж1 составляет 0(Ж1пЖ).
Однако в пространствах большей размерности даже применение теоремы о свертке не позволяет снизить сложность вычислений до приемлемого уровня.
5. Случай Ж2. Прямое и обратное дискретные преобразования Фурье в двумерном случае выглядят так:
Ртп = ^ еХР ( - (кт + 3П) ) > /¡У = 2 Етп еХР \ (кт + ^ / '
Сложность такого алгоритма 0(Ж4), поскольку нам необходимо найти Ы2 значений в каждом узле квадратной сетки. Однако формулы (6) можно преобразовать:
Ргпп }кз ехР ( - ~^кгп
к=0 j=0 ^
т=0 *-п=0 ^
ехр ( ——:1»
2т .
V'
(7)
2т . ехр | — ]п
Сравнивая (7) с (5), заметим, что мы свели вычисление двумерного преобразования Фурье к вычислению одномерного. Таким образом, использование алгоритма быстрого преобразования Фурье позволит получить вычислительную сложность О (Ж3 1п Ж). Далее попытаемся ускорить алгоритм.
Вернемся к непрерывному преобразованию Фурье:
к2
Будем считать, что функция / является радиально-симметричной, т.е. /(ж,у) = /г(г), г = \/х2 + у2. Данное предположение не выводит нас за рамки задачи, поскольку мы рассматриваем только радиально-симметричные ядра и решения, исходя из биологических соображений и рассуждений, проведенных выше.
Перейдем к полярным координатам:
= рсозф, Ш2 = рвтгф, х = гсо$(р, у = гзтср. В терминах новых координат преобразование Фурье имеет вид
Однако
"ТОО 7Г
I Гк(г) ¿Г I е-"-рсо8«>-¥>) ¿(рш
0 -7Г
7Г 7Г
I е-ггрсо8(^-^) 4ср= I е-*грсо^<1ср = 2жМгр),
где </о(ж) — функция Бесселя первого рода нулевого порядка. В таком случае получаем следующий вид для преобразования Фурье в полярных координатах для радиально-симметричной функции
/(ж, у) = h(r):
Hf}{PiÍJ) = 27r I rh(r)J0(rp)dr.
Остается лишь заметить, что образ функции / также не зависит от угла, т. е. является радиально-симметричной функцией. Таким образом, мы получили преобразование Ханкеля нулевого порядка. В общем случае данное преобразование имеет вид
+ ОС
Ш(р)= I г/(гШгр)ёг, (8)
о
где За — функция Бесселя первого рода порядка а. Обратное к этому преобразование задается формулой
+ ОС
П l[g]{r)= j pg{p)Ja{rp)dp.
Итак, преобразование Фурье в Ж2 в случае радиально-симметричных функций сводится к преобразованию Ханкеля в Ж1. Сделав в (8) экспоненциальную замену переменных вида г = г0еж, Р = Pqcv; получим следующую формулу:
+ ОС
Щ)(Роеу) = I r2Qe2xf(rQex)JQ(rQpoex+y) dx. (9)
—сю
Видно, что (9) это свертка двух функций, причем всего лишь одномерная; сложность ее вычисления 0(N liaN). Итак, мы можем вычислить двумерную свертку, используя теорему о свертке и двумерные преобразования Фурье, каждое из которых с помощью преобразования Ханкеля сводится к вычислению одномерной свертки. Таким образом, в условиях, рассматриваемых в статье, сложность вычисления двумерной свертки такая же, что и у одномерной.
6. Случай Ж3. ДПФ в трехмерном случае задается аналогично двумерному, поэтому мы опустим соответствующие формулы. Алгоритм также можно ускорить, сведя его к одномерному БПФ, получив сложность 0(N4\iaN). Покажем, что в случае радиально-симметричных функций можно получить значительный выигрыш в скорости. Введем функции
где Pn(t) — присоединенные полиномы Лежандра:
/ i\n jn+k
Заметим, что замечательным свойством полиномов Лежандра, является их ортогональность. Теперь воспользуемся разложением функции ег(ш,г\ где ш = {wi, шг, шз}, г = {ж, у, z}, по функциям Y* в сферических координатах [9]:
+ СЮ п I-
= (10) п=0 — п '
Здесь
X = Г COS sin 0, у = г sin sin 0, z = r COS0,
ш\ = р eos ф sin г], Ш2 = ps'mi/js'mr], u)^ = pcosr¡.
После подстановки (10) в теорему о свертке для радиально-симметричных функций /(г, (р, в) = = 1г(г) и д(г, (р, в) = и(г) получим
атвЛгЛр Лв
и(г)е «т в (1г ¿ср в,в )
(Лш-Г)р2 хтг}(1р (1/ф ¿7] =
+ ОС П I-
/г(г)4тг X £(-1)" ф)ГЦ{9, V)г2 нт в Лг Лр Лвх
п — П — п. V
Iя Ж3хЖ3
+ос п
п = 0 — и
Т V. и I-
х«(г)4тг X £(-1)" чУ нт в Лг Лр Лв
п = 0 —п '
ег(ш'г)р2 нт г] йр (Ц) в,г].
Меняя порядки суммирования под интегралами, приходим к выводу, что ненулевое слагаемое остается только одно. В итоге формула сворачивается к следующему виду:
[/ * <7]ж* = 4тг[г/г(г) * «(г)]Ж1.
Итак, мы вновь получили вычислительную сложность О {И 1пЖ).
7. Результаты численных экспериментов. Численные решения были получены для ядер-куртозиан, связанных следующими условиями: ,в'0" = = ,ч0, = = ,Ч1. Параметр замыкания а был выбран равным 0.4, как самый оптимальный, что было показано в [7]. В данной статье мы не будем подробно рассматривать зависимость решения от варьирования параметра а.
Приведем решения для случая ,ч0 = ,Ч1 = 1 в одномерном, двумерном и трехмерном случаях для Ь = = 1 (см. рис. 1). В силу радиальной симметричности функций второго пространственного момента приведем только зависимость плотности парных корелляций от расстояния до индивида. Видно, что с увеличением размерности пространства "пиковость" второго момента уменьшается. С биологической точки зрения это означает увеличение радиуса взаимодействия особей с уменьшением при этом активности данного взаимодействия. Такая тенденция сохраняется и при других наборах параметров.
Далее рассмотрим зависимость первого момента (средней плотности индивидуумов) от параметров (куртозианов) ядер. Для визуализации результатов построим поверхности, задаваемые
6
Рис. 2
функциями N = N(.40..Ч1), в пространствах трех различных размерностей, случай Ш (рис. 2,а), 2Б (рис. 2,6), ЗБ (рис. 2,в).
Видно, что чем выше размерность пространства, тем круче поднимается поверхность вблизи начала координат, однако, во всех случаях при увеличении параметров ядер в совокупности график функции N выходит на стационарное "плато" (в некоторых случаях до полного вымирания сообщества). Таким образом, можно сделать вывод, что количество особей рассматриваемого вида после выхода на стационар тем больше, чем меньше по величине параметры ,ч0 и ,41, т.е. чем меньше коэффициент эксцесса ядер.
Благодарим У. Дикмана за постановку задачи, а также внимание к нашей работе. Авторы благодарят также А. Нестеренко, А. Бодрова и А. Савостьянова за неоценимую помощь в подготовке рукописи и ценные обсуждения.
СПИСОК ЛИТЕРАТУРЫ
1. Dieckmann U.. Law R. Moment approximations of individual-based models // The Geometry of Ecological Interactions: Simplifying Spatial Complexity / Ed. by U. Dieckmann. R. Law. J. Metz. Cambridge University Press. 2000. P. 252 270.
2. Dieckmann U.. Law R. Relaxation projections and the method of moments // The Geometry of Ecological Interactions: Simplifying Spatial Complexity / Ed. by U. Dieckmann. R. Law. J. Metz. Cambridge University Press. 2000. P. 412 455.
3. Bolker В.. Pacala S. Using moment equations to understand stochastically driven spatial pattern formation in ecological systems // Theor. Population Biol. 1997. 52. N 3. P. 179 197.
4. Бодров А. Г.. Никитин А. А. Качественный и численный анализ интегрального уравнения, возникающего в модели стационарных сообществ // Докл. АН. 2014. 455. № 5. С. 507 511.
5. Кали с трат о ва А. В.. Никитин А. А. Исследование уравнения Дикмана с интегральными ядрами, имеющими переменное значение коэффициентов эксцесса // Докл. АН. 2016. 470. № 6. С. 628 631.
6. Seier Е.. Bonett D. Two families of kurtosis measures // Metrika. 2003. 58. P. 59 70.
7. Murrell D. J.. Dieckmann U. On moment closures for population dynamics in continuous space //J. Theor. Biology. 2004. 229. P. 421 432.
8. Бодров А.Г., Никитин A.A. Исследование интегрального уравнения плотности биологического вида в пространствах различных размерностей // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2015. № 4. С. 7-13. (Bodrov A. G., Niki tin A.A. Examining the biological species steady-state density equation in spaces with different dimensions // Moscow Univ. Comput. Math, and Cybern. 2015. 39. N 4. P. 157-162.)
9. Baddour N. Operational and convolution properties of three-dimensional Fourier transforms in spherical polar coordinates //J. Opt. Soc. Am. A. 2010. 10. P. 2144-2155.
Поступила в редакцию 31.01.18
УДК 519.233.24, 519.233.5 А. Г. Белов1
МОДЕЛИРОВАНИЕ СОВМЕСТНОЙ ДОВЕРИТЕЛЬНОЙ ПОЛОСЫ СРЕДНЕГО ЗНАЧЕНИЯ ПОВТОРНЫХ ОТКЛИКОВ С ПРЯМОУГОЛЬНОЙ ОБЛАСТЬЮ ДЛЯ ПРЕДИКТОРОВ
В статье рассмотрена задача моделирования совместных доверительных интервалов для среднего значения повторных откликов в линейной множественной нормальной регрессионной модели с предикторными переменными, определенными в интервалах. Для ее решения применен численный метод вычисления критического значения, определяющего совместный доверительный интервал заданного уровня. Проведено численное моделирование и сравнительный анализ совместных доверительных интервалов для регрессии, среднего значения повторных откликов и отдельного наблюдения.
Ключевые слова: совместные доверительные интервалы, нормальная регрессия, повторные отклики.
1. Постановка задачи. Рассмотрим линейную множественную нормальную регрессионную модель наблюдений:
у = Х/3 + е,
где у = (?/1,... ,уп)т — вектор-столбец случайных величин (с. в.) у^ откликов, описывающих результаты г-го опыта, е = (в1,... ,еп)т — вектор-столбец случайных "ошибок" с нормальным законом распределения £(е) = Л/"„(0, с21п), не зависящий от вектора параметров /3 = (/?!,... ,/3/;)т; X = Цх^,... , х^|| € Дпх/г — регрессионная матрица из вектор-столбцов х^ = (жц,... ,Хщ)т, оказывающих влияние только на среднее значение отклика Еу^, при этом 1п = diag(l,..., 1) € ейпхп, гапкХ = &, к^п.
Пусть имеется т повторных наблюдений ут = ..., ут)Т\ соответствующих фиксированным значениям регрессоров х = (ж1,..., у-} = х; /3 + г/. 1 ^ ] ^ т, где вектор-столбец случайных "ошибок" ет = (в1,... ,ето)т не зависит от е и £(ето) = Л/"то(0, ст21то), а £(етоо) = Л/"то(0, 1то) для ето0 = £то/ст.
Для среднего значения повторных откликов ут = е^ут/т = хт/3 + (те^£тоо/т, где ето = = (1,..., 1)Т € Дто, используется 100(1 — а)%-ъ доверительный поточечный интервал [1]
у =Fii-f ^ + xTA"1xj , (1)
где А = ХТХ, у = хт/3 — оценка отклика у для х, /3 = А_1Хту = /3 + ctA_1Xt£q — оценка вектора параметров /3, найденная по выборке у с помощью метода наименьших квадратов (МНК),
1 Факультет ВМК МГУ, ст. науч. сотр., к.ф.-м.н., e-mail: belovQcs.msu.ru