CC BY
24
52. Абакумов М.В., Жданов А.А., Мухин С.И., Попов Ю.П., Чечеткнн В.М. Аккреционные диски в двойной звездной системе с переменным расстоянием между компонентами. Математическое моделирование. Препринт ИПМ им. М. В. Келдыша РАН. 2003. № 65.
53. Абакумов М.В., Мухин С. П., Попов Ю.П. О некоторых задачах гравитационной газовой динамики // Матем. моделир. 2000. 12. № 3. С. 110-120.
Поступила в редакцию 20.06.04
Д. П. Костомаров
ОТ ЧИСЛЕННОГО РЕШЕНИЯ КЛАССИЧЕСКИХ ЗАДАЧ МАТЕМАТИЧЕСКОЙ ФИЗИКИ К МОДЕЛИРОВАНИЮ РЕАЛЬНОГО МИРА
(кафедра автоматизации научных исследований факультета ВМиК)
1. Введение. Первые работы, посвященные компьютерному моделированию различных физических явлений, относятся к началу 50-х гг., когда военная наука потребовала быстрого решения однотипных вычислительных задач, относящихся, как правило, к математической физике. Потребности практики стимулировали, с одной стороны, развитие теории дифференциальных уравнений, с другой — создание вычислительных алгоритмов, ориентированных на компьютерные расчеты. Важное значение приобрела "сертификация" моделей, степень соответствия результатов расчета экспериментальным данным. Так в науке появились новое направление — вычислительная физика и новый термин — "вычислительный эксперимент".
На первом этапе своего развития вычислительная физика рассматривалась как раздел вычислительной математики, в котором основное внимание уделялось численному решению различных задач для уравнений математической физики. При этом, поскольку компьютеры обладали весьма низкой производительностью, исследователи обычно ограничивались решением одномерных задач на довольно грубых сетках. Это привело к распространению точки зрения, что вычислительная физика лишь позволяет уточнить результаты, полученные аналитическими методами для упрощенных моделей. В некотором смысле вычислительная физика рассматривалась как служанка теоретиков без права на собственное мнение. Но даже в рамках такого ограниченного подхода были получены впечатляющие результаты. В качестве примера можно привести предсказание существования Т-слоя в высокопрово-дящей жидкости, сделанное Тихоновым и Самарским и их учениками. Отметим еще один интересный факт. Вычислительная практика стала существенно опережать математическую теорию вычислительных алгоритмов. Не умаляя достижений теоретической вычислительной математики, следует признать, что реально используемые и доказавшие на практике свою высокую эффективность алгоритмы строгое теоретическое обоснование либо получили с опозданием, либо не получили вовсе. Буквально на глазах одного поколения исследователей вычислительная физика стала в значительной степени экспериментальной наукой.
Тем временем компьютеры развивались в соответствии с законом Мура. Вычислительные машины все больше и больше превращались из "вычислительных мельниц" в универсальные устройства обработки информации самого общего характера, не только числовой, но и графической, текстовой и т.д. Это открыло блестящие перспективы и перед вычислительной физикой, которая из Золушки превратилась в Принцессу. Бурное развитие компьютеров поставило на повестку дня совершенно новые и часто даже неожиданные задачи. Во-первых, стало ясно, что можно моделировать окружающий мир исходя из первичных принципов, а не ограничиваться численными решениями дифференциальных уравнений, которые сами были получены в рамках достаточно существенных ограничений. Фактически мы можем рассматривать некоторую искусственную "компьютерную" материю, поведение которой удается проследить достаточно точно и подробно. Именно сейчас мы можем проводить настоящий вычислительный эксперимент и сравнивать его результаты как с теоретическими представлениями, так и с данными реального опыта. При этом возникает очень важная методологическая и практическая проблема интерпретации результатов вычислительного эксперимента. Во-вторых, потребовалось развить саму технологию проведения вычислительного эксперимента. В частности, это
подразумевает создание удобного интерфейса для ввода исходных данных не только в числовой, но и в аналитической форме. Но еще более важным является представление результатов расчетов. Оно потребовало использования всего современного арсенала компьютерной графики, а также разработки уникальных методов вывода графической информации, что является весьма актуальным особенно для трехмерных задач. В-третьих, необходимо как можно более эффективно использовать имеющиеся компьютеры, что, естественно, ставит вопрос о широком использовании параллельных вычислений и суперкомпьютеров. В-четвертых, требуется постоянно развивать и совершенствовать сами математические модели, а также численные методы их исследования, в частности предназначенные для параллельных вычислений. Наконец, в-пятых, создание нового инструментария позволило несколько по-иному взглянуть на задачи из тех областей, которые на первый взгляд имеют мало общего с физикой.
На кафедре автоматизации научных исследований факультета ВМиК ведутся работы по многим аспектам современной вычислительной физики. Традиционными являются работы по математическому моделированию высокотемпературной плазмы, ориентированные на решение проблемы управляемого термоядерного синтеза (УТС). У истоков этого направления в Московском университете стояли академики И. В. Курчатов и А.Н. Тихонов. Оказалось, что плазма является весьма благодатным объектом для физиков-вычислителей. Наличие большого разнообразия явлений и в то же время достаточно простые физические принципы описания этих явлений привели к созданию целой иерархии моделей и исследованию с их помощью весьма тонких и интересных эффектов. Работа по этой теме на кафедре проводится в тесном сотрудничестве как с ведущими российскими учеными, так и с коллегами из крупнейших научных центров США, Великобритании, Италии и других стран.
В этом обзоре, подготовленном по случаю 250-летнего юбилея Московского государственного университета, мы опишем лишь некоторые, наиболее значимые результаты, полученные сотрудниками кафедры за последние 5 лет. Обзор состоит из четырех разделов. Первый раздел посвящен моделированию высокотемпературной плазмы в связи с проблемой управляемого термоядерного синтеза [1]. Этим направлением мы начали заниматься задолго до создания кафедры. Направление, описанное во втором разделе, первоначально развивалось в рамках первого, но потом "отпочковалось" в качестве самостоятельного. Третий и четвертый разделы посвящены исследованиям в нетрадиционных для нас областях, которые на кафедре только начинаются. Описанные работы были поддержаны грантами РФФИ 02-01-00299, 99-01-01167, НШ-1349.2003.1, "Поддержка научных школ" № 00-15-96052, программой "Университеты России" № УР 03.03.035, 015.03.002.008, INTAS-00233. В конце обзора приведен список наиболее важных публикаций. При подготовке обзора мне активно помогали многие сотрудники кафедры. С чувством глубокой благодарности к коллегам перечислю их имена: Е.Ю. Еч-кина, Ф. С. Зайцев, И. Н. Иновенков, A.B. Леоненко, A.M. Попов, А. П. Смирнов, Д.Ю. Сычугов, А. Б. Шмелев.
2. Моделирование высокотемпературной плазмы в установках токамак
2.1. Трехмерное моделирование неоклассических неустойчивостей в токамаке. Важной характеристикой эффективности удержания плазмы в установках токамак является параметр /3, равный отношению давления плазмы к давлению магнитного поля. Желательно сделать его как можно больше, однако увеличению ß выше некоторого порогового значения препятствуют неустойчивости. Среди них особенно выделяются неоклассические тиринг-моды (NTM). При подъеме давления плазмы выше порогового тиринг-неустойчивость приводит во многих случаях к срыву разряда. Неустойчивость является нелинейной: она связана с взаимодействием винтовых мод и требует для своего исследования решения полной системы уравнений магнитной гидродинамики (МГД). Был разработан трехмерный МГД-код NFTC (Nonlinear Full Toroidal Code), позволяющий моделировать нелинейные неустойчивости плазмы в токамаке в реальной геометрии [2-6]. Математическая модель основана на МГД-уравнениях, включающих неоклассические эффекты: возмущение бутстрэп-тока, транспортную модель с анизотропной теплопроводностью, ток поляризации. Неоклассические члены включены в основные уравнения для магнитного поля и давления. Эффективная неявная схема позволяет транспортным профилям самосогласованно эволюционировать в присутствии нелинейных МГД-неустойчивостей и внешних источников магнитного поля и давления.
С помощью кода NFTC найдены условия возбуждения пороговых неустойчивостей неоклассических тиринг-мод для профилей плазмы, полученных в ряде экспериментов на токамаке DIII-D (Düblet III — D-Shape, США). Построена диаграмма нелинейных неустойчивостей в зависимости от значений
коэффициента запаса q и обнаружены области параметров, позволяющие получать устойчивые разряды (щели устойчивости) [2]. Результаты моделирования подтверждены прямым сравнением с экспериментами на установке DIII-D [2-4]. Исследованы условия стабилизации тиринг-неустойчивости с помощью дополнительных токов, которые создаются инжекцией электронно-циклотронных волн (Electron Cyclotron Current Drive — ECCD-токи) [3, 7]. Определены оптимальные параметры радиально локализованного ECCD-тока и найдены согласованные режимы генератора волн, обеспечивающие подавление неустойчивости на нелинейной стадии. Проведено моделирование различных разрядов как с монотонным профилем тока, так и в перспективных разрядах с отрицательным магнитным широм.
Проведено моделирование возникновения неустойчивости NTM в проектируемом токамаке-реак-торе ITER (International Thermonuclear Experimental Reactor). Найден эффект двухпорогового возбуждения неоклассических тиринг-мод в плазме с реакторными параметрами [8, 9]. Высокая температура плазмы в ITER приводит к необходимости проведения вычислений с экстремально большим магнитным числом Рейнольдса Rem ~ Ю10 и требует высокой разрешающей способности кода при моделировании узких внутренних слоев. Показано, что малый порог магнитных возмущений не является опасным для ITER и приводит к длительным стационарным состояниям с магнитными островами, которые могут быть стабилизированы обратными связями. Получены условия стабилизации нелинейных МГД-неустойчивостей в ITER с помощью ECCD-тока.
Для проектирования устойчивых профилей и проведения вычислений для плазмы с реакторными параметрами создан параллельный комплекс программ на основе кода NFTC для компьютера IBM 690 p-series Regatta, установленного на факультете ВМиК МГУ. Проведены исследования начальной стадии турбулентности при экстремально больших магнитных числах Рейнольдса на вычислительных сетках с мелким шагом.
2.2. Разработка стандартных кодов равновесия и устойчивости плазмы в установках токамак. В связи с ведущимся в настоящее время проектированием установок токамак нового поколения особую важность приобретают вопросы доведения до уровня стандартных программ классических кодов расчета МГД-равновесия, МГД-устойчивости и магнитной диагностики плазмы. Такая работа проводится с кодами, которые на протяжении нескольких десятилетий были написаны на кафедре автоматизации научных исследований. В основном стандартизация проводится по следующим направлениям. Коды должны позволять детально учитывать все элементы конструкции установки, что важно как при проектировании, так и при сопоставлении результатов расчетов с экспериментальными данными. Например, при выборе магнитной системы внешних полей приходится принимать во внимание инженерные ограничения на расположение катушек, их размеры, предельную плотность тока. При анализе устойчивости плазмы приходится более детально учитывать геометрию и материал стабилизирующих элементов. Поскольку проектирование является итерационным процессом, то оно требует проведения многократных серийных расчетов. Поэтому важной задачей является ускорение алгоритмов. Еще более жесткие требования предъявляются к алгоритмам диагностики, которые являются частью комплекса программ по управлению экспериментом в режиме реального времени. Следующим важным аспектом является разработка таких версий кодов, которыми могли бы пользоваться не только математики, но также инженеры и экспериментаторы. Для этого необходимо сводить количество "внутренних" вычислительных параметров к минимуму. Одной из важных проблем является проблема интеграции. В настоящее время в мире существует уже немало кодов, которые всеми признаны и используются как стандартные. Поэтому новый код может быть принят сообществом как стандартный только при наличии совместимости с уже имеющимся "программным полем". И наконец, остается актуальной разработка средств сервиса. Сервис включает в себя наличие руководств по применению, а также встроенной документации и средств помощи. Обязательным является либо наличие собственных средств визуализации, либо возможность визуализации данных при помощи стандартных графических пакетов. На стадии конструирования огромную помощь может оказать встроенный графический редактор входных данных. В начале 2000-х гг. сотрудниками кафедры совместно с РНЦ "Курчатовский институт" был разработан и внедрен в качестве стандартного современный высокоэффективный код TOKAMEQ (TOKAMak EQuilibrium) для расчета МГД-равновесия и МГД-устойчивости плазмы, а также контроля границы плазменного шнура по результатам магнитных измерений. Код ориентирован на физиков, проектировщиков, экспериментаторов и инженеров и обладает следующими достоинствами:
1) возможностью расчета МГД-равновесия с любым распределением плотности тока по сечению
шнура. Это обстоятельство позволяет в ходе расчета сценария разряда обмениваться информацией с транспортными кодами;
2) высокими скоростными характеристиками. Время расчета МГД-равновесия и МГД-устойчиво-сти на ПЭВМ Pentium IV составляет порядка 3-5 с, восстановление границы плазменного шнура по результатам магнитных измерений требует нескольких секунд, что достаточно для контроля плазмы в режиме реального времени;
3) высокоэффективным текстово-графическим редактором данных, позволяющим вести проектирование установки и интерпретацию эксперимента в интерактивном режиме;
4) возможностью использования как встроенных средств визуализации, так и обращения к стандартным графическим пакетам;
5) возможностью автоматической настройки вычислительных параметров;
6) эффективным сервисом, включающим в себя удобный интерфейс, средства самодокументации (помощь, встроенные описания кода, инструкции по его эксплуатации), средства обмена информацией с другими кодами, средства поддержки диалога на русском и английском языках.
Для поддержки кода не требуется никаких специальных программных средств. С помощью кода были просчитаны основные сценарии разряда для проектируемой в России установки токамак Т15-М [10, 11]. Проведен анализ МГД-устойчивости плазмы, в том числе и при конечной проводимости стабилизирующих элементов. На основе расчетов была выбрана конструкция витков пассивной обратной связи для подавления вертикальной МГД-неустойчивости плазмы. Исследована управляемость плазменного шнура при помощи системы активных обратных связей. Совместно с сотрудниками UKAEA (Великобритания) проводился расчет магнитной системы для проектируемой в Великобритании установки CTF (Components Test Facility). Удалось сконструировать систему проводников, удовлетворяющую одновременно инженерным ограничениям и требованиям сценария разряда [12].
2.3. Математическое моделирование эволюции равновесия плазмы в установках токамак. В тороидальной плазме развивается большое количество процессов с различными характерными временами. Важную роль играют относительно медленные процессы перехода одного квазиравновесного состояния в другое, вызванные самосогласованным взаимодействием плазмы и электромагнитного поля. Такие процессы определяют ключевые характеристики плазмы, параметры режимов улучшенного удержания, величину возникающих в плазме дополнительных электрических токов. Характерные времена этих процессов имеют масштаб длительности разряда в установке.
Модель эволюции равновесия тороидальной плазмы подробно описана в работе [13]. В качестве исходных рассматриваются уравнения Максвелла, в которых опущен ток смещения, закон Ома и уравнение баланса газодинамического и магнитного давления
^^хВ. (1)
Закон Ома характеризует свойства среды, его записывают с разной степенью детализации. Для плазмы токамака используется закон Ома в следующем виде:
Е + Еаси + V х В = — + —, (2)
(Гц <7_|_
где Е — электрическое поле; <7ц, ^ и <т_|_, — соответственно проводимость плазмы и компоненты тока,], параллельные и перпендикулярные магнитному полю В; V — скорость движения плазмы. С помощью члена Е^ в формуле (2) явно выделены дополнительные токи:
_ Ja.dc!,11
-Ь<ас1с1 =--1--•
(Гц <7_|_
Они генерируются градиентом давления р (бутсрэп-ток, ток Пфирша-Шлютера, диамагнитный ток), ВЧ-волнами, инжекцией нейтралов, меняющимся магнитным полем и т.д.
В случае аксиальной симметрии задача эволюции (1), (2) сводится к системе двух сильнонелиней-
ных уравнений для двух неизвестных функций ip(t,R,Z) и F(t,ip),
„3/1 дф\ д2ф
RdR \RdR Г = (3)
jr,{RA) =
RMhÉÏ i i 9F* м) D (f]
л дф ^ 2n0R Оф ь
Ji(t)6(R - Ri)6(Z - Zi) вне Dp(t),
i = 1
(4)
1 -^F -д£дф\__( RJL( _
dt\J )dZ~dt ~ fi0(T\\\dRdR + dZ ÔZ I ~ fi0(T\\ V dR \R dRj + dZ2 ) ~
о v 7
R2
--jadd-B В Dp(t), (5)
<71|
,M) = f ÏM&v/f JL. ,6,
J pol / J ф=const n\>o\
ф=const
Уравнения (3)-(6) записаны в цилиндрических координатах R, г], Z\ Dp — область, занятая плазмой; Гр — ее граница. Искомые функции фи F определяют тороидальную и полоидальную составляющие магнитного поля В:
Btor = (F/R)^, Вро1 = [Уф х ij. (7)
Через Ji(t) в формуле (4) обозначены токи, протекающие вне плазмы: в соленоиде, в катушках поло-идального поля, в стенках камеры. Усреднение в (6) проводится по линии пересечения тороидальной магнитной поверхности ф = const с меридиональной плоскостью г] = const. Уравнения (3)-(6) следует дополнить начальными и граничными условиями
if)(0,R,Z) = il)o(R,Z), F(0, ф) = Fo(V0, lim ^{t,R,Z) = lim ф{Ь, R,Z) = 0,
R-> 0 й->со
Z—>oo
F(t,R,Z) =^/tor(i), (8) rp(t) Zn
где /tor(i) — электрический ток, протекающий через центральный сердечник и создающий тороидальное магнитное поле, z) — решение уравнения (3) с заданными начальными функциями р(0,ф) и Fo(ij)). Система уравнений (3), (5) достаточно своеобразна. Она включает стационарное уравнение (3), представляющее условие равновесия тороидального плазменного шнура (уравнение Греда-Шафранова), и эволюционное уравнение (5), содержащее временные производные обоих искомых функций. Граница плазмы является свободной. Она определяется в процессе решения задачи и изменяется со временем. Модель эволюции (3)-(6) реализована в виде кода SCoPE (Self Consistent Plasma Evolution), код имеет интерфейс на языке JAVA, позволяющем задавать входные данные, выполнять расчеты, отображать их результаты графически, преобразовывать данные в формат NetCDF.
С помощью кода SCoPE исследованы условия возникновения и способы стабилизации режимов улучшенного удержания плазмы. Теоретические предсказания подтверждены экспериментально на токамаках START и MAST (Mega-Amper Spherical Tokamak, UKAEA, Великобритания) [14-16]. Методом математического моделирования впервые изучена проводимость плазмы в сферическом токамаке MAST [17]. Выполнен объемный вычислительный эксперимент для различных моделей проводимости. Рассчитаны основные характеристики плазмы. Проведено их сравнение с измерениями в натурном эксперименте. Близость экспериментальных и численных данных, физическая интерпретация наблюдаемых эффектов убедительно свидетельствуют о том, что проводимость плазмы в установке MAST соответствует неоклассической теории. Относительная универсальность математической модели позволяет экстраполировать полученный результат и на другие сферические токамаки.
Многочисленные трудности, возникающие при решении столь сложных задач удалось преодолеть благодаря применению методики нейронных сетей, что позволило существенно упростить анализ и визуализацию многомерных данных. В результате были решены задачи об оптимизации длительности разряда и величины дополнительного тока, вызванного градиентом давления; задачи о классификации разрядов по степени устойчивости, по видеоизображению плазмы; задачи отображения многомерных
областей устойчивости плазмы на плоскость [18]. Разработанные подходы на базе нейронных сетей оказались полезными и в других областях, особенно в тех случаях, когда применение классических методов невозможно.
2.4. Разработка кода GENRAY для расчета распространения и поглощения электромагнитных волн в тороидальной плазме в приближении геометрической оптики. В установках УТС используются различные типы электромагнитных волн для нагрева плазмы, генерации в ней тока, диагностики. Одним из эффективных и широко используемых методов расчета соответствующих электромагнитных полей является метод геометрической оптики, основанный на описании электромагнитных волн с помощью лучевых траекторий. Для их расчета используется система га-мильтоновых уравнений для координат и компонент волнового вектора. Гамильтонианом является дисперсионная функция, равная детерминанту уравнений Максвелла для амплитуды периодического электромагнитного поля. Код GENRAY может использовать несколько видов диэлектрического тензора: тензор холодной плазмы, тензор тепловой нерелятивистской плазмы и релятивистский тензор для электронной плазмы. Это позволяет рассчитывать лучевые траектории для всех типов электромагнитных волн в плазме, включая электронно-циклотронные, нижнегибридные, быстрые магнитно-звуковые, бернштейновские, ионно-циклотронные.
При прохождении электромагнитной волны около зоны циклотронного резонанса антиэрмитова часть диэлектрического тензора становится сравнима с его эрмитовой частью. При этом дисперсионная функция, полученная с помощью только эрмитовой части диэлектрического тензора, может давать неверные лучевые траектории. В этом случае, как показано в ряде работ, в качестве дисперсионной функции можно использовать реальную часть собственного числа дисперсионного тензора. В коде GENRAY реализована возможность использовать такую дисперсионную функцию.
Магнитное поле задается в виде, позволяющем описывать реальную геометрию установок тока-мак. Профили плотности и температуры компонент плазмы могут задаваться в виде, соответствующем экспериментальным данным. Вдоль лучевых траекторий вычисляются поляризация и поглощение энергии электромагнитной волны и плотность тока, создаваемая в плазме поглощенной электромагнитной волной. Код может вычислять эмиссию электронно-циклотронных и электронных бернштей-новских волн из плазмы. Для этого решается уравнение переноса излучения волн лучевых траекторий. В расчетах может использоваться произвольная релятивистская функция распределения электронов.
Код GENRAY разрабатывался в сотрудничестве с фирмами Дженерал Атомикс, СотрХ, Вискон-синским и Принстонским университетами (США). С его помощью проведены расчеты по моделированию электронных бернштейновских волн в установках DIII-D (Дженерал Атомик) [19], MST (Madison Symmetric Torus, Висконсинский университет) [20, 21], NSTX (National Spherical Torus Experiment, Принстонский университет) [22] и ионно-циклотронных волн в DIII-D [23].
2.5. Усредненное кинетическое уравнение и его применение к моделированию процессов в плазме токамака. Одной из основных математических моделей плазмы является система кинетических уравнений для функций распределения ионов и электронов с операторами кулоновских столкновений. Кинетическое уравнение записывается в шестимерном фазовом пространстве, объединяющем геометрическое пространство и пространство скоростей. В общем виде такая модель сложна своей многомерностью и наличием эффектов с сильно различающимися характерными временами (до 10 порядков). Однако в конкретных случаях можно перейти от декартовых координат к специально подобранным фазовым переменным, связанным с интегралами движения частиц в рассматриваемых электромагнитных полях, разделить движения на "медленные" и "быстрые" и усреднить по "быстрым" переменным. В результате снижается размерность задачи, из рассмотрения исключаются "быстрые" осциллирующие процессы и выделяются "медленные" эволюционные изменения функции распределения, которые, собственно говоря, и представляют основной интерес.
В случае высокотемпературной тороидальной плазмы шесть степеней свободы распадаются на три "медленные" и три "быстрые" переменные. При этом зависимость от азимутального угла в геометрическом пространстве отсутствует. В результате после усреднения по "быстрым" движениям получается кинетическое уравнение с тремя фазовыми переменными. Различные модификации метода усреднения для кинетических уравнений рассматривались во многих публикациях. Их результаты обобщены и развиты в работе [24].
В исходном кинетическом уравнении оператор кулоновских столкновений обладает следующими фундаментальными свойствами [1]: 1) обращается в ноль для максвелловских распределений частиц
с одинаковой температурой и средней переносной скоростью; 2) сохраняет число частиц; 3) сохраняет суммарный импульс частиц; 4) сохраняет суммарную энергию частиц; 5) не увеличивает Н-функцию (аналог Н-теоремы Больцмана); 6) является эллиптическим.
Процедура усреднения не является тождественным математическим преобразованием. Несмотря на это, удалось показать, что усредненное кинетическое уравнение сохраняет в некоторой ослабленной форме все шесть перечисленных выше свойств неусредненного уравнения [25], подтвердив тем самым пригодность усредненных кинетических моделей для физически корректного описания тороидальной плазмы.
Усредненное кинетическое уравнение является математической моделью, пригодной для компьютерных расчетов. Были разработаны алгоритмы его решения с помощью неявных и явных итерационных схем [26, 27]. Недавно проведено исследование свойств модифицированной схемы локальных итераций ЛИ-М. Для схемы ЛИ-М и для неявной схемы расщепления предложены параллельные MPI-алгоритмы. Исследования показали, что в некоторых случаях схема ЛИ-М имеет преимущества перед неявными схемами. Использование схемы ЛИ-М в комбинации со схемой расщепления позволяет проводить устойчивые вычисления со вторым порядком аппроксимации по времени.
Разработана и реализована в коде FPP-3D (Fokker-Plank Package) трехмерная кинетическая модель для расчета поведения частиц высокой энергии в плазме установок токамак. Проведено моделирование кинетики термоядерных альфа-частиц как в ныне действующих токамаках, так и в проектируемых реакторах. В частности, исследовано поведение альфа-частиц в установке TFTR (Test Fusion Tokamak Reactor, США) [28]. Изучена роль упругих ядерных столкновений в термоядерных дейтериево-тритиевых экспериментах на крупнейшей в мире установке JET (Joint European Torus, Euroatom) [29]. Рассчитанные данные сопоставлены с результатами измерений анализатором нейтральных частиц. Показано хорошее соответствие. Найдено теоретическое объяснение некоторых экспериментально наблюдаемых эффектов.
3. Перезамыкание силовых линий магнитного поля и общая проблема устойчивости векторных полей
3.1. Физические и математические особенности задачи. Как было отмечено в предыдущем разделе, компьютерное моделирование процессов в высокотемпературной плазме подразумевает использование моделей различного типа: магнитогидродинамических, кинетических, транспортных, адекватно описывающих те или иные особенности поведения плазмы как сплошной среды. При этом если вначале основным побудительным мотивом было желание разобраться в особенностях функционирования установок, предназначенных для решения проблемы управляемого термоядерного синтеза, то в дальнейшем детальное аналитическое и численное исследование соответствующих моделей приводило к необходимости более глубокого изучения физических и математических особенностей возникающих при этом задач. С одной стороны, такой подход потребовал более абстрактной, рафинированной их постановки, с другой стороны, зачастую оказывалось, что многие задачи физики плазмы тесно связаны с другими прикладными проблемами. Одной из них является проблема перезамыкания или пересоединения силовых линий магнитного поля.
Явление магнитного перезамыкания наблюдается при исследовании космической и лабораторной плазмы и заключается в изменении топологии магнитного поля В в сравнительно небольшой области. Его проявлениями являются быстрая диссипация магнитной энергии с ее превращением в энергию ускоренных частиц, энергию электромагнитного излучения и тепловую энергию плазмы.
Суть явления магнитного перезамыкания можно пояснить следующим образом. В идеально проводящей плазме силовые линии магнитного поля "вморожены" в плазму и движутся вместе со средой, так что изменение топологии магнитного поля невозможно. Реальная плазма является очень хорошим, но не идеальным проводником. В целом для нее почти всюду имеет место эффект "квазивморожен-ности". Исключение составляют окрестности точек, где поле обращается в нуль и "вмороженность" нарушается. В результате становится возможным перезамыкание силовых линий и как следствие изменение их топологии. Таким образом, нулевые точки магнитного поля оказываются критическими по отношению к процессу магнитного перезамыкания, причем характер процесса и его конечный результат зависят от типа критической точки.
Переведем описание процесса перезамыкания с физического языка на математический. Система МГД-уравнений, описывающая поведение плазмы в магнитном поле, является сингулярно возмущенной. Члены, содержащие вторые производные вектор-потенциала, входят в уравнение для поля с ко-
эффициентом, равным магнитной вязкости vm = с2/Аттст. Если положить проводимость плазмы а равной бесконечности, то vm =0, и система вырождается. Для изучения ряда процессов можно пользоваться упрощенной, вырожденной системой. Однако в окрестности критических точек приближение идеальной магнитной гидродинамики неприменимо и нужно решать полную МГД-систему. Учет сингулярного члена, несмотря на малое значение коэффициента vm, позволяет адекватно описать эффект перезамыкания.
3.2. Математическая постановка задачи. Рассмотрим произвольную точку х, перенесем ее для удобства в начало координат и разложим магнитное поле в окрестности этой точки по формуле Тейлора
Bi(x,t) = Bi{0,i) + djBi(x.,t)xj + djkBi(-x.,t)xjXk + ... . (9)
Введем следующие обозначения: Ь{ = Bi(Q,t) — однородная составляющая магнитного поля, Aij = = д3Вг\Хк=0 и Aijk = djkBi \Хк=о — матрицы его первых и следующих производных. Здесь и ниже dj = = d/dxj. Предполагается суммирование по повторяющимся индексам. Если однородная составляющая поля равна нулю, то начало координат оказывается критической точкой:
t>i — AijX j + Aij^X jXk + .... (10)
Компьютерное моделирование процесса перезамыкания основано на численном решении смешанной задачи для системы уравнений магнитной гидродинамики [30]. Это позволяет корректно и весьма подробно описать самосогласованную эволюцию плазмы и магнитного поля в окрестности критических точек. Для упрощения задачи предполагается, что коэффициенты переноса, характеризующие плазму: магнитная вязкость z/m, электропроводность а, коэффициент теплопроводности к, — являются постоянными. Магнитное поле описывается с помощью вектор-потенциала В = rot А, что автоматически обеспечивает его соленоидальность.
В рамках одножидкостной модели плазмы система МГД-уравнений имеет вид
^ + сМН = о, (и)
<9v \ в
—+ (vV)vj =-|Vp-[AAxrotA], (12)
дА
— = [v X rot А] + z/m А А, (13)
7-i V^ + (vV)T) +Pdivv = div(KVT) + ^r(AA)2' (14)
р=\рт. (15)
Выше мы отмечали, что при vm = 0 уравнение (13) вырождается: из него выпадают члены со вторыми производными вектор-потенциала.
Для системы уравнений (11)—(15) формулировалась начально-краевая задача в кубической области G = { — 1 ^ х ^ 1; —1 ^ у ^ 1; —1 ^ z ^ 1}. Начальные условия выбирались следующим образом. Начальное магнитное поле представлялось потенциальным, что означает отсутствие в начальный момент тока:
rot В = rot rot А = 0, j = 0, (16)
а плазма считалась однородной и неподвижной:
v = 0, р= 1, Г = 1, р=1. (17)
Функции A, v, р, Т, р (16), (17) являются стационарным решением системы уравнений (11)—(15).
Перейдем к обсуждению граничных условий. Для вектор-потенциала А они выбирались такими, чтобы на боковой поверхности куба G возбуждался электрический ток, параллельный оси z. Ток создавал МГД-волну, сходящуюся к оси z, которая выводила систему из равновесного состояния.
Для функций v, р, Т, р граничные условия ставились следующим образом: на участках границы, где плазма втекает в область G, полагалось Т = 1, р = 1, р = 2/3. На остальных участках границы использовались условия свободного вытекания.
Численное решение нелинейной задачи с четырьмя независимыми переменными оказалось на грани возможности доступной вычислительной техники, но нам повезло! В НИВЦе появился кластер, состоящий из 18 двухпроцессорных узлов, на базе Pentium III/500 Mhz, и мы получили возможность совместить приятное с полезным — освоить технологию параллельных вычислений и использовать ее для решения нашей задачи. Расчет на многопроцессорном кластере тут же поставил проблему аппроксимации временных производных. В однопроцессорном режиме мы использовали сложные полунеявные схемы с итерациями. Для параллельных вычислений, как выявили тестовые расчеты, более эффективными оказались условно устойчивые явные схемы. Необходимость расчетов с малым временным шагом окупалась простотой распараллеливания и минимизацией обмена между процессами. Останавливаться подробнее на этой теме не будем. Отметим лишь, что соответствующие вопросы мы докладывали на крупной международной конференции [31].
Двумерные МГД-возмущения в холодной плазме в линейном приближении расщепляются на две независимые моды разной поляризации — магнитозвуковую и альфвеновскую. Благодаря этому они имеют достаточно простой вид. В трехмерном нелинейном случае структура решения намного сложнее и, как показали расчеты, допускает много разных вариантов. В оригинальных публикациях для обсуждения результатов расчетов мы широко использовали рисунки [30-36]. Здесь опишем некоторые результаты словами.
Начнем с наиболее простого и распространенного случая, когда магнитное поле В имеет в начале координат нуль первого порядка. Пусть, например,
Такая конфигурация имеет сепаратрисную плоскость г = 0, перпендикулярную возбуждаемому электрическому току. В этом случае МГД-волна вызывает вращение плазмы: силовые линии закручиваются и плазма выбрасывается из расчетной области через торцевые поверхности куба С. Токовые слои при этом не формируются.
Теперь рассмотрим другую конфигурацию магнитного поля:
Здесь сепаратрисной поверхностью является плоскость у = 0, параллельная току на границе. Это изменение существенно. Оно приводит к тому, что МГД-волна, достигнув центра куба С, изменяет конфигурацию поля и формирует токовые слои, параллельные оси г.
Эти примеры показывают зависимость процесса формирования токового слоя от взаимного расположения сепаратрисной поверхности в начальной конфигурации магнитного поля и направления граничного тока.
В заключение рассмотрим магнитную конфигурацию, которая наблюдается на Солнце и порождает "спорадические" солнечные вспышки. Она имеет вид
В (ж, у, г) = (-2 ху - г2 + ж2 + 0,16)ех + (-2 уг - х2 + у2 + 0,16)еу + (~2гх - у2 + г2 + 0,16)ег. (20) У этого поля имеются две критические точки, соединенные сепаратрисой: М\ = и
В результате возбуждения на границе расчетной области электрического тока структурно-неустойчивая магнитная конфигурация (20) трансформируется в структурно-устойчивую. Происходит образование сепаратрисной поверхности, разделяющей две нулевые точки. Электрический ток локализуется в окрестности этих точек.
Исследование процесса магнитного перезамыкания связано с абстрактной математической проблемой устойчивости векторных полей. Как известно, топологию векторного поля можно описать на языке динамических систем. Из теории динамических систем следует, что структурно-неустойчивыми являются поля с вырожденными критическими точками, в которых первая отличная от нуля матрица Aij разложения (10) имеет нулевые собственные значения. В то же время необходимо отметить, что задача магнитного перезамыкания является более сложной, чем классическая задача устойчивости векторного поля, зависящего от параметра. В нашем случае переход структурно-неустойчивой магнитной конфигурации в устойчивую определяется самосогласованным нелинейным взаимодействием двух векторных полей — магнитного поля В и поля скоростей V. Здесь много чисто теоретических вопросов, требующих дополнительных исследований.
В(ж, у, z) = 1,3хех + 0,7уеу - 2zez.
(18)
В(ж, у, z) = 0,25жеж — 0,75уеу + 0,5zez.
(19)
3.3. Визуализация. Система уравнений (11)—(15) содержит две векторные и три скалярные функции. По результатам ее решения нужно дополнительно рассчитать магнитное поле В и электрический ток ,]. Таким образом, на каждом временном шаге в машине формируется 15 трехиндексных массивов. Чтобы проанализировать и осмыслить такой огромный объем информации, необходимо использовать все средства, среди которых важное место занимает визуализация. В нашей работе было применено два метода: для изображения токовых слоев используются изоповерхности плотности тока, а для передачи структуры трехмерных трехкомпонентных векторных полей — их силовые линии. В то время как для изображения изоповерхностей предложено много стандартных методов, качественное представление пространственных кривых потребовало от нашей группы разработки новых программных средств.
Метод передачи трехмерности силовых линий за счет использования проекционных преобразований для создания объемного эффекта оказался неэффективным. Кроме того, перекрытие силовых линий на изображениях сравнительно невелико, поэтому удаление невидимых точек, а также изображение теней не дает заметного эффекта. В результате для передачи трехмерности изображения был выбран один из наиболее мощных механизмов — реалистичное освещение.
Классическая модель Фонга была разработана для моделирования освещенности поверхности точечным источником света. В этой модели предполагается, что в расчете интенсивности света каждой точки поверхности принимают участие три самостоятельные компоненты, две из которых предполагают задание нормали. Применение данной модели к одномерным многообразиям, которыми являются интегральные линии, затруднено невозможностью однозначного определения нормали. Для преодоления этой трудности предлагается использовать два подхода.
Первый из них основан на замене линии тонкой трубкой с помощью поверхности свертки. Поверхность свертки — это изоповерхность скалярного поля, полученного пространственной сверткой ядра полиномиального вида /г(г) вдоль изображаемой кривой с индикатором д(г):
Ф(г) = ¿г(г) ® /г(г) = J д(в)/г(г - в) я3
Данный тип поверхностей характеризуется бесконечной гладкостью, что существенно улучшает визуальное восприятие и необходимо для создания анимации. Кроме того, варьирование параметров ядра свертки позволяет отображать дополнительные параметры, в частности величину поля, кручение интегральной линии, различные дополнительные параметры, а также связанные скалярные поля. Существенным недостатком подхода является то, что получающиеся поверхности характеризуются довольно большим значением кривизны. Для качественного изображения таких поверхностей необходимо использовать сложные методы типа трассировки лучей, поскольку полигональная аппроксимация с последующей закраской полигонов по Гуро предполагает малые изменения нормали поверхности в рамках одного многоугольника, что неверно для поверхностей с большой кривизной.
Второй подход состоит в изменении модели освещенности. Новая модель основана на следующем замечании: если вместо кривой рассматривать бесконечно тонкую трубку и на ней всевозможные единичные нормали, то наибольшая освещенность будет соответствовать той нормали, которая лежит в плоскости, образованной касательным вектором и вектором к источнику света, и направлена в полупространство источника света. Таким образом, если предположить, что наблюдатель в силу малых размеров трубки видит лишь наиболее яркую ее точку, то можно определить новую модель, в которой вместо нормали к поверхности будет использоваться касательный вектор.
Такая модель прекрасно подходит для эффективной реализации с использованием текстур. Действительно, в предположении, что сцена освещается лишь одним источником света, освещенность каждой точки будет определяться двумя скалярными произведениями, которые можно интерпретировать как текстурные координаты, лежащие в диапазоне от —1 до 1. Таким образом, если предварительно рассчитать текстуру освещенности на некоторой сетке, то с помощью линейной интерполяции текстурных координат можно смоделировать линейную интерполяцию нормалей для каждого звена ломаной, аппроксимирующей изображаемую кривую. К недостаткам данного метода можно отнести невозможность передачи дополнительных параметров иными способами, кроме цветопередачи.
4. Математическое моделирование электрической активности головного мозга. В последнее время резко возрос интерес к математическому моделированию электрической активности головного мозга. Помимо медицинских и нейрофизиологических приложений этот интерес связан с
развитием нового направления в информатике — созданием интерфейса "мозг-компьютер", устанавливающего связи между ощущениями человека и их внешним проявлением, которое может быть зафиксировано приборами. В работах [37-40] методами математического моделирования исследовалась связь между источником электрического поля внутри человеческого мозга и результатами измерения потенциала электрического поля в процессе электроэнцефалографии (ЭЭГ).
Нейронные источники тока в мозге, генерирующие ЭЭГ-сигнал, могут быть разделены на две компоненты: первичные электролитические токи, которые выражают нейронные микроскопические пассивные клеточные токи, и вторичные или объемные токи, появляющиеся в результате макроскопического электрического поля. Первичные токи представляют основной интерес в электроэнцефалографии, так как они представляют области нейронной активности, ассоциированной с процессами чувств, моторным или мыслительным процессом.
В используемой модели мозг считается объемным проводником, в котором расположены генераторы стороннего электрического тока, представляемые в виде набора токовых диполей. Прямая задача ЭЭГ сводится к решению уравнения Пуассона в сложной области с сильно меняющимися коэффициентами электропроводности. Обратная задача ЭЭГ состоит в определении положений и моментов токовых диполей по измеренному потенциалу на скальпе.
Предложен численный подход к решению прямой задачи, основанный на представлении численного решения в виде разложения по сферическим гармоникам для произвольной формы области [37]. Коэффициенты разложения находятся из граничных условий и условий сшивки методом наименьших квадратов. Показано, что метод позволяет моделировать сложные многослойные структуры тканей мозга с анизотропной проводимостью с коэффициентами, меняющимися на несколько порядков.
Обратные задачи рассмотрены в двух постановках [38]. В первой постановке задачи предполагаются известными положения источников на некоторой разностной сети, а разыскиваются моменты диполей (линейная задача компьютерной ЭЭГ-томографии). Показано, что регуляризация задачи по Тихонову позволяет локализовать активные зоны внутри мозга, ответственные за создание регистрируемой ЭЭГ. Оценивается точность локализации и зависимость квазирешения от параметров регуляризации. Во второй постановке решается нелинейная задача локализации источников и находятся как положения, так и моменты эквивалентных диполей.
Рассмотрены как градиентные, так и стохастические методы решения задачи. Так как минимальное пространство параметров равно шести (один диполь), то для локализации нескольких диполей используется стохастический метод — генетический алгоритм, развитие которого применительно к таким задачам было выполнено ранее [41].
Разработана эволюционная модель обратной задачи ЭЭГ, основанная на использовании генетического алгоритма для поиска оптимального расположения диполей [39]. Предложенный метод позволил определять необходимое число диполей и их стабильные комбинации для описания измеренного пространственно-временного поведения электрического потенциала [39, 40].
Работа выполнена при финансовой поддержке АО СИМЕНС, контракт А-2003 с лабораторией нейроинформатики СИМЕНС (Мюнхен, Германия).
5. Эконофизика. Попытка применения математических методов в экономике выявила схожесть в описании многих экономических и физических явлений. Это дало толчок к развитию нового научного направления, которое получило название эконофизики и в рамках которого экономические процессы описываются с помощью математических моделей, заимствованных из физики.
5.1. Математические модели урбанистики. Проблемы городов весьма усложнились в результате технологического прогресса и изменения поведения людей. Городские системы характеризуются возрастанием пространственного и временного разнообразия протекающих в них процессов. Централизация городов наблюдается в развитых и развивающихся странах, но в крупных мегаполисах начали проявляться процессы децентрализации. Образцами таких городов являются Нью-Йорк, Париж, Токио, Москва. В экономической географии построено довольно много моделей образования и развития городов. В настоящее время преобладают три подхода: экономико-урбанистический, основанный на исследовании равновесных состояний; регионально-дискретный, представляющий некоторую разновидность клеточного автомата, и пространственно-динамический, рассматривающий проблему временной эволюции городов с учетом пространственной зависимости исследуемых величин. Последний подход стал развиваться относительно недавно. Его основу составляет обобщенное уравнение
непрерывности [42-44]
ад
— + <ИуЗ = Р(А), (21)
где А — плотность некоторой величины, J — плотность потока, а -Р(А) — интенсивность источников. В случае когда J = А\~, уравнение описывает конвективный перенос, в случае J = —кУА уравнение представляет собой хорошо известное уравнение диффузии.
В качестве примера рассмотрим результаты компьютерного моделирования рынка недвижимости в Москве. Соответствующее уравнение для цены жилья А имеет вид
<9А
— = <Цу(к(Х)£таЛ\)-г(\), (22)
где к — коэффициент диффузии, г(А) = А2(1 — А) — локальный избыток предложения. Подобный выбор функции г соответствует статистическим данным.
В качестве области исследования рассматривался правильный шестиугольник, который довольно хорошо описывает границы Москвы. На его границе ставились условия, соответствующие свободной миграции. При компьютерном моделировании использовалась пространственная сетка 300 X 300, в которой одна ячейка имела площадь 0,013 км2, а шаг по времени соответствовал двум месяцам.
В нашем примере городское пространство рассматривается как среда, в которой происходит диффузия [45]. Действительно, в Москве существуют престижные для застройки районы, например исторический центр. Однако ряд нормативных актов ограничивает там жилищное строительство. Наиболее "открытыми" для строительства жилья являются районы города, где есть пустоши или территории, которые освобождаются в результате реорганизации промышленных предприятий. Эти факторы и определили выбор коэффициента диффузии.
За начальное условие принималась ценовая функция первичного жилья г. Москвы на 2003 г. С учетом диффузионного члена строился прогноз развития ценовой политики на ближайшие годы. Результаты исследования такой модели показывали, что в ближайшие пять лет наличие жилья и связанная с ним плотность населения в центре города уменьшится, в то время как плотность застройки и плотность населения Юго-Запада, Севера и Востока Москвы увеличится.
5.2. Моделирование финансовых рынков. Для описания процессов в плазме широко используется кинетический подход, о котором мы коротко говорили в разд. 2.5. Рассмотрим теперь возможность его применения для моделирования флуктуации валютного рынка (FOREX). Исследование реальных баз тиковых данных с валютных рынков (данных о сделках на рынке, с указанием реального времени совершения сделки) показало, что курсовые разницы на разных временных шкалах удовлетворяют уравнению Чепмена-Колмогорова и, следовательно, система обладает свойством эргодичности. Такие системы могут быть описаны с помощью уравнения Фоккера-Планка
dp(x,t) _ d2A(x,t)p(x,t) dB(x,t)p(x,t)
dt ~ дх2 + дх ' 1 j
где p(x,t) — плотность вероятности курсовой разницы х в момент времени t. Уравнение (23) рассматривается на бесконечной прямой. Если мы знаем функцию p(x,t) в некоторый момент времени t = to и примем эту информацию за начальное условие, то, решая задачу Коши, можно спрогнозировать дальнейшее поведение этой функции при t > to.
Решить численно такую задачу несложно. Основная проблема заключается в определении коэффициентов уравнения A(x,t) и B(x,t). Их подбирают на основании анализа реальной базы данных сделок. Исследования показывают, что коэффициенты имеют значения одного порядка, т.е. оба одинаково важны для описания процесса эволюции функции p(x,t) со временем.
Анализ уравнения (23) позволяет сделать выводы о свойствах функции распределения курсовых разниц. Показано, что гипотеза о гауссовом распределении, корректность которой долго обсуждалась в научном сообществе, неверна. Численное моделирование позволило рассмотреть переходный процесс установления функции распределения курсовых разниц, а аналитическое исследование дало возможность изучить асимптотику поведения функции распределения. Распределение вероятностей
— а2
при Ai —т- 0 имеет главный член, равный |Дж| . Таким образом, функция распределения вблизи нуля будет иметь существенно более "толстые хвосты", что хорошо известно из реальных биржевых
данных. Результаты расчетов показали хорошее количественное совпадение с реальными рыночными данными [46].
Другим активно развиваемым направлением является исследование рынка опционов. С точки зрения биржевого брокера, очень важно знать "справедливую" цену опциона в каждый момент времени. На основании данных о рыночной и "справедливой" ценах брокер принимает решения о покупке или продаже. Для определения этой цены существует несколько методик, которые основываются на различных предположениях о характере рыночных процессов. Например, одна из наиболее известных методик, так называемая формула Блэка-Шолза, условием применимости имеет требование, чтобы случайный процесс колебания цены был винеровским. Анализируя реальные рыночные данные, легко убедиться, что практически все накладываемые ограничения, к сожалению, не выполняются в реальной жизни.
Исследование цены американских опционов в постановке проблемы, близкой к реальности, приводит к задачам, которые не имеют аналитического решения. Для их численного исследования необходимо имитационное моделирование (разыгрывание) случайного процесса флуктуации цены базового актива и многокритериальная минимизация риска для портфеля, содержащего базовый актив и опцион. С математической точки зрения задача формулируется следующим образом. Рассмотрим американский опцион на некий базовый актив со сроком исполнения N дней и ценой исполнения С. В начальный момент времени известна цена единицы базового актива. Пусть задано распределение, которому подчиняется колебание цены базового актива. Задача состоит в том, чтобы построить инвестиционный портфель, который имел бы минимальный риск в каждый из дней от 1 до N при условии, что колебание базового актива подчиняется заданному распределению. В общем случае данный портфель будет иметь каждый день различное количество единиц базового актива на один опцион. Зная портфель с минимальным риском, можно определить справедливую цену опциона как сумму риска и стоимости всех единиц базового актива, входящих в портфель. Таким образом, задача состоит в минимизации функционала риска R(<p), где <р = (ípi, ipi,..., <^дг) > Vfc G ( —оо,+оо) — количество единиц базового актива, входящих в портфель. Отрицательное значение соответствует короткой позиции. Единица базового актива предполагается дробящейся. Такая постановка задачи соответствует бесконечной ликвидности базового актива, что не всегда соответствует реальным условиям. В случае ограниченной ликвидности необходимо рассматривать задачу на некотором ограниченном множестве ф G Ф.
Очевидно, что цена опциона будет меняться в зависимости от меры риска, т. е. от того, каков будет вид функционала R. В качестве наиболее распространенных мер риска обычно используется дисперсия совокупной цены портфеля, вероятность изменения цены портфеля больше, чем некая критическая величина (VAR), и тому подобное. На практике для определения цены опциона берется взвешенное среднее значение по нескольким мерам риска.
Сложность задачи усугубляется очень жесткими временными ограничениями, накладываемыми на скорость получения решения, так как для многих реальных ситуаций на бирже время реакции системы не должно превышать секунды, а для некоторых приложений (статистический арбитраж) должно измеряться сотыми долями секунды. Данный подход требует разработки специальных кластерных решений и соответствующих алгоритмов для достижения заданных временных параметров.
СПИСОК ЛИТЕРАТУРЫ
1. Днестровский Ю.Н., Костомаров Д. П. Математическое моделирование плазмы. М.: Наука, 1982. (Издание на англ. языке: Dnestrovskij Yu. N., Kostomarov D.P. Numerical simulations of plasmas. N.Y.: Springer-Verlag, 1986.)
2. Popov A.M.,La Haye R. J., M u r ak am i M.,Liu Y. Q., P о p о v a N. N., T u r n b u 11 A.D.Simulation of neoclassical tearing modes (NTMs) in the DIII-D tokamak. Part I. NTM Excitation // Physics of Plasmas. 2002. 9. N 10. P. 4205-4228.
3. Popov A.M., La Haye R. J., M u r ak am i M., Liu Y. Q., P o p o v a N. N., T u r n b u 11 A.D. Simulation of neoclassical tearing modes (NTMs) in the DIII-D tokamak. Part II. Suppression by radially localized electron cyclotron current drive // Physics of Plasmas. 2002. 9. N 10. P. 4229-4240.
4. Popov A.M., Chan V.S., Chu M.S., Liu Y. Q., Turn bull A.D. Nonlinear three-dimensional self-consistent simulations of negative central shear discharges in the DIII-D tokamak // Physics of Plasmas. 2001. 8. N 8. P. 3605-3619.
5. Popov A.M., Liu Y.Q., Popova N.N., La Haye R. J., T u r n b u 11 A. D. NFTC code for simulation of nonlinear instabilities in toroidal Plasmas // International Conference on Computatiolal Physics. Granada, Spain. 1998. P. 200.
6. Popov A.M. Simulation of neoclassical tearing modes coupling in tokamak discharges using 3D nonlinear code NFTC // Bull. International Sherwood Fusion Theory Meeting. Rochester; N.Y., 2002. P. 22-24.
7. Popov A.M., La Haye R.J. Simulations of ECCD suppression of neoclassical tearing modes in sawteething discharges using 3D nonlinear code NFTC // Bull. 28th EPS Conference on Plasma Physics. Madeira, 2001. P. 329.
8. Popov A.M., Pustovitov V.D. Simulation of neoclassical tearing Modes in ITER Using 3D Nonlinear code NFTC // 30th EPS Conf. on Contr/Fusion and Plasma Phys. St. Petersburg. 7-11 July 2003. ECA, 27A. P. 3-140.
9. Popov A.M. 3D simulations of neoclassical tearing modes in ITER // 31st EPS Conf. on Plasma Physics. 2004. London: ECA, 28A. P. 4-144.
Ю.Бондарчук Э.Н., Вознесенский В. А., Д н e с т p о в с к и й Ю.Н., Леонов В.М., Максимова И. И., Сычугов Д.Ю., Цаун С. В. Вертикальная МГД-устойчивость плазмы в Т-15М // 7-я Международная конференция по инженерным проблемам термоядерных реакторов. Санкт-Петербург. 28-31 октября 2002. ИПТР, 2002. С. 197.
11. В о n d аг с h u k E.N., Dnestrovskii Yu.N., Leonov V.M., Maksimova 1.1., Sychu-gov D.Yu., Tsaun S.V., Voznesensky V.A. Vertical MHD stability of the T-15M tokamak plasma // Plasma Devices and Operations. 2003. 11. N 4. P. 219-227.
12. Wilson H.R., Voss G.M., Akers R.J., Appel L., С h r i s t i an s e n J. P., D n e s t г о v s к i j A., Keating O., Hender T.C., Hole M.J., Huysmans G., Kirk A., Knight P.J., Louglin M., McClemens K.G., O'Brien M.R.,Sychugov D. Yu., V al о v i с M. The physics basis of spherical tokamak components test facility // 31st EPS Conf. on Plasma Physics. 2004. London: ECA, 28 June-2 July. 28A. P. 4-196.
13. Zaitsev F.S., Shishkin A.G., Kostomarov D.P., O'Brien M.R., Akers R. J., Gryazne-vi с h M., T r ef i 1 о v А. В., Yelchaninov A.S. The numerical solution of the self-consistent evolution of plasma equilibria // Сотр. Phys. Comm. 2004. 157/2. P. 107-120.
14. Kostomarov D. P., Zaitsev F.S., Shishkin A.G., O'Brien M.R., Gryaznevich M. Access to "advanced" regimes in tight aspect ratio plasmas // 25th European Conference on Controlled Fusion and Plasma Physics. Praha, 1998. P. 4.
15. Kostomarov D.P., Zaitsev F.S., Shishkin A.G., O'Brien M.R., Gryaznevich M., Akers R. J., Krastylev A.V. Access to optimised shear equilibria in spherical tokamaks // 26th European Conference on Controlled Fusion and Plasma Physics. Maastricht: ECA, 1999. 23J. P. 17451748.
16. Akers R.J., Bond A., Buttery R.J., Carolan P.G., Counsell G.F., Cunningham G., Fielding S.J.,Gimblett C.G.,Gryaznevich M.,Hastie R. J.,Helander P.,Hender T.C., Knight P., Lashmore-Davies C. N., M ad d i s о n G.P., Martin T.J., McClements K.G., Morris A.W., O'Brien M.R., Ribeiro C., Roach C.M., Robinson D.C., Sykes A., Voss G.M., Walsh M.J., Wilson H.R., Zaitsev F.S. Steady state operation of spherical tokamaks // Nucl. Fusion. 2000. 40. N 6. P. 1223-1244.
17. Костомаров Д. П., Зайцев Ф.С., Akers R. J., Шишкин А.Г. Исследование электрической проводимости плазмы в сферическом токамаке // Докл. РАН. 2004. 396. № 6. С. 762-765.
18. Зайцев Ф. С., Лукьяница А. А., Шишкин А. Г. Применение нейронных сетей в задачах математического моделирования тороидальной плазмы // Докл. РАН. 2003. 393. № 2. С. 173-175.
19. Forest С. В., Harvey R. W., Smirnov А. P. Power deposition by mode converted electron Bernstein waves in the DIII-D "heat-pinch" experiments // Nuclear Fusion. 2001. 41(5). P. 619-623.
20. Forest C.B., С h a 11 о p ad h у ay P.K., Harvey R.W., Smirnov A. P. Off-midplane launch of electron Bernstain waves for current drive in overdense plasmas // Physics of Plasmas. 2000. N 7. P. 1352-1355.
21. Chattopadhyay P. K., Anderson J.K.,Biewer Т., Craig D., Forest С. В., Harvey R. W., Smirnov A. P. Electron Bernstein wave emission from an overdensed reversed field pinch plasma // Physics of Plasmas. 2002. 9. N 3. P. 752-755.
22. Taylor G., Efthimion Р.С., Jones D., Le Blanck B.P., Wilson J.R., Wilgen J.B., Bell G.L., Bigelow T.S., Maingi R., Rasmussen D.A., Harvey R.W., Smirnov A.P., Paoletti F., Sabbagh S.A. Enhanced conversion of thermal electron Bernstein waves to the extraordinary electromagnetic mode on the National Spherical Torus Experiment (NSTX) // Physics of Plasmas. 2003. 10. N 10. P. 1395-1401.
23. Bernabei S., Bundy R.V., Fredrickson E.D., Gorelenkov N.N., Hosea J.C., Phillips C.K., White R.W., Wilson J.R., Petty C.C., Pinsker R.I., Harvey R.W., Smirnov A. P. The combined effect of EPM and TAEs on energetic ion confinement and sawtooth stabilization // Nuclear Fusion. 2001. 41(5) P. 513-518.
24. Zaitsev F.S., O'Brien M.R., Cox M. Three dimensional neoclassical non-linear kinetic equation for low collisionality axisymmetric tokamak plasmas // Physics of Fluids B. 1993. 5. N 2. P. 509-519.
25. Костомаров Д. П., Зайцев Ф. С. Фундаментальные свойства оператора кулоновских столкновений в усредненном кинетическом уравнении // Докл. РАН. 2002. 384. № 6. С. 747-750.
26. Zaitsev F.S., Longinov V.V., O'Brien M.R., Tanner R. Difference schemes for the time evolution of three-dimensional kinetic equations // J. Comput. Phys. 1998. 147. P. 239-264.
27. Зайцев Ф.С., Костомаров Д. П., Курбет И. И. Применение явных итерационных схем для решения кинетических задач // Матем. моделир. 2004. 16. № 3. С. 13-21.
28. O'Brien M.R., Сох М., Gardner С.A., Zaitsev F.S. 3D Fokker-Planck calculation of alpha-particle distributions: a TFTR simulation // Nucl. Fusion. 1995. 35. N 12. P. 1537-1547.
29. Zaitsev F.S., Akers R. J., O'Brien M. R. Perturbations to deuterium and tritium distributions caused by close collisions with high-energy alpha-particles // Nucl. Fusion. 2002. 42. P. 1340-1347.
30. Bulanov S.V., Echkina E.Yu., Inovenkov I.N., Pegoraro F., Pichushkin V. V. On the structural stability of magnetic configurations with two null lines // Physics of Plasmas. 1999. 6. P. 802.
31. Kostomarov D. P., Echkina E.Y., Inovenkov I.N., Leonenko A.V., Pavlova O.A. Application of parallel progromming method for the 3D MHD computer simulations of magnetic reconnection in plasma // Proc. of Parallel Computational Fluid Dynamics. ELSERVIER. 2004. P. 373-379.
32. Буланов С.В., Ечкина Е.Ю., Иновенков И.Н., Пегораро Ф., Пичушкин В.В. Формирование токовых слоев в структурно-устойчивых и структурно-неустойчивых магнитных конфигурациях с двумя нулевыми линиями // Физика плазмы. 2000. 26. № 7. С. 600-614.
33. Bulanov S.V., Echkina E.Yu., Inovenkov I.N., Pegoraro F. Current sheet formation in three-dimensional magnetic configuration // Physics of Plasma. 2002. 9. N 9. P. 3835-3850.
34. Костомаров Д.П., Буланов С.В., Ечкина Е.Ю., Иновенков И.Н., Леоненко А.В., Пичушкин В. В., Пегораро Ф. Структурно-неустойчивые магнитные конфигурации в трехмерной геометрии // Докл. РАН. 2003. 390. № 2. С. 178-182.
35. Bulanov S. V., Inovenkov I.N., Echkina E.Y., Leonenko A. V., Pegoraro F. On the form of the MHD singularity in the 3D and 2D structurally unstable magnetic configurations // 30th EPS Conference on Controlled Fusion and Plasma Physics. 2003. 27A. P. 1-36.
36. Буланов С.В., Ечкина Е.Ю., Иновенков И.Н., Пегораро Ф., Пичушкин В. В. Трехмерное МГД-моделирование вынужденного магнитного перезамыкания // Физика плазмы. 2001. 27. № 4. С. 334-342.
37. Хоффманн К., Попов A.M., Федулова И.А., Певцов С. Е. Численное решение прямых математических задач электроэнцефалографии // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2004. № 1. С. 19-26.
38. Хоффманн К., Попов A.M., Федулова И. А., Певцов С.Е. Численное решение обратных математических задач электроэнцефалографии // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2004. № 4. С. 16-27.
39. Хоффманн К., Попов А.М., Федулова И.А., Певцов С.Е. Моделирование пространственно-временной электрической активности нейронных источников // Прикл. матем. и информ. № 17. М., 2004. С. 55-71.
40. Hoffmann С., Popov A.M., Fedulova I. A., Pevtsov С.Е. EEG data interpretation method for determination of possible underlying source. Submited to IEEE Trans on Biomed. Eng. 2004.
41. Малютина Э. Э. Применение стохастических методов для расчета управления динамической системой // Информационные технологии и вычислительные системы. 2000. 1. С. 52.
42. Kostomarov D.P., Echkina E.Yu., Inovenkov I.N., Bulanov S.V., Leonenko A.V. The nonlinear dynamics of the business center in Beckmann's model // Abstracts in International Conference on Applications of Physics in Financial Analysis 4. 2003. 27C. P. A67.
43. Bulanov S. V., Echkina E. Yu., Inovenkov I. N. The evolution of the spatial economic model / / Abstracts Europhysics Conference on Computational Physics. Germany, Aachen, 2001. P. A115.
44. Bulanov S. V., Echkina E. Yu., Inovenkov I.N. The transformation of structurally unstable configurations into structurally stable ones in the Beckmann model // Abstracts in International Conference on Applications of Physics in Financial Analysis 1. 2000. 25C.
45. Bulanov S. V., Echkina E.Yu., Inovenkov I.N. The nonlinear dynamics of the business center in Beckmann's model // Physics Letters A. 2004. 344. P. 104-107.
46. Smirnov A.P., Shmelev А.В., Sheinin E. Y. Analysis of Fokker-Planck approach for foreign exchange market study // Abstracts in International Conference on Applications of Physics in Financial Analysis 4. 2003. 27C. P. A85.
Поступила в редакцию 15.09.04
Е. И. Моисеев, И. С. Ломов
СПЕКТРАЛЬНЫЕ РАЗЛОЖЕНИЯ И ГРАНИЧНОЕ УПРАВЛЕНИЕ В ЗАДАЧАХ МАТЕМАТИЧЕСКОЙ ФИЗИКИ
(кафедра общей математики факультета ВМиК)
В данной заметке изложены основные направления научных исследований преподавателей кафедры общей математики. На кафедре проводятся исследования по ряду фундаментальных проблем математической физики, теории оптимального управления, теории функций, вычислительных методов линейной алгебры, дифференциальной геометрии и компьютерной алгебры, актуарной математики. В частности, изучается общая задача граничного управления колебательными процессами, в том числе и задача о полном успокоении колебательного процесса. Исследуются спектральные характеристики различных задач математической физики. Эти задачи возникают при расчете ядерных реакторов, при исследовании устойчивости плазмы, колебательных процессов, в том числе и процессов, протекающих с диссипацией, т.е. с затуханием, а также при изучении сверхзвуковых и дозвуковых течений газа. Рассматривается устойчивость в процессах с нелинейными колебаниями, возникающими, например, при вращении гироскопов, движении частиц в ускорителях, изучаются вопросы о глобальном во времени существовании и об асимптотическом поведении решений нелинейных эволюционных уравнений и систем уравнений. Развиваются как аналитические, так и численные методы исследования дифференциальных уравнений в частных производных. Изучаются линейные уравнения смешанного типа, системы уравнений типа Стокса. Введены в расмотрение и успешно применяются для решения ряда задач новые функциональные пространства — усиленные и ослабленные пространства Соболева.
Кафедра общей математики создана в 1973 г. по инициативе академика А. Н. Тихонова. Организатор и бессменный руководитель кафедры — лауреат Государственной премии академик РАН В. А. Ильин. Владимир Александрович Ильин создал широкий круг единомышленников, развивающих математическую науку в России и за ее пределами. Более половины сотрудников кафедры принадлежат к большой научной школе академика Ильина.
На кафедре работают декан факультета ВМиК МГУ академик РАН Е. И. Моисеев, профессора В. В. Власов, A.A. Дезин, Е.Г. Дьяконов, Х.Д. Икрамов, И. С. Ломов, М.М. Хапаев, Е.В. Шикин, И. А. Шишмарев, 14 доцентов, 1 старший преподаватель, 3 ассистента.
Опишем ряд конкретных полученных результатов.
1. Вопросы спектральной теории самосопряженных дифференциальных операторов.
В конце 1960-х—начале 1970-х гг. В. А. Ильиным был разработан универсальный метод изучения спектральных разложений, отвечающих произвольным самосопряженным расширениям эллиптических операторов [1]. Хотя вопросам локализации и равномерной сходимости спектральных разложений
