автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.16, диссертация на тему:Использование систем компьютерной алгебры при построении эффективных численных алгоритмов решения некоторых задач математической физики

ОБЪЕДИНЁННОЙ ИНСТИТУТ ЯДЕРНЫХ ИССЛЕДОВАНИЙ ЛАБОРАТОРИЯ ВЫЧИСЛИТЕЛЬНОЙ ТЕХНИКИ И АВТОМАТИЗАЦИИ

На правах рукописи УДК 519.6

БОГОЛЮБСКАЯ Алла Анатольевна

ИСПОЛЬЗОВАНИЕ СИСТЕМ КОМПЬЮТЕРНОЙ АЛГЕБРЫ ПРИ ПОСТРОЕНИИ ЭФФЕКТИВНЫХ ЧИСЛЕННЫХ АЛГОРИТМОВ РЕШЕНИЯ НЕКОТОРЫХ ЗАДАЧ МАТЕМАТИЧЕСКОЙ ФИЗИКИ

05.13.16 - применение вычислительной техники, математического моделирования и математических методов для научных

исследований

Диссертация на соискание учёной степени кандидата физико-математических наук

Научный руководитель -

доктор физико-математических наук

СЕРДЮКОВА С.И.

Дубна - 1998

Оглавление

Введение .............................................................1

1 Построение явных С-устойчивых схем максимального нечётного порядка точности с применением системы REDUCE 16

1.1 Постановка задачи. Алгоритм построения схем порядка

q = 2k-l............................. 16

1.2 Проверка точности и устойчивости построенных разностных схем. Особенности использования системы REDUCE

в данной задаче ......................... 19

1.3 Схемы максимального нечётного (3,5,7,9) порядка точности 21

1.3.1 Случай к = 2....................... 21

1.3.2 Случай к = 3 ...................... 23

1.3.3 Случай к = 4....................... 24

1.3.4 Случай к = 5....................... 26

1.4 Случай схемы максимального чётного порядка точности . 27

2 Исследование устойчивости разностных краевых задач с

применением систем компьютерной алгебры 34

2.1 Необходимые сведения из GKS-теории ........................34

2.2 К исследованию спектра одной разностной краевой задачи 37

2.3 Возможность полного исследования устойчивости разностных краевых задач с применением СКА ......................51

2.3.1 Примеры полного исследования устойчивости ... 58

3 Математическое моделирование двумерных солитонов в

моделях магнетиков Гейзенберга (МГ) 61

3.1 Двумерные топологические солитоны в модели анизотропного МГ со стабилизирующими членами...........61

3.2 Двумерные топологические солитоны в калибровочно-инвариантной модели легкоосного антиферромагнетика

Гейзенберга ........................... 71

Заключение ........................................................84

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

Введение

Численные исследования и компьютерная алгебра в настоящее время одинаково успешно используются в цикле вычислительного эксперимента [1, 2]. Численные методы, развиваясь по мере развития теории и совершенствования ЭВМ, остаются важнейшим, а иногда и единственным средством решения сложных задач математической физики, например, нелинейных уравнений. В связи с этим возрастает роль теоретического обоснования корректности численных методов. Для этой цели, а также для построения численных методов с заранее заданными свойствами всё чаще используются системы компьютерной алгебры (СКА) [3]-[5].

Построение разностных аппроксимаций заданного качества — „один из главных вопросов теории численных методов," — отмечает в [1] академик А.А.Самарский. Построение разностных схем связано с громоздкими выкладками и поэтому является естественной областью для приложений систем компьютерной алгебры (СКА). Перечислим наиболее известных авторов ([б]-[16]): Вирт М.С., Ворожцов Е.В., Ганжа В.Г., Лиска Р., Мазурик С.И., Стейнберг С., Шапеев В.П., Шашков М.Ю., Щенков И.Б. (Отметим, что почти все эти авторы занимаются также исследованием устойчивости разностных схем с применением СКА). В приведённых работах для построения разностных схем используются различные методы - метод конечных элементов [7], метод разностных аппроксимаций [11], метод неопределённых коэффициентов [9] и др. Иногда авторы разрабатывают специальные методы и языки [6, 8], удобные для заданных классов дифференциальных уравнений. Программы, связанные с построением и исследованием устойчивости разностных схем, включаются в библиотечные пакеты универсальных СКА (см., например, [17]).

Построенные разностные аппроксимации используются, как правило,

для проведения дальнейших численных экспериментов, поэтому возникает необходимость представления их в формате языков численного программирования. Одним из важнейших требований к СКА при построении разностных схем (и вообще при построении эффективных численных алгоритмов) является хорошо организованный символьно-численный интерфейс [18], и, в частности, возможность генерации эффективных программ на языках численного программирования (а в идеале — и для современных параллельных суперкомпьютеров).

Одним из первых, кто в нашей стране применил СКА для построении разностных схем, был В.В.Русанов — еще в 1969 г. в [23] он использовал язык РЕФАЛ [24] для вычисления коэффициентов своей знаменитой явной схемы, широко используемой в газодинамических расчётах. Схема Русанова имеет третий порядок точности. С увеличением порядка точности необходимость использования СКА возрастает.

К схемам высокого порядка точности относятся так называемые схемы максимального нечётного порядка точности. Эти схемы были предложены в работе [26] Стрэнгом (для уравнения щ ~ их) и рассматривались Цинь Мэн-чжао в работе [27], где изучались конечноразностные методы решения задачи Коши для уравнения щ + и^ = 0, t = пт, х — иЬ. Накладывалось естественное ограничение на отношение шагов сетки (т/Н) < 1. Было показано, что [27] на схемах порядка точности 1пД-1 достигаются оптимальные в Ь2 оценки сходимости.

С.И.Сердюковой в работах [28]-[30] рассматривался подкласс С-устой-чивых схем максимального нечётного порядка точности (2к — 1). При начальных данных из С£ для схем порядка точности 1п/г-1 установлена оценка погрешности решения в С порядка точности 1п 1п/г-1). Кроме того, доказано, что число точек, на которые „размывается" изолированный разрыв, обратно пропорционально порядку схемы. Зона размывания разрыва имеет ширину порядка 1п /г-1.

Тем самым было доказано, что схемы максимального нечётного порядка точности 2к — 1, к = 0( 1п/г-1), обладают свойствами, близкими к оптимальным в С, и хорошо приспособлены для счёта разрывных решений, которые нередко возникают в практических расчётах.

Построение аналога схем Стрэнга для двумерного случая производит-

ся в диссертации. Эту задачу для схем высоких порядков точности невозможно осуществить без применения СКА. Отметим, что рассматриваемое нами простейшее гиперболическое уравнение щ = сих обычно применяется при исследовании численных методов, разрабатываемых, например, для решения задач газовой динамики. Система уравнений газовой динамики в эйлеровых координатах

дю х, 2)

т

дх

/

т =

V

г

\1/

—г

\

-(7-1)/ +

(7 ~ 3> 2 р

- 7/

(7 — 1)г2\ г

2 р

Р

Р

+

7 — термодинамическая константа, р

(где г = ру, / = - (7_1} , 2 , плотность, V — скорость, р - давление) является квазилинейной гиперболической системой. Это означает, что в рассматриваемой области изменения матрица Тш имеет различные вещественные собственные значения. При фиксированных значениях переменных такая матрица приводится к диагональному виду, то есть, получается система

дщ дщ г дх'

Поэтому простейшее гиперболическое уравнение щ = сих является хорошей моделью для исследования явлений неустойчивости, которые наблюдаются в газодинамических расчётах. Таким исследованиям посвящены работы С.К.Годунова и В.С.Рябенького [32, 33], В.В.Русанова [23], А.А.Самарского [34, 39], Р.П.Федоренко [35], И.Н.Яненко [36], а также работы Х.-О.Крайса [54, 55], Дж.Олигера [59], Р.Рихтмайера и К.Морто-на [37], Г.Стрэнга [26], В.Томэ [38].

Анализ устойчивости РКЗ является ключевым вопросом теории численных методов. Для исследования устойчивости РКЗ применяются различные методы - метод энергетических неравенств [1, 39, 40]; метод разделения переменных [41, 44]; методы, основанные на принципе максимума [1, 41, 44] и др. Однако применение этих методов на практике

2

часто оказывалось сложным - во многом из-за необходимости проведения громоздких алгебраических преобразований. Поэтому при исследовании устойчивости часто приходилось „обращаться к эмпирике, к здравому смыслу, к интуиции" [45].

Появление и развитие СКА стимулировало появление работ, в которых устойчивость исследуется с помощью компьютера ( [46]-[53]) с применением, например, спектральных методов в [46, 47], метода дифференциальных приближений в [48] и других. Отметим, что в перечисленных работах по исследованию устойчивости речь идет в основном о задаче Коши. Математический аппарат исследования краевых задач, о которых пойдет речь далее, на порядок сложнее.

В диссертации рассматриваются линейные разностные краевые задачи (РКЗ), аппроксимирующие системы гиперболических [83] дифференциальных уравнений (ДУ). Гиперболические ДУ описывают широкий класс нестационарных процессов и часто встречаются в приложениях. Общая теория устойчивости для РКЗ этого класса изложена в работах [39]-[44], [54]-[58]. Получены необходимые и достаточные условия устойчивости полубесконечных разностных линейных РКЗ.

Теоретической основой проведённых нами исследований является GKS -теория, предложенная шведскими математиками Б.Густаффсоном, Х.-О.Крайсом и А.Сандстрёмом и в работах [54, 55] и развитая С.И.Сер-дюковой в работах [56]-[58]. Основные сведения из GKS-теории представлены в разделе 2.1 диссертации. GKS-теория опирается на важнейшие результаты спектральной теории Годунова-Рябенького [42]: для устойчивости полу бесконечной РКЗ необходимо, чтобы была устойчива соответствующая задача Коши и чтобы спектр оператора перехода от слоя к слою лежал в единичном круге \г\ < 1. Если есть точки спектра вне единичного круга, наблюдается сильная неустойчивость экспоненциального типа. Если спектр расположен сторого внутри единичного круга \z\ < 1, задача устойчива. Случай точек спектра на единичной окружности исследован С.И.Сердюковой в работах [56]-[58]. В этом случае задача может быть как устойчивой, так и неустойчивой.

Метод исследования устойчивости РКЗ, основанный на GKS-теории, носит название NMA — "Normal Mode Analysis". Анализ устойчивости

в NMA сводится к задаче определения точек спектра оператора G (оператора перехода от слоя к слою) и, тем самым, к алгебраической задаче: выводу и решению системы полиномиальных уравнений. Ряд исследователей решает эту задачу численно; рассмотрим кратко некоторые работы.

В работе [59] Дж.Олигер исследует устойчивость РКЗ для схемы 4-го порядка точности по пространственной переменной и 2-го порядка точности по времени для уравнения щ = —сих. В качестве дополнительных граничных условий используется аппроксимация 3-го порядка точности по пространственной переменной. Система нелинейных уравнений решается методом редукции - сведением к одному полиномиальному уравнению (очень высокой степени) от одного неизвестного. Для этого производится последовательное исключение переменных и многократное вычисление результанта (методом Коллинза вычисления полиномиальных результантов от многих переменных [60], [61]). Корни полинома находятся затем стандартными численными методами. Методика требует больших машинных ресурсов.

В работе [62] Дж.Гари обобщил предложенную Олигером аппроксимацию введением скалярного параметра в дополнительные граничные условия, что позволило управлять устойчивостью.

Д.М.Слоун [63] определил область устойчивости рассматриваемой РКЗ. Система полиномиальных уравнений в данной работе решается численно (методом, описанным в [60]), для фиксированных значений параметра а = т/h - отношения шагов сетки. Аналогичная методика применена в работе [64].

Наиболее эффективной из существующих численных методик исследования устойчивости с использованием GKS — теории представляется методика М.Сунэ [65, 66]. Во-первых, этот алгоритм предназначен не только для одномерного, но и для двумерного по пространственным переменным гиперболического уравнения, во-вторых, он специализирован: определяются не все решения системы полиномиальных уравнений, а только те, которые лежат в кольце 1 < \z\ < 1 + где е = 1/3A", N— число переменных в аппроксимируемой системе ДУ. (Необходимость нахождения только таких решений доказывается в [66]). Методика реализована

в виде системы автоматического исследования устойчивости гиперболических РКЗ определённого класса — ШвТАВ [65].

В этих работах все вычисления проводятся, как правило, для фиксированных значений параметра а (а - отношение шагов сетки), а это не гарантирует, что в области \г\ > 1 будут найдены все точки спектра РКЗ. В случае численного исследования окончательный вывод об устойчивости данной РКЗ сделать невозможно, поэтому в особенности важным представляется аналитическое исследование.

Оказалось, однако, что проведение такого исследования аналитически до конца требует значительных усилий даже для простых модельных задач и аппроксимаций невысокого порядка [76, 82], так как только вывод основной системы уравнений СКБ-теории (не говоря о её решении) является достаточно трудоёмким и связан с громоздкими аналитическими выкладками, выполнить которые без применения СКА не представляется возможным. В настоящей диссертации предложена методика исследования устойчивости РКЗ с помощью СКА, включающая вывод системы полиномиальных уравнений, описывающей спектр данной задачи, и решение этой системы с помощью специального алгоритма.

Разработка этой методики началась с исследования явления неустойчивости, которое наблюдалось при численном моделировании [68]- [70] движения флюксонов в протяжённых структурах с микронеоднородно-стями [67]. Использовалась явная схема Русанова [71] третьего порядка точности, имеющая 5 точек в основании и требующая дополнительных граничных условий. В зависимости от способа аппроксимации дополнительных граничных условий рассматриваемая РКЗ оказывалась устойчивой или неустойчивой (при счёте наблюдалась сильная неустойчивость экспоненциального типа, распространяющаяся от границ). Хотелось бы такого рода явления неустойчивости исключить заранее.

Как обычно при исследовании устойчивости, использовалось модельное уравнение - как правило, если удаётся выяснить причины возникновения неустойчивости на подходящем модельном уравнении, то удаётся найти устойчивый алгоритм решения задачи. В диссертации проводится теоретическое исследование граничных условий двух типов (определяющих как устойчивость, так и неустойчивость). Выводится система

полиномиальных уравнений, описывающая спектр задачи. При решении этой системы мы столкнулись со следующими трудностями. Попытки применить уже обсуждавшуюся программу IBSTAB [65], а также метод базисов Грёбнера [72, 73] успеха не имели. Применение формул Кардано и Феррари [74], а также общих символьных алгоритмов, реализованных в системе REDUCE, приводило к громоздким выражениям с радикалами, не поддающимся преобразованиям. С.И.Сердюковой и Н.Е.Мазепой был разработан специальный алгоритм [75]-[77] для решения систем полиномиальных уравнений, описывающих спектры РКЗ. Сложная алгебраическая задача вычисления спектра была сведена к решению одного полиномиального уравнения. Этот алгоритм был положен в основу программы SPECTR [79] для вычисления спектров разностных краевых задач. С помощью этой программы была исследована устойчивость ряда известных модельных краевых задач [75, 77], в частности, устойчивость дополнительных граничных условий для разностных схем максимального нечётного порядка точности [80].

В дальнейшем „алгоритм сведения" был усовершенствован [81]. Усовершенствованный вариант алгоритма включён в методику полного исследования устойчивости линейных разностных гиперболических краевых задач, которая обсуждается в диссертации. По сравнению с предыдущей новая версия алгоритма более автоматизирована; все этапы исследования проводятся на компьютере с помощью систем REDUCE и MAPLE [22]. Расширен круг исследуемых задач вплоть до схем 6 порядка точности. Алгоритм может быть использован для решения более сложных задач - так, в работе [82] данный алгоритм используется С.И.Сердюковой и М.Сунэ для исследования устойчивости схемы Русанова и схемы Гари с начальными условиями, дополнительными граничными условиями и условиями стыковки на областях составной структуры; используются СКА REDUCE и MAPLE. В случае схемы Гари была получена система уравнений седьмого порядка, которая была сведена к полиному 900 степени. Этот полином был вычислен и факторизован с помощью системы MAPLE; полученные уравнения (от одной переменной), порядки которых не превосходили 100, были решены численно с помощью REDUCE.

Усовершенствованный алгоритм „работает" следующим образом: сначала исследуемая нелинейная система уравнений, описывающая спектр, алгебраическими методами сводится к двум полиномиальным уравнениям. После вычисления результанта [72] получаем одно полиномиальное уравнение, которое решается численно. (Отметим, что на вычисление полиномиального результанта опираются многие современные алгоритмы компьютерной алгебры для решения нелинейных систем уравнений [72, 86, 87], развиваемые наряду с алгоритмами, основанными на использовании базисов Грёбнера [5, 72]).

Таким образом, рассма�