авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

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


Pages:   || 2 | 3 | 4 | 5 |   ...   | 9 |
-- [ Страница 1 ] --

Е.А. Новиков, Ю.В. Шорников

КОМПЬЮТЕРНОЕ

МОДЕЛИРОВАНИЕ

ЖЕСТКИХ ГИБРИДНЫХ

СИСТЕМ

Е.А. Новиков, Ю.В. Шорников

КОМПЬЮТЕРНОЕ

МОДЕЛИРОВАНИЕ

ЖЕСТКИХ ГИБРИДНЫХ

СИСТЕМ

НОВОСИБИРСК

2012

УДК 004.9

Н 731

Рецензенты:

Заслуженный деятель науки РФ,

д-р техн. наук, профессор В.И. Денисов;

д-р физ.-мат. наук, гл. науч. сотр. ИВТ СО РАН Л.Б. Чубаров Утверждено к печати Редакционно-издательским советом Новосибирского государственного технического университета и Научно-издательским советом СО РАН Новиков Е.А.

Н 731 Компьютерное моделирование жестких гибридных систем :

монография / Е.А. Новиков, Ю.В. Шорников. – Новосибирск :

Изд-во НГТУ, 2012. – 451 с.

ISBN 978-5-7782-2023- Монография посвящена проблеме построения оригинальных численных методов решения задачи Коши для жестких систем обыкновенных дифферен циальных уравнений. Особое внимание уделяется контролю точности вычисле ний и устойчивости численной схемы, а также созданию алгоритмов интегри рования переменного порядка и шага. Подробно рассматривается методология гибридных систем и приведена их классификация. Описаны возможности ин струментальной среды машинного анализа гибридных моделей. На ряде прак тических задач продемонстрированы особенности использования разработан ного программного комплекса.

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

УДК 004. ISBN 978-5-7782-2023-5 © Новиков Е.А., Шорников Ю.В., © Институт вычислительного моделирования СО РАН, © Новосибирский государственный технический университет, E.A. Novikov, Yu.V. Shornikov COMPUTER SIMULATION OF STIFF HYBRID SYSTEMS Monograph NOVOSIBIRSK UDC 004. N Reviewers:

Prof. V.I. Denisov, D. Sc. (Eng.), RF Honoured Scientists;

L.B. Chubarov, D. Sc. (Phys. & Math.), Chief Scientist, SB RAS Computing Technology Institute Approved for publishing by Editorial publishing Council Of Novosibirsk State Technical University (NSTU) and Scientific publishing Council of SB RAS Novikov E.A.

N 731 Computer Simulation of Stiff Hybrid Systems : monograph / E.A. Novikov, Yu.V. Shornikov. – Novosibirsk : NSTU Publishers, 2012. – 451 p.

ISBN 978-57782-2023- The monograph is concerned with designing unconventional numerical methods to solve the Cauchy problem for stiff systems of ordinary differential equations. Spe cial attention is given to the control of calculation accuracy and numerical scheme stability as well as to the development of variable order and step integration algo rithms. The methodology of hybrid systems is described in detail and their classifica tion is given in the monograph. The capabilities of instrumental environment of hy brid model computer analysis are presented. Peculiarities of the application of the proposed program complex are demonstrated by a number of practical tasks.

The book is intended for a wide range of specialists in applied mathematics and numerical analysis as well as for experts involved in computer simulation of physi cal, chemical, biological and other processes.

UDC 004. ISBN 978-57782-2023-5 © Novikov E.A., Shornikov Yu.V., © SB RAS Computer Simulation Institute, © Novosibirsk State Technical University, ВВЕДЕНИЕ В о многих приложениях возникает проблема численного ре шения задачи Коши для жестких систем обыкновенных диф ференциальных уравнений. Основные тенденции при построении чис ленных методов связаны с расширением их возможностей для решения жестких систем все более высокой размерности [16, 125, 135, 163, 176, 177]. Проблема создания новых эффективных численных методов ре шения жестких задач в связи с этим остается актуальной [16, 17, 19, 40, 54, 87, 163, 169, 176, 207, 227, 228, 282, 283, 286, 287, 303, 305]. При построении эффективных алгоритмов интегрирования требуется раз решить ряд вопросов, которым и посвящена данная монография. Пре жде всего нужно выбрать методы, которые соответствовали бы классу решаемых задач. Здесь в основу алгоритмов интегрирования положены одношаговые безытерационные численные схемы [163, 169, 176, 177, 207]. Такие численные формулы обладают определенными преимуще ствами перед многошаговыми методами, характеризующимися в неко тором смысле осреднением решения – срезанием экстремумов. При моделировании некоторых динамических объектов этот факт делает их неприемлемыми. Кроме того, если правая часть исходной задачи раз рывная, то применение многошаговых методов невозможно. Доста точно полный обзор работ по многошаговым численным методам при веден в [176, 177, 207, 238, 239, 242, 271, 272, 274, 276, 278, 282, 283, 286, 287, 303, 305]. Требование безытерационности позволяет оценить затраты на шаг интегрирования до проведения расчетов и упрощает программную реализацию алгоритмов интегрирования [40, 121, 125, 135, 149, 163, 176, 177].

Важность указанных задач породила за последнее время огромное количество методов интегрирования жестких систем [1–4, 9–12, 16,17, 19, 21–23, 25, 60, 62–64, 71, 77, 81–83, 93, 94, 96, 97, 100–107, 109–130, 132–142, 147, 149, 152, 158–159, 161, 162, 165, 169, 176, 177, 207–209, 239, 252, 258–260, 264, 265, 267–269, 279, 300, 304, 338, 343, 348, 353, 358]. Однако для того, чтобы от идеи метода перейти к его алгоритми ческой реализации, необходимо решить большой круг важных вопро сов, связанных с изменением величины шага интегрирования и оцен кой точности получаемых численных результатов, что и делает метод экономичным и надежным [121, 122, 135, 163, 169, 176, 177, 207]. Со временные способы управления шагом основаны обычно на контроле точности численной схемы [163, 176, 177]. Такой подход хорошо заре комендовал себя во многих случаях и представляется наиболее естест венным, поскольку основным критерием при практических расчетах является точность нахождения решения. Многие алгоритмы интегри рования для выбора величины шага используют оценку локальной ошибки [163]. Это оправдано тем, что если на каждом шаге контроли ровать некоторый минимальный уровень локальной ошибки, то гло бальная ошибка будет ограничена. В настоящее время можно выделить три практических способа оценки данной ошибки.

Классическим способом оценки локальной ошибки одношаговых методов считается способ, основанный на экстраполяционной формуле Ричардсона [163, 176, 177, 330, 331]. Его еще называют правилом Рун ге, и он заключается в следующем. В каждой сеточной точке интервала интегрирования приближенное решение вычисляется с шагом h и h / 2, а искомая оценка определяется через разность приближений к решению. Недостаток такого способа – в необходимости дважды вы числять решение в каждой точке.

Более дешевым выглядит многошаговый способ, предложенный в [242]. Он заключается в том, что одношаговой формуле p-го порядка точности в соответствие ставится многошаговая схема (p + 1)-го по рядка. Затем данная схема преобразуется таким образом, чтобы после подстановки в нее приближений получилась оценка локальной ошибки n, p одношагового метода.

Недостатком данного способа является многошаговость оценки, что приводит ко всем недостаткам многошаговых методов.

В последнее время все чаще оценку локальной ошибки вычисляют с помощью вложенных методов. Приближение к решению в каждой точке вычисляется двумя методами вида p-го и (p + 1)-го порядков точности. Затем локальная ошибка метода p-го порядка оценивается через разность полученных приближений к решению. Обычно такой способ используется тогда, когда для расчетов по методу p-го порядка не требуется дополнительных вычислений правой части и матрицы Якоби. Следует отметить оперативность и относительную дешевизну оценки локальной ошибки с помощью вложенных методов, по затра ВВЕДЕНИЕ там на шаг лежащей между оценкой ошибки с помощью правила Рунге и многошаговой оценкой. В то же время по отношению к многошаго вой оценке в ней при вычислении ошибки используется информация только с данного шага, что повышает ее надежность. Данный способ успешно применялся в работах [43–46, 93, 94, 96, 97, 100–107, 109–130, 132–138, 163, 176, 177, 264 – 266, 276, 293–295, 320–324] и будет ис пользоваться здесь.

Применение оценки локальной ошибки при выборе величины шага интегрирования и для контроля точности вычислений в ряде случаев приводит к успеху. Однако с целью повышения надежности расчетов необходимо найти оценку глобальной ошибки. Наиболее известный способ определения данной ошибки основан на предположении о ли нейном характере накопления глобальной ошибки из локальных по грешностей [271, 272, 274–275, 286, 287].

В последнее время при численном исследовании некоторых жест ких задач все большее внимание привлекают явные методы [31, 81, 121, 131, 135, 160, 163, 176, 177, 207, 266]. Это связано с тем, что при решении ряда задач абсолютно устойчивыми методами возникает про блема с размещением элементов матрицы Якоби в оперативной памяти ЭВМ и, что более существенно, с ее обращением, вернее, декомпози цией. В то же время явные методы не нуждаются в вычислении матри цы Якоби, и если жесткость задачи не слишком велика, то они будут предпочтительнее. Появление многопроцессорных ЭВМ позволяет иначе рассматривать явные методы, которые легко распараллеливают ся [81, 135].

Можно выделить две основные причины, которые приводят к трудностям при применении явных методов для решения жестких за дач. Первая причина связана с противоречием между точностью и ус тойчивостью численной схемы на участке установления. Следствием этого является раскачивание шага интегрирования, что в ряде случаев заканчивается аварийной остановкой вычислений. В лучшем случае раскачивание шага существенно снижает эффективность алгоритма интегрирования. Этого недостатка можно избежать предложенным в [101] способом контроля устойчивости. Алгоритмы интегрирования такого типа построены, например, в [101, 103, 110, 118, 119, 135, 137, 159, 161, 162, 321–324].

Вторая причина ограниченного применения явных методов связана с тем, что области устойчивости известных численных схем слишком малы [135]. В настоящий момент имеется ряд работ, посвященных во просам построения явных методов с расширенными областями устой чивости [83, 96, 134, 136, 138, 215, 218, 220–224, 231, 243–251, 253, 255, 256, 267, 268, 270, 288–292, 296–298, 301, 306, 310, 312–313, 319, 328, 332, 339–341, 343–351, 354–358, 360]. Ясно, что расширение об ласти устойчивости связано с ростом вычислительных затрат на шаг интегрирования. Поэтому, если шаг ограничен по точности, такие схе мы будут малоэффективны. Если же шаг ограничен в силу устойчиво сти, что имеет место на участке установления, то за счет применения численных схем с расширенными областями устойчивости удается значительно повысить эффективность алгоритма интегрирования [71, 81, 97, 102, 111, 112, 135, 158, 325]. В качестве критерия выбора эф фективной численной формулы используется неравенство для контро ля устойчивости [101].

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

Интенсивное исследование неявных методов началось после рабо ты [254], в которой было введено понятие A -устойчивости. Однако требование оказалось слишком обременительным для линейных мно гошаговых методов, и поэтому были введены менее ограничительные определения устойчивости [163, 176, 177, 272]. Среди многошаговых методов наибольшее распространение получили формулы дифферен цирования назад [271, 272, 274, 275, 286, 287], обладающие свойством жесткой устойчивости [163].

Понятие A -устойчивости привело к рассмотрению неявных мето дов типа Рунге–Кутты. Наиболее полное исследование этих методов освещалось в работах [229–230, 232–237, 301], а позднее – в моногра фиях [176, 177, 207, 237]. В [234, 235] доказана теорема о том, что для каждого m существует неявная m -стадийная схема, порядок точности которой равен 2m. Заметим, что аналогичная теорема для явных мето дов типа Рунге–Кутты отсутствует. Согласно [207] функция роста или функция устойчивости схемы максимального порядка является диаго нальной аппроксимацией Падэ для функции exp( x). Если отказаться от максимального порядка, то можно построить схемы с лучшими свойствами устойчивости [176, 177].

Несмотря на хорошие свойства точности и устойчивости, практи ческое использование неявных методов типа Рунге–Кутты ограничено ВВЕДЕНИЕ вследствие больших вычислительных затрат на шаг интегрирования.

При реализации требуется итерационный процесс. Метод простой ите рации не эффективен при решении жестких задач, потому что он при водит фактически к такому же ограничению на размер шага, что и яв ный метод. Поэтому возникает необходимость использования метода Ньютона или каких-либо его модификаций. Это приводит к необходи мости обращения матрицы размерности ( m N m N ), где m есть число стадий, а N – размерность системы. Значительно сократить вы числительные затраты при реализации можно за счет использования одной и той же матрицы на нескольких шагах интегрирования. Замо раживание матрицы становится возможным благодаря тому, что это не влияет на порядок точности численной схемы, а определяет только скорость сходимости итерационного процесса. Поэтому необходи мость в ее пересчете возникает при значительном замедлении сходи мости итерационного процесса.

Трудности с реализацией неявных методов типа Рунге–Кутты, при вели к поискам более простых их модификаций. В работах [78, 163, 176, 177, 237] был рассмотрен класс полуявных формул типа Рунге– Кутты, т. е. таких методов, для которых имеет место ij = 0 при i j.

В этом случае итерационная матрица является блочно-диагональной, причем число блоков совпадает с числом стадий, а размерность каждого блока – с размерностью вектора решения. В результате вместо обра щения матрицы размерности ( m N m N ) надо обратить m матриц размерности N каждая. Исследование полуявных методов содержится в работах [163, 176, 177]. Заметим, что если матрица, составленная из коэффициентов ij, 1 i, j m, имеет одно m-кратное собственное число, то неявный метод можно реализовать с теми же затратами, что и полуявный метод [163]. Вопросам построения диагонально-неявных методов посвящена работа [78]. Дальнейшего сокращения вычисли тельных затрат можно добиться, если положить равными все ii, 1 i m, и аппроксимировать все диагональные матрицы одной. Тогда на шаге требуется обратить только одну матрицу размерности N. Для этого случая в [207] доказана теорема о том, что порядок точности (m + 2) не может быть достигнут ни для какого m-стадийного полуяв ного метода при 11 = … = mm.

В [333] предложены два метода второго и третьего порядков точ ности. Они отличались от явных методов типа Рунге–Кутты регуляри зацией правой части дифференциальной задачи. В дальнейшем такие схемы стали называть методами типа Розенброка.

Свойства методов исследовались в работах [1–4, 9–12, 41, 42, 44, 45, 140, 141, 163, 165, 176, 177, 211, 293, 295, 334, 352]. В [9, 10] приведены явные формулы для определения параметров методов с первого по чет вертый порядок включительно. В [1–4, 165] изучаются методы типа Ро зенброка с комплексными коэффициентами. Наиболее эффективные реализации методов типа Розенброка возникают тогда, когда все коэф фициенты равны между собой, а ij = 0. В этом случае на каждом вре менном шаге требуется вычисление и обращение всего лишь одной мат рицы размерности N. Для улучшения свойств точности методов типа Розенброка в [359] предлагается вычислять стадии по определенным эффективным формулам. Данная модификация получила название ROW-методов. В [211, 293, 295, 333, 334, 352] изучаются частные мето ды такого типа, а также рассматриваются вопросы их реализации. В [30, 127, 129] изучаются возможности реализации таких методов на ЭВМ кластерной архитектуры.

Методы типа Розенброка и ROW-методы относятся к одношаговым безытерационным методам. Они существенно проще с точки зрения реализации, чем алгоритмы на основе численных формул, в которых стадии вычисляются с использованием итераций. Однако в данных ме тодах матрица Якоби влияет на порядок точности численной схемы.

Поэтому недостатки таких методов в основном связаны с проблемой замораживания матрицы, декомпозиция которой при больших значе ниях N фактически определяет общие вычислительные затраты. Если вопрос об использовании одной матрицы на нескольких шагах оста вить нерешенным, то в случае использования безытерационных мето дов заведомо нужно ограничиваться задачами небольшой размерности.

Хотя в [114] на основе схемы типа Розенброка второго порядка точно сти построен алгоритм интегрирования с замораживанием матрицы Якоби, уже при построении такого алгоритма на основе формулы третьего порядка возникают принципиальные трудности. В работах [82, 216] численные схемы, в которых матрица Якоби включена непо ВВЕДЕНИЕ средственно в вычислительные формулы, приведены для решения дифференциально-алгебраических задач.

В [20–23] исследуется класс методов, основанный на дробно рациональном представлении приближенного решения, а в [25, 162] рассмотрены численные формулы, основанные на использовании ап парата цепных и ветвящихся дробей. Еще один подход к построению вычислительных алгоритмов заключается в конструировании чис ленных схем, учитывающих специфику исходной задачи. Здесь мож но выделить схемы экспоненциальной подгонки, а также методы, ба зирующиеся на обращении главной части дифференциального оператора [21].

Новый численный метод интегрирования жестких систем, в основе которого лежит принцип последовательной фильтрации членов, соот ветствующих наибольшим собственным значениям матрицы Якоби системы, предложен в [47]. Согласно [47] этот метод позволяет без по тери устойчивости увеличить шаг интегрирования даже в случае про стейших численных схем. В [93, 104, 108, 128] изучаются адаптивные методы на основе явных и L -устойчивых численных формул. В них эффективная схема выбирается с применением неравенства для кон троля устойчивости. В [100, 132, 142, 147] одношаговые методы рас смотрены применительно к моделированию химических реакций.

В [12] рассмотрены вопросы реализации методов интегрирования на Фортране.

При численном решении задач механики сплошной среды после дискретизации по пространственным переменным методом конечных разностей или конечными элементами возникает проблема решения задачи Коши для системы обыкновенных дифференциальных уравне ний. Для решения таких задач в [107, 123, 124, 133] изучаются адди тивные методы. Там рассмотрены неоднородные численные формулы, часть стадий которых имеет вид явных схем типа Рунге–Кутты, а часть стадий взята из L -устойчивых методов. Кроме задач механики сплош ной среды аддитивные методы можно применять для решения локаль но-неустойчивых задач.

Во второй половине ХХ века в ряде областей техники появились сложные технические системы [18, 74, 75, 180], к которым относятся прежде всего сложные системы управления динамическими объекта ми. Можно выделить их характерные особенности: непрерывную и дискретную составляющие поведения системы;

разнородные физиче ские принципы действия (электрические, механические, гидравличе ские, оптические и др.);

множество информационных и физических связей;

иерархическую многоуровневую структуру и т. д. При анализе таких систем необходимо учитывать наличие быстрых и медленных процессов;

высокую размерность;

существенную нелинейность мате матических моделей, возможно, с разрывными функциями;

жесткость непрерывных моделей;

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

В современной терминологии системы с перечисленными особен ностями принято называть гибридными системами (ГС) [108, 202, 329].

Понятие гибридных систем появилось в 1990-е годы при исследовании различных технических приложений, связанных с многорежимными динамическими процессами в исследовательских лабораториях уни верситета Беркли [65, 66, 75, 307, 308, 318]. В литературе также ис пользуются термины: «непрерывно-дискретные системы», «системы с переменной структурой», «реактивные, событийно-управляемые и со бытийно-непрерывные системы». Специфические особенности ГС ог раничивают использование традиционных аналитических способов анализа таких систем, и в связи с этим методы численного моделиро вания приобретают ведущую роль.

Теоретические вопросы модельного представления гибридных сис тем в классе дифференциально-алгебраических уравнений (ДАУ) об щего вида с индексом 1 и выше и согласованием алгебраических и фа зовых переменных исследуются в [75, 156, 157]. В работах [168, 170, 172] рассматриваются решения, которые проходят по поверхности раз рыва в системах автоматического управления, и исследованы свойства решений таких систем. Теоретические основы дифференциальных уравнений с разрывной по времени правой частью можно найти, на пример, в [153, 171]. Современному состоянию теории уравнений с разрывной правой частью полностью посвящена книга [171]. Основной математической моделью как для гибридных систем, так и для систем релейного и импульсного управления служат дифференциальные уравнения с разрывами первого рода в первой производной фазовых переменных. В таких системах, если решение принадлежит поверхно сти разрыва (случай неодносторонних событий в ГС), то после попада ВВЕДЕНИЕ ния на нее изображающей точки дальнейшее движение может оказать ся невозможным, так как правая часть не определена на поверхности разрыва. В этом случае правая часть требует доопределения. В работе А.Ф. Филиппова было предложено формально доопределять правую часть способом, который сейчас носит название доопределения Фи липпова [171]. В этой работе дается определение, что считать решени ем системы уравнений с разрывной правой частью, рассчитываются условия существования решения и выясняется, какие свойства класси ческих динамических систем сохраняются при появлении разрывов в правой части [156, 157].

В.И. Уткин [168] анализирует различные способы доопределения, вводит понятие эквивалентного управления и показывает, в каких слу чаях оно совпадает с доопределением Филиппова. Уткин применяет его для анализа скользящих режимов при скалярном и векторном управлении. Предложенный В.И. Уткиным и А.Ф. Филипповым способ относится к аналитическим методам исследования гибридных систем, который практически не применим для численного исследования ГС, потому что движение по поверхности скольжения сопровождается эф фектом Зенона [156].

Введение нового состояния с уравнениями движения, аппроксима ция разрывных функций в правой части, обусловленная реальным по ведением любого инерционного объекта, часто устраняют причину возникновения эффекта Зенона в окрестности линии скольжения при численном интегрировании. Однако выбор аппроксимирующей функ ции не всегда очевиден и во многих случаях является нетривиальной задачей. В классическом численном методе без алгоритма корректного обнаружения событий в ГС моменты переключения находятся неточ но, что позволяет изображающей точке переходить из одной области допустимых значений в другую. Такие переходы недопустимы для од носторонних гибридных систем. Для устранения этого в ГС решается задача корректного обнаружения событий, связанных с точкой пере ключения [74, 75, 156, 157, 262, 307, 327].

Изученные в работах отечественных ученых явления скольжения в динамических системах с разрывными правыми частями – аналог эф фекта Зенона при численном анализе гибридных систем. Теоретиче ские результаты, полученные для непрерывных систем с разрывами, обосновывают существование и единственность решения и коррект ность вычислительных экспериментов при численном анализе гло бального поведения гибридной системы.

Отметим, что моделями, которые идентифицированы разными дифференциальными уравнениями в различных областях фазового пространства, занимаются достаточно давно в теории управления [27, 68, 168, 170–172]. В указанных работах исследуется поведение слож ных нелинейных задач, которые подпадают под современное опреде ление гибридных систем. Примеры таких систем [18, 27, 68, 74, 75, 155–157, 214, 299, 307, 309] подробно проанализированы в моногра фии [157] и будут рассмотрены ниже. Несомненно, ценные теоретиче ские результаты связаны с определенными ограничениями на порядок системы уравнений моделей и вид правой части. Однако анализ моде лей в перечисленных работах не касается современных методов чис ленного моделирования и использования инструментальных средств машинного анализа. В [261–263] показано, что в общем случае не су ществует аналитического общего решения таких систем. Этот тезис поддерживается и другими специалистами по анализу гибридных сис тем [75], которые отмечают, что попытки применить классические подходы к анализу гибридных систем [219] пока дают весьма ограни ченные результаты. В связи с этим анализ таких систем необходимо проводить с использованием аппарата численных или численно аналитических методов [29] в окружении инструментальных средств [85, 213].

Несмотря на то что в любой физической системе время t R ап риори является величиной непрерывной, можно говорить о дискрет ных моментах времени t * как о подмножестве значений непрерывного времени. Такие моменты принято называть временной щелью (time gap) [318]. Задача определения момента времени t = t *, когда событий ная функция g ( x(t ), t ) = 0, является достаточно сложной задачей кор ректного обнаружения событий [74, 75, 157, 262, 307].

Существует набор контекстов, в которых область определения пра вой части системы дифференциальных уравнений ограниченна. В од них задачах векторное поле определено не везде, поскольку модель спроектирована для применения только в некоторых областях фазово го пространства. Например, нельзя применять сверхзвуковую аэроди намическую модель летательного аппарата на дозвуковых скоростях [75]. В других ситуациях область неопределенности является следст ВВЕДЕНИЕ вием физической природы проблемы. Так, в задачах химической кине тики фазовые переменные соответствуют концентрациям различных компонентов. Их отрицательные значения не имеют физического смысла, но могут быть получены вследствие численных ошибок моде лирования или неадекватности модели. Подобные ситуации возникают в проблемах с фазовыми переходами. Такие модели характеризуются тем, что имеют области неопределенности в фазовом пространстве. В этом случае управляющий компьютер должен переключиться на дру гой закон управления в этой области. Пример – контроллер, который вырабатывает управляющее решение в зависимости от результатов расчета инверсной кинематики для руки робота [261]. Если требуемая конфигурация управляющего органа недопустима, то решение уравне ний инверсной кинематики не существует.

Очевидно, что, если в процессе моделирования гибридных систем неправильно сработает алгоритм обнаружения события, это может па губно отразиться на глобальном решении. Впервые требования о необ ходимости корректного обнаружения дискретных событий для гиб ридных моделей были сформулированы в работах [318, 326].

Большинство гибридных симуляторов [157, 227, 311] разделяют за дачу обнаружения на два этапа – фазу обнаружения, или детекции, за которой следует фаза локализации. Фаза локализации обычно основана на методе деления отрезка пополам, методе дихотомии или использо вании метода Ньютона–Рафсона для поиска корней событийной функ ции на границе режима. В отдельных случаях прибегают к методу ус тановления для поиска корней. Но эти методы не всегда эффективны, если их применять к так называемым односторонним событиям гиб ридной системы. Эта проблема рассматривается в работах [65, 74, 75, 157, 261–263, 281, 307, 308]. Однако не все вопросы решаются в этих работах. В частности, проблема становится еще более актуальной для режимов ГС с повышенной жесткостью [191–193, 197].

Наиболее опасна для моделирования ситуация, когда переходный участок решения лежит вблизи границы области неопределенности и якобиан событийной функции резко возрастает вблизи границы. Это может привести к «проскакиванию» точки переключения с большей вероятностью, чем в гладких режимах. И в этом случае ситуация наи более опасна для ГС с односторонними событиями. Как только собы тие локализовано, интегратор останавливается и происходит переход в другое состояние. Хотя этот базовый метод, впервые описанный в [241], оказывается пригодным для многих задач, возможны ситуации, когда он склонен к сбоям. Эти ситуации связаны с многократными пе ресечениями траектории на границе событийной функции либо с «ост рой» геометрией границ.

В качестве второго класса задач, для которых стандартные методы работают неудовлетворительно, можно рассмотреть случай, когда пра вая часть дифференциального уравнения не определена для некоторых значений x, при которых g ( x) 0. Следует также отметить, что неко торые проблемы могут возникать при поиске точек переключения в так называемых «мультиагентных» моделях гибридных систем [26, 75, 361], в которых независимо и параллельно функционирует большое число взаимодействующих объектов.

Ф. Келлер (F. Cellier) был первым [241], кто отметил, что события состояний требуют особого рассмотрения, и предложил остановку при приближении к точке разрыва. К. Гир (C.W. Gear) [273] продемонстри ровал неэффективность, которая возникает, если не использовать спе циальные методы. М. Карвер (M.B. Carver) [240] первым отметил, что интенсивность изменений событийной функции вдоль поля движения есть важнейшая характеристика для обнаружения события.

Идея диф ференцирования событийной функции и введения ее как дополнитель ной переменной состояния была предложена им же. В каждом из этих случаев события обнаруживались простым наблюдением за сменой знака g (t, x ) на каждом шаге интегрирования. При таком подходе со бытия не обнаруживаются в случае многократной смены знака за один шаг. Работая над данной проблемой, Л. Шампайн (L.F. Shampine) и его коллеги [335] использовали тот факт, что интерполяционные полино мы могут быть сгенерированы для событийной функции и применены для правильного определения расположения событий с помощью по следовательностей Струма (Strum sequences). Однако они не использо вали эту информацию для выбора величины шага интегрирования.

Несколько простых алгоритмов для обнаружения событий в ДАУ описаны в [157]. Эти методы могут обнаруживать многократные пере ходы, однако они имеют склонность к большим вычислительным за тратам. Т. Парк (T. Park) и П. Бартон (P.I. Barton) [326] объединили некоторые из этих идей и использовали интервальную арифметику для создания эффективных тестов для определения интервала, где, воз можно, произошло событие. Этот метод обнаружения событий выгля дит наиболее надежным, он рациональный и хорошо приспособлен к ВВЕДЕНИЕ противостоянию проблеме. Несмотря на то что все эти методы исполь зуют остановку при приближении к точке разрыва, ни один из них не обеспечивает алгоритма выбора величины шага таким образом, чтобы состояние никогда не пересекало границу события. Поэтому они не смогут локализовать событие, произошедшее в окрестности особой области модели, которую принято [101] также называть погранслоем.

С учетом обозначенных проблем Дж. Эспозито (J. Esposito) [261] предложил применить аналогию с управлением и рассматривать моде лируемую систему как систему управления: величина шага интегриро вания есть вход, а событийная функция – выход. Проблема заключает ся в выборе такого закона обратной связи, чтобы система приближалась к границе события g ( x) = 0 асимптотически, никогда не пересекая ее. Так как решение приближается к погранслою асимптоти чески, то увеличиваются шансы обнаружить событие без пересечения границы. В этом случае нет риска попадания в особые точки, где гиб ридная модель может быть не определена.

Однако предложенный алгоритм по выбору шага для асимптотиче ского приближения решения к погранслою не учитывает критерия ус тойчивости методов численного интегрирования, что особенно важно при анализе гибридных систем с жесткими режимами. Задача состоит в корректном определении точки переключения с учетом ограничений на событийную функцию и жесткости режима. В настоящее время су ществуют два конструктивных подхода для построения алгоритма, ко торый гарантирует правильное решение задачи Коши с ограничения ми. Первый алгоритм основан на идее линеаризации [261], по аналогии с методом линеаризации в САУ. Второй алгоритм основан на методе установления [17] и связан с поиском корней событийной функции g ( x) = 0 на решении классической задачи Коши без ограничений. Од нако решение, полученное методом линеаризации, в общем случае не совпадает с решением, полученным методом установления. По мне нию автора этой идеи [157], оба решения верные, и в этом легко убе диться для линейной правой части в исходной задаче Коши. Совпаде ние решений, полученных разными методами, может происходить только при определенном значении соответствующих параметров, что противоречит общему случаю из области определения этих пара метров.

Для анализа жестких режимов ГС разработан адаптивный алгоритм DISPF1_RADAU5, в основе которого лежит оригинальный явный метод переменного порядка и шага DISPF и известный неявный метод RADAU5 [177]. Эффективность этого алгоритма повышена путем вве дения контроля жесткости и выбора эффективной численной схемы в зависимости от текущего решения. Ограничение времени работы неяв ного метода позволяет сократить вычислительные затраты.

Для анализа жестких режимов ГС практически во всех инструмен тальных средствах используются неявные методы. Вопросы эффектив ного использования явных методов переменного порядка и шага для исследования жестких задач обсуждаются в [99, 131, 179, 200, 337].

Следует также отметить, что во многих практических гибридных сис темах при анализе жестких режимов эффективными оказываются ком бинированные явно-неявные методы или, как их еще называют, неод нородные схемы [93, 320]. Основой алгоритма является осуществление контроля устойчивости [184], а приближение к границе устойчивости является критерием переключения схем. Отметим, что комбинирован ные схемы могут быть наиболее эффективно применены для решения ГС с разнородными жесткими и нежесткими режимами [164].

Новые алгоритмы анализа сложных динамических систем оказы ваются более эффективными для прикладного специалиста, если они окружены специализированными инструментальными средствами [85, 86, 91, 95, 192, 204, 213]. Индустрия создания таких средств моделиро вания становится сама по себе фундаментальной задачей исследова ния. Сейчас масштаб и объем трудностей при создании инструмен тальных средств настолько вырос, что можно говорить, что задача их преодоления сама стала задачей науки и представляет собой проблему фундаментального значения [210, 277]. Универсальные передовые оте чественные MVS [74, 75, 155–157, 219], AnyLogic [65, 226] и зарубеж ные DYMOLA [http://www.dynasim.se/Dymola, 164, 227, 257], Ptolemy II [308] и HyVisual [http://ptolemy.eecs.berkeley.edu/hyvisual], HyTech [284], Charon [261–263, 307], Simulink/Stateflow [56, 342] программные комплексы моделирования гибридных систем широко используются для анализа сложных динамических процессов. Однако с помощью этих программных комплексов в отдельных случаях не удается полу чать качественные результаты при решении важных практических за дач [202]. Рассмотренные здесь расчеты, редактирование, структурные преобразования и анализ результатов проведены с применением среды ИСМА (инструментальные средства машинного анализа) [61]. Тексто вые модели представлены на оригинальном языке LISMA [203].

ВВЕДЕНИЕ Важным этапом компьютерного анализа ГС является композиция компьютерных моделей. В современных инструментальных программ ных комплексах это связано с двумя аспектами – визуальное и сим вольное представления математической модели, которые непосредст венно связаны с вопросами формализма и спецификации ГС.

Несомненными передовыми технологиями стали графические языки спецификаций предметных категорий. Для программных моделей ГС это карты состояний Харела (statechart) [280, 281], канонизированные в проекте UML [24, 28, 225] и успешно развитые в системах HyVisual, MVS и др. Они показали себя удобным и наглядным изобразительным средством представления визуальной модели ГС в части дискретного поведения. Узлами диаграмм Харела являются локальные состояния ГС. Направленные дуги с предикатами событийных функций показы вают переходы из локальных состояний. В интерфейсах карт поведе ния ГС программные модели содержат общепринятые декларации всех фазовых, алгебраических и булевых переменных, что не относится по существу к компьютерной модели, а является необходимым атрибутом программирования. Для систем высокой размерности [164, 198, 317] сектор описания типов переменных может быть соизмерим с матема тическим описанием. Поэтому бездекларативный язык [202] в этом смысле более лаконичен и доступен предметному пользователю. Сле дует отметить, что многие современные графические оболочки ис пользуют вместе с тем и другие формализмы с соответствующим гра фическим языком. Например, сети Петри [88, 145, 146, 309] в системе DYMOLA, структурные схемы в системах HyVisual и Simulink обла дают своими функциональными преимуществами с точки зрения предметной ориентации пользователя. Структурные схемы отражают операторные преобразования Лапласа [151] и являются традиционным графическим языком представления моделей динамических систем с огромным накопленным арсеналом аналитических методов исследова ния локальных непрерывных поведений ГС [7, 8, 29, 89, 168, 170–172], унаследованных в основном из теории автоматического управления [33, 49, 57, 68, 73]. Представление непрерывного поведения структур ными моделями позволяет более детально представить непрерывную часть модели и эффективно организовать активный вычислительный эксперимент [90], что весьма важно при отладке моделей и параметри ческой верификации.

Символьный язык [69, 84] – неотъемлемый атрибут спецификации, он сопровождает графические конструкции, например, язык MODELICA [315, 316] описывает гибридную модель в целом. Выбор соответствующего символьного языка и средств его эффективной реа лизации также является важной проблемой разработки программных систем [59].

Способы визуальной интерпретации результатов вычислительного эксперимента в современных зарубежных и особенно в передовых оте чественных системах моделирования ограничены в части манипуляции графическими и числовыми данными, полученными в результате ре шения. В частности, ограничен режим катенации окон с графическими данными, импорт данных из внешних приложений, трассировка точеч ных решений, интерполяция графических данных, например, с помо щью вейвлет-преобразований [32, 50, 55, 181, 186, 217]. В то же время все перечисленные вопросы широко востребованы в практике анализа результатов вычислительного эксперимента и недостаточно решены в современном инструментарии.

Еще одним важным аспектом при разработке программного обес печения (ПО) является его унификация. Унификация программного обеспечения предполагает использование разработанного языкового процессора к предметным задачам из разных областей с настройкой предметного языка нового приложения как подмножества базового языка.

Предложенная вниманию читателей монография состоит из 14 глав.

Глава 1 посвящена определению основных понятий, которые рас сматриваются в дальнейшем. Без доказательства приведены теоремы существования и единственности решения, зависимости решения от начальных данных, устойчивости и асимптотической устойчивости.

Основные определения и теоремы взяты из книг [40, 144, 152, 163, 169, 176, 177, 182, 183, 207].

В главе 2 рассматриваются вопросы нахождения приближенного решения задачи Коши для систем обыкновенных дифференциальных уравнений одношаговыми численными методами. Предложены два способа оценки аналога глобальной ошибки, которые не приводят к значительному увеличению вычислительных затрат: с вычислением матрицы Якоби и без вычисления. Предложен способ получения нера ВВЕДЕНИЕ венства для контроля устойчивости методов с ограниченной областью устойчивости, что существенно повышает надежность расчетов [101– 103, 128]. Данное неравенство основано на оценке максимального соб ственного числа матрицы Якоби, которое определяется с использова нием ранее вычисленных стадий.

В главе 3 построены алгоритмы переменного шага на основе явных методов второго, третьего и четвертого порядков точности. Для кон троля точности вычислений и при выборе величины шага интегриро вания используется оценка аналога глобальной ошибки.

В главе 4 построены неравенства для контроля устойчивости явных методов типа Рунге–Кутты со второго по пятый порядок точности.

Глава 5 посвящена построению алгоритмов интегрирования пере менного порядка и шага. На основе явных схем типа Рунге–Кутты с расширенными областями устойчивости разработаны алгоритмы пере менного порядка и шага для решения умеренно жестких задач.

В главе 6 изучаются некоторые проблемы, которые возникают при замораживании матрицы Якоби в методах типа Розенброка. Данные методы получили широкое распространение благодаря достаточно хо рошим свойствам точности и устойчивости и простоте реализации на ЭВМ. Доказана теорема о максимальном порядке точности методов типа Розенброка с замораживанием матрицы Якоби и построен метод максимального (второго) порядка точности.

В главе 7 описан класс ( m, k ) -методов решения жестких задач [109, 116, 120, 125]. Основное отличие (m, k ) -методов от существующих численных формул заключается в том, что в данных численных схемах стадия метода не связывается с обязательным вычислением правой части исходной системы дифференциальных уравнений.

В главе 8 исследованы (m, 3) -методы решения жестких задач. До казана теорема о максимальном порядке точности. Получены коэффи циенты A - и L -устойчивых численных схем пятого порядка.

В главе 9 введены основные определения гибридных систем [202], режима ГС, событийной функции и т. д., проводится классификация событий. Здесь же определяется локальное и глобальное поведение гибридных систем в виде совокупности дифференциально-алге браических уравнений, анализируются системы с разрывами первого рода, приводятся необходимые и достаточные условия существования и единственности решения в таких системах.

В главе 10 рассматриваются вопросы корректного обнаружения со бытий в гибридных системах с учетом жестких режимов и в условиях односторонних событий. Доказывается теорема о выборе шага с уче том жесткости и гибридности модели, определенной задачей Коши с ограничениями на событийную функцию.

В главе 11 для задач повышенной жесткости разработан адаптив ный алгоритм DISPF1_RADAU [184]. В адаптивном алгоритме DISPF1_RADAU использован оригинальный явный метод DISPF, рас смотренный ранее и известный L -устойчивый метод RADAU5 [176].

Из результатов вычислительных экспериментов следует, что эффек тивность разработанного адаптивного алгоритма в 5–10 раз выше, чем в лучших мировых инструментах анализа жестких режимов ГС.

Глава 12 посвящена программному обеспечению исследования ди намических и гибридных систем. Приводится структурная специфика ция [194] для описания непрерывной части моделей гибридной систе мы. Введены макроструктуры [143, 201] графического или визуального языка непрерывных моделей ГС. Показана организация нелинейной функции с режимом импорта данных внешних приложений [5, 6, 190, 217], которая выгодно отличает возможности пограммного приложе ния от известных аналогов [61]. При анализе структурных моделей ус танавливается корректность связей, контролируются алгебраические петли [70, 75, 157] рекурсивным алгоритмом в глубину [154], форми руется орграф программной модели непрерывной части [36]. Диаграм мы Харела рассматриваются с позиции теории графов [51, 142]. Дис кретная составляющая гибридной системы специфицирована также символьным описанием. Символьные компьютерные модели могут специфицировать и непрерывную часть ГС [314], если непрерывное поведение задано дифференциальными или дифференциально-алгеб раическими уравнениями. Языковые конструкции разрабатываются по математическим моделям ГС с требованием максимального сохране ния естественной формы и доступности для непрофессиональных пользователей в программировании. Язык LISMA [188, 203, 205] раз рабатывается на бездекларативной основе. По языку строится контек стно-свободная порождающая грамматика [34, 52]. Кроме порождаю щих грамматик, механизмом порождения лексемных конструкций языка являются регулярные выражения [14, 206]. Грамматика арифме тических и логических выражений преобразуется путем устранения ВВЕДЕНИЕ левой рекурсии [15, 66, 150] к эквивалентной однозначной грамматике LL(1) [13, 14, 37], которая эффективно анализируется методом рекур сивного спуска [39, 206]. Наконец, рассматривается наиболее общая форма структурно-символьной спецификации [195] или визуально лингвистическое описание ГС.

Графическая интерпретация вычислительного эксперимента вы полнена во временной и в фазовой плоскостях [199] встроенным гра фическим интерпретатором решений, который, в отличие от мировых аналогов, позволяет проводить вейвлет-преобразования [32, 50, 53, 55, 181, 185, 209], катенацию окон и другие манипуляции с числовыми и графическими данными результатов моделирования.

В главе 13 рассмотрены вопросы унификации программного обес печения. Предметом унификации является возможность исследования новых приложений, модели которых идентифицированы в обозначен ном классе систем дифференциально-алгебраических уравнений.

В качестве примера такого приложения выбрана химическая кинетика, в которой математической моделью для анализа концентрации реаген тов химических реакций являются системы ОДУ. Для спецификации задач химической кинетики [166] разработан язык LISMA+, являю щийся расширением базового языка спецификации сложных динами ческих и гибридных систем LISMA. Для грамматики справедливо включение GLISMA + GLISMA, что обеспечивает преемственность языков спецификации новых приложений с наследуемыми методами обработки входного текста, при этом не требуется изменение интер претатора ПО. При разработке языка химических реакций LISMA+ рассмотрены практические примеры химических реакций пиролиза этана, дифференциации растительной ткани и др. [48, 58, 72, 176, 302].

Здесь же учтены недостатки используемых методик расчета концен траций в существующем программном обеспечении [13, 14, 68, 178].

Глава 14 посвящена инструментально-ориентированному анализу практических моделей гибридных систем разной природы. Модель импульсной системы автосопровождения баллистических и космиче ских ракет [7, 8, 38] наглядно и просто идентифицирована символьной программной моделью в идеологии ГС и успешно проанализирована разработанными инструментальными средствами. При анализе собы тийно-управляемой системы кольцевого модулятора [285] показана эффективность алгоритма MK22 по сравнению с известным методом Гира в пакете MAPLE [35, 273], причем особенности представленной гибридной модели резко ограничивают использование лучших миро вых инструментов. При исследовании живой билиарной системы ис пользовано описание потоков, взятое из работ [18, 92, 174, 175, 187, 189, 196, 212, 336] по аналогии с гемодинамикой [76]. Сконструиро ванная гибридная модель эффективно исследована во временной и фа зовой плоскости. Конструктивно доказана эффективность явных мето дов при исследовании гибридной системы реакции-диффузии [98] высокой размерности. Разработанная компактная программная модель из 400 уравнений на языке LISMA проанализирована и неявной схемой [198], но при этом получено неверное глобальное решение, что прак тически доказывает преимущество явных методов при исследовании жестких гибридных систем высокой размерности. Успешно проанали зированы в инструментально-ориентированной идеологии гибридных систем и другие системы из области химической кинетики [166].

Основные результаты монографии получены при поддержке Российского фонда фундаментальных исследований (проект № 11-01-00106).

Глава НЕКОТОРЫЕ СВЕДЕНИЯ ИЗ ТЕОРИИ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ Д анная глава посвящена определению тех понятий, которые будут рассматриваться в дальнейшем. Без доказательства приведены теоремы существования и единственности решения, зави симости от начальных данных, устойчивости и асимптотической ус тойчивости.

1.1. ПОСТАНОВКА ЗАДАЧИ Рассмотрим систему обыкновенных дифференциальных уравне ний первого порядка вида y = f (t, y ), (1.1) где f : R R N R N, y : R R N – вектор-функции, t – независимая переменная.

Определение 1.1. Решением системы (1.1) на интервале (a, b) называется непрерывно дифференцируемая на (a, b) функция y : (a, b) R N, обращающая (1.1) в тождество.

Важное место в теории дифференциальных уравнений занимает задача Коши, или начальная задача. Для системы (1.1) задача Коши ставится следующим образом. Среди всех решений системы (1.1) тре буется найти такое решение y (t ), что имеет место y (t0 ) = y0, (1.2) где y0 = ( y1, y2,…, y N )T – заданный вектор, верхний индекс T означа 00 ет транспонирование. Вектор y0 называется начальным вектором, или начальными данными задачи Коши, а условие (1.2) – начальным усло вием.


Г л а в а 1. НЕКОТОРЫЕ СВЕДЕНИЯ ИЗ ТЕОРИИ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ 1.2. СУЩЕСТВОВАНИЕ И ЕДИНСТВЕННОСТЬ РЕШЕНИЯ Существование решения задачи (1.1), (1.2) регулируется теоре мой.

Теорема Пеано. Если правые части системы (1.1) определены и непрерывны в некоторой замкнутой области G, то через каждую внутреннюю точку (t, y1, y2,…, y N )T области G проходит хотя бы одна интегральная кривая.

Теорема Пеано гарантирует существование решения задачи Коши, однако решение может быть не единственным. Например, для задачи Коши y = 3 y 2, y (0) = 0, можно выписать не менее двух решений, т. е. y1 (t ) = 0 и y2 (t ) = t 3 / 27. Единственность решения задачи Коши можно гарантировать при некоторых дополнительных предпо ложениях.

Локальная теорема Коши–Пикара. Пусть для системы уравне ний (1.1) поставлена задача Коши с начальным условием (1.2) и пусть правая часть системы (1.1) определена в замкнутой облас ти E : | t t0 | a, || y y0 || b, причем точка (t0, y0 ) лежит внутри области E. Пусть выполнены условия:

1) вектор-функция f (t, y ) непрерывна в области E, т. е.

существует число M 0 такое, что || f (t, y ) || M в любой точке об ласти E ;

2) функция f (t, y ) удовлетворяет условию Липшица, т. е.

|| f (t, y1 ) f (t, y2 ) || L || y1 y2 ||, где L – некоторое положительное число, называемое постоянной Липшица, а y1 и y2 – любые точки области E.

Тогда задача (1.1), (1.2) имеет единственное решение y = y (t ).

Это решение определено и непрерывно дифференцируемо на интерва ле | t t0 | h, h = min(a, b / M ) и не выходит из области E, т. е.

|| y (t ) y0 || b при | t t0 | h.

1.3. Зависимость решения от начальных данных Замечание. Условие Липшица выполнено, если все компоненты вектор-функции f (t, y ) имеют в области E ограниченные частные производные по всем компонентам вектора y, т. е.

| fi (t, y1, y2, …, y N ) / y j | C, 1 i, j N, где C – некоторое положи тельное число.

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

Во многих практических задачах начальные данные получены экспериментальным путем. Погрешность измерения начальных дан ных не должна существенно влиять на решение задачи (1.1), (1.2).

В связи с этим возникает естественный вопрос: будет ли малому изме нению начальных данных соответствовать малое изменение решения?

Другими словами, насколько устойчиво решение к малым возмущениям начальных данных?

1.3. ЗАВИСИМОСТЬ РЕШЕНИЯ ОТ НАЧАЛЬНЫХ ДАННЫХ Рассмотрим задачу Коши для системы дифференциальных урав нений, в которой правая часть и начальные данные зависят от парамет ров z, z = ( z1, z2,…, zm )T. Такая задача имеет вид y = f (t, y, z ), y (t0, z ) = y0. (1.3) Теорема о непрерывной зависимости и дифференцируемости решения по параметру. Пусть в некоторой области G пространст ва (t, y, z ) вектор-функция f имеет непрерывные производные по t и y порядка p, p 1. Тогда решение y (t, z ) задачи Коши (1.3) непре рывно и имеет p непрерывных производных по t и z.

Теорема о непрерывной зависимости и дифференцируемости решения по начальным данным. Пусть в задаче Коши y = f (t, y ), y (t0 ) = y0, (1.4) Г л а в а 1. НЕКОТОРЫЕ СВЕДЕНИЯ ИЗ ТЕОРИИ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ вектор-функция f имеет непрерывные производные порядка p, p по t и y. Тогда решение y (t, t0, y0 ) непрерывно и имеет p непрерыв ных производных по t, t0 и y0.

Начальные данные t0 и y0 можно рассматривать как параметры.

С помощью замены переменных переходим к задаче Коши, правая часть которой зависит от начальных данных, как от параметров. В ре зультате выполнены условия предыдущей теоремы.

1.4. УСТОЙЧИВОСТЬ ПО ЛЯПУНОВУ Непрерывная зависимость от начальных данных гарантирует, что полученные при достаточно близких начальных данных решения бу дут близки, по крайней мере, в начальный момент. С течением време ни, т. е. при возрастании независимой переменной, возможны две си туации: решения при различных начальных данных сближаются или решения все больше расходятся. Теория устойчивости Ляпунова по зволяет выяснить характер поведения решений, соответствующих ма лым возмущениям начальным данных.

Рассмотрим автономную систему вида y = f ( y ), (1.5) где y (t ) и f ( y ) – N -мерные вектор-функции. Обозначим через y (t, y0 ) решение этой системы, соответствующее начальным условиям y (t0 ) = y0. (1.6) Нетрудно убедиться, что при условии f ( y0 ) = 0 система (1.5) с на чальными данными (1.6) имеет стационарное решение y (t, y0 ) y0.

С точки зрения механики можно сказать, что система находится в со стоянии равновесия.

Определение 1.2. Положением равновесия системы (1.5) называ ется такая точка a R N, что имеет место f (a ) = 0.

Если при небольшом отклонении от положения равновесия полу ченное возмущенное решение остается вблизи точки равновесия, будем говорить, что положение равновесия устойчиво.

1.5. Функция Ляпунова Определение 1.3. Положение равновесия a называется устойчи вым по Ляпунову, если выполнены два условия:

1) существует 0 0 такое, что если || y0 a || 0, то решение y (t, y0 ) существует при всех t, 0 t ;

2) для любого 0 существует = () 0 такое, что если || y0 a ||, то || y (t, y0 ) a || при всех t, 0 t.

Определение 1.4. Положение равновесия a называется асим птотически устойчивым по Ляпунову, если оно устойчиво по Ляпунову и имеет место lim y (t, y0 ) = a при достаточно малых || y0 a ||.

t + Асимптотическая устойчивость означает, что решение, получен ное при небольших возмущениях начальных данных, с течением вре мени стремится к стационарному решению.

1.5. ФУНКЦИЯ ЛЯПУНОВА Устойчивость системы (1.5) исследуется с помощью некоторой функции, обладающей наперед заданными свойствами.

Определение 1.5. Непрерывно дифференцируемая в окрестности U точки y = a функция V ( y ) : R N R называется положительно оп ределенной в области U, если V ( y ) 0, y a;

V (a ) = 0, и отрицатель но определенной, если V ( y ) 0, y a;

V (a ) = 0.

Определение 1.6. Производной в силу системы (1.5) функции V ( y ) : R N R, непрерывно дифференцируемой в окрестности U N V ( y ) точки y = a, называется величина V ( y ) = fi ( y ).

i =1 yi Определение 1.7. Положительно определенная в окрестности U точки a функция V ( y ) называется функцией Ляпунова системы (1.5), если V ( y ) 0, y U, где V ( x ) – производная в силу системы (1.5).

Г л а в а 1. НЕКОТОРЫЕ СВЕДЕНИЯ ИЗ ТЕОРИИ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ Теорема Ляпунова об устойчивости. Если в некоторой окрест ности U положения равновесия a существует функция Ляпунова V ( y ), то это положение равновесия устойчиво по Ляпунову.

Теорема Ляпунова об асимптотической устойчивости. Пусть в некоторой окрестности U положения равновесия a существует функция Ляпунова V ( y ) такая, что функция V ( y ) отрицательно оп ределена в U. Тогда положение равновесия a асимптотически ус тойчиво.

1.6. УСТОЙЧИВОСТЬ ПОЛОЖЕНИЙ РАВНОВЕСИЯ ЛИНЕЙНОЙ СИСТЕМЫ Рассмотрим линейную однородную систему с постоянными ком плексными коэффициентами y = Ay, (1.7) где A – матрица размерности ( N N ) с постоянными коэффициентами.

Критерий устойчивости линейной задачи. Положение равнове сия y = 0 системы (1.7) асимптотически устойчиво тогда и только тогда, когда вещественные части всех собственных значений матри цы A отрицательны. При этом имеет место соотноше ние || y (t ) || C || y0 || et, где C – некоторая вещественная постоян ная, y (t ) – решение системы (1.7) с начальными данными y (0) = y0, = max Re(i ) +, i, i = 1, 2,…, N, – собственные числа матрицы 1i N A, а число 0 можно выбрать сколь угодно малым.

В случае, когда элементы матрицы A вещественны, все коэффи циенты характеристического многочлена a0 N + a1 N 1 + …+ a N = 0 (1.8) матрицы A также вещественны. Тогда можно исследовать систему (1.7) на асимптотическую устойчивость с помощью критерия Рауса– Гурвица, не находя корни уравнения (1.8).

1.7. Устойчивость положений равновесия нелинейных систем Критерий Рауса–Гурвица. Пусть все коэффициенты уравнения (1.8) вещественны и a0 0. Для того чтобы все корни уравнения (1.8) имели отрицательные вещественные части, необходимо и достаточ но, чтобы все главные миноры определителя a1 a0 0 0 a3 a2 a1 0 0 a5 a4 a3 a2 0 a aN 2 N 1 a2 N 2 a2 N 3 a2 N 4 были положительны.

1.7. УСТОЙЧИВОСТЬ ПОЛОЖЕНИЙ РАВНОВЕСИЯ НЕЛИНЕЙНЫХ СИСТЕМ Исследование устойчивости нелинейной системы будем сводить к исследованию устойчивости линейной системы. Для этого применим метод линеаризации. Пусть a – положение равновесия автономной системы y = f ( y ), (1.9) где y (t ) и f ( y ) – N -мерные вектор-функции, причем f ( y ) дважды непрерывно дифференцируема в некоторой окрестности U точки a.

Разложим вектор-функцию f ( y ) по формуле Тейлора в окрестности точки a с учетом соотношения f (a) = 0, т. е.

f ( x) = f ( a)( x a ) + g ( x), (1.10) где f (a) = {fi (a) / y j }, 1 i, j N, есть матрица Якоби системы (1.9), а для функции g ( y ) выполняется неравенство || g ( y ) || T || y a ||2, y U. Если точки y и a достаточно близки, то имеет место соотно Г л а в а 1. НЕКОТОРЫЕ СВЕДЕНИЯ ИЗ ТЕОРИИ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ шение f ( y ) f (a )( y a ). Оставляя в разложении (1.10) только ли нейную часть, получим линейную однородную систему с постоянными коэффициентами z = Az, (1.11) где z = y a, A = {fi (a ) / y j } – матрица Якоби системы (1.9).

Система (1.11) называется линеаризованной для системы (1.9) в окрестности положения равновесия a, а переход от нелинейной сис темы (1.9) к линейной системе (1.11) называется линеаризацией. Связь между устойчивостью положения равновесия линеаризованной и ис ходной системы устанавливает следующая теорема.


Терема Ляпунова об устойчивости по линейному приближе нию. Пусть a положение равновесия автономной системы y = f ( y ) и вектор-функция f ( y ) дважды непрерывно дифференци руема в некоторой окрестности точки a. Если вещественные части всех собственных значений матрицы Якоби f (a ) = {fi (a ) / y j } отрицательны, то положение равновесия a асимптотически устой чиво.

Другая формулировка теоремы Ляпунова об устойчивости по линейному приближению. Если положение равновесия линеаризован ной системы асимптотически устойчиво, то асимптотически устой чиво положение равновесия нелинейной системы.

1.8. УСТОЙЧИВОСТЬ НЕАВТОНОМНЫХ СИСТЕМ Рассмотрим неавтономную систему y = f (t, y ), (1.12) где y (t ) и f (t, y ) – N -мерные вектор-функции, причем f (t, y ) дваж ды непрерывно дифференцируемая функция в некоторой области : || y || r, t 0, и при этом выполняется соотношение f (t, 0) 0, 0 t, (1.13) 1.8. Устойчивость неавтономных систем т. е. система (1.12) имеет решение y (t ) 0. В этом случае говорят об устойчивости по Ляпунову или об асимптотической устойчивости ну левого решения y (t ) 0. Эти понятия формулируются аналогично со ответствующим понятиям для автономных систем.

Пусть x = (t ) – решение некоторой автономной системы x = g ( x). (1.14) Сделаем замену переменных, полагая x(t ) = (t ) + y (t ). В результате для переменной y (t ) получим систему уравнений вида y = g ( (t ) + y (t ) ) g ( (t ) ) f ( t, y (t ) ). (1.15) Очевидно, что имеет место тождество f (t, 0) 0 при 0 t, т. е.

условие (1.13) выполняется. Решение (t ) системы (1.14) называется устойчивым по Ляпунову или асимптотически устойчивым, если тако вым является нулевое решение y (t ) 0 системы (1.15). Асимптотиче (t ) ская устойчивость решения означает, что если величина || x(0) (0) || достаточно мала, то имеет место соотношение lim || x(t ) (t ) || = 0. Проведем линеаризацию системы (1.15) в окре t + стности решения (t ). Считая || y (t ) || малой величиной, получим g ( (t ) + y (t ) ) g ( (t ) ) = g ( (t ) ) y (t ) + O (|| y (t ) ||2 ), где g ( (t ) ) – матрица Якоби системы (1.14). Линеаризованная систе ма уравнений имеет вид y = g ( (t ) ) y. (1.16) Это линейная система уравнений, но с переменными коэффициентами.

Система (1.16) называется первой вариацией системы (1.14) относи тельно решения (t ).

Г л а в а 1. НЕКОТОРЫЕ СВЕДЕНИЯ ИЗ ТЕОРИИ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ Определение 1.8. Функция V (t, x) называется функцией Ляпуно ва системы (1.12), если:

1) V (t, x) определена и непрерывно дифференцируема при x, t 0 ;

2) V (t, 0) 0 при t 0 ;

3) существует положительно определенная в области функ ция W ( x) такая, что V (t, x) W ( x) при всех x, t 0 ;

4) V (t, x) 0 при всех x, t 0.

Итак, теоремы Ляпунова об устойчивости положения равновесия и об асимптотической устойчивости полностью переносятся на случай неавтономной системы.

Глава КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ Д анная глава посвящена вопросам нахождения приближенно го решения задачи Коши одношаговыми численными мето дами. Одной из основных проблем при численном интегрировании яв ляется вопрос о точности получаемых численных результатов. В по следнее время при получении неравенства для контроля точности при меняется в основном оценка главного члена глобальной ошибки, вычисленная тем или иным способом. По сравнению с использованием локальной ошибки это повышает надежность расчетов.

Ниже предложены два способа оценки аналога глобальной ошибки, которые не приводят к значительному увеличению вычислительных затрат. В первом способе требуется вычисление матрицы Якоби сис темы дифференциальных уравнений, поэтому оценку имеет смысл ис пользовать для контроля точности неявных или полуявных методов, где данная матрица участвует в вычислительном процессе. Второй способ не нуждается в вычислении матрицы Якоби, такую оценку можно применять в алгоритмах интегрирования на основе явных чис ленных схем.

При применении явных методов для решения жестких задач шаг интегрирования ограничен в основном устойчивостью численной схе мы. Поэтому соответствующие алгоритмы на основе явных формул обладают малой эффективностью. Здесь предложен способ получения неравенства для контроля устойчивости методов с ограниченной об ластью устойчивости, что повышает надежность расчетов. Данное не равенство основано на оценке максимального собственного числа мат рицы Якоби, которое определяется с использованием ранее вычисленных стадий.

Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ 2.1. ОСНОВНЫЕ ОПРЕДЕЛЕНИЯ Ниже будет рассматриваться задача Коши для неавтономной сис темы обыкновенных дифференциальных уравнений y = f (t, y ), y (t0 ) = y0, t0 t tk, (2.1) где y и f – вещественные N -мерные вектор-функции;

t – независи мая переменная, изменяющаяся на заданном интервале [t0, tk ]. Пред положение о гладкости правой части системы (2.1) влечет выполнение условий локальной теоремы Коши–Пикара, откуда, в свою очередь, следует существование единственного решения задачи (2.1);

[t0, tk ] – область определения решения. Поэтому всюду ниже это условие особо оговариваться не будет.

Введением дополнительной переменной систему (2.1) можно при вести к автономному виду. Определив x = (t, yT )T, g = (1, f T )T, x0 = (t0, y0 )T и перейдя снова к обозначениям y, f и y0, получим T y = f ( y ), y (t0 ) = y0, t0 t tk. (2.2) Поэтому часть рассуждений проведем для автономной системы. Слу чаи, когда различие задач (2.1) и (2.2) принципиально, будут оговари ваться отдельно. Ниже будем предполагать, что задача (2.2) является жесткой.

Определение 2.1. Задача Коши (2.2) называется жесткой на не котором интервале I [t0, tk ], если для t I имеет место Re(i ) 0, 1 i N ;

1) 2) max1i N Re(i ) / min1i N Re( i ) 1, где i есть собственные числа матрицы Якоби, вычисленной на ре шении y (t ).

Класс задач, называемых жесткими, весьма разнообразен, и поэто му в литературе встречаются другие определения жесткости, разли чающиеся степенью строгости. Поясним кратко, в чем заключается суть жесткости. Отметим, что речь идет в основном об устойчивых задачах, что следует, по крайней мере для автономных систем, из пер вой части определения 2.1. В некоторых определениях допускается наличие небольших по величине положительных собственных чисел 2.1. Основные определения матрицы Якоби и это достаточно существенно расширяет класс решаемых практических задач.

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

На других участках, называемых участками установления, произ водные решения относительно малы, а их длина велика по сравнению с длиной переходного участка. Отметим, что переходные участки могут отсутствовать вовсе. Однако в случае жесткой задачи даже небольшие возмущения (например, за счет ошибок округлений) приводят к их по явлению.

Трудности решения жестких задач классическими явными метода ми связаны с наличием участков установления. Здесь следует отме тить, что если весь интервал интегрирования представляет собой пере ходный участок, то такая задача в ряде случаев к жестким не относится, хотя в смысле определения 2.1 она может быть жесткой.

Так как на переходном участке производные от решения велики, шаг интегрирования из условия точности выбирается небольшим. Можно сказать, что на данном участке выполняется неравенство hL(t ) C, (2.3) где h – шаг интегрирования, L(t ) – классическая константа Липшица функции f ( y ), C – постоянная. Поэтому решение может быть найде но явным методом даже с небольшой областью устойчивости. До не давнего времени полагалось, что явные методы можно использовать при C ( C 10 ). В связи с развитием явных методов с расширенными областями устойчивости их возможности определяются более широко.

Ниже мы еще вернемся к этому вопросу.

На промежутке установления производные от решения невелики, поэтому шаг по точности может быть выбран большим. Можно ска зать, что характерным для данного участка является выполнение нера венства hL(t ) C. (2.4) Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ Потому применение классических явных методов, для которых усло вие (2.3) необходимо для устойчивости, практически невозможно на современных ЭВМ.

Итак, ниже под жесткими задачами будем понимать такие задачи, которые являются жесткими в смысле определения 2.1 и которые удовлетворяют условиям (2.3) и (2.4) на интервале интегрирования.

Отметим, что на практике жесткие задачи распознаются хорошо. Лю бая физическая система, которая моделируется системой обыкновен ных дифференциальных уравнений и имеет компоненты с сильно раз личающимися временными константами, приводит к жесткой задаче.

Далее будут рассматриваться одношаговые безытерационные ме тоды, которые в векторной форме можно записать в виде yn+1 = yn + h f (tn, yn, h), n = 0, 1, 2…, (2.5) с заданным начальным условием y0. Здесь h – шаг интегрирования, f – заданная вещественная N -мерная вектор-функция, зависящая от правой части исходной задачи. В форме (2.5) можно представить яв ные методы типа Рунге–Кутты, а также неявные и полуявные методы типа Рунге–Кутты с фиксированным итерационным процессом, в ко тором число итераций не зависит от номера шага интегрирования.

Рассмотрим источники погрешностей, которые могут возникать при реализации методов (2.5). Это ошибки округлений и ошибки, свя занные с применением приближенных формул интегрирования. Одной из основных проблем при построении эффективных алгоритмов чис ленного решения задачи Коши (2.2) является получение надежных оценок полной погрешности метода. Пусть приближенное решение y0, y1, …, yM, Mh = tk t0, задачи (2.2) вычислено по формуле (2.5), y (t ) – точное решение (2.2).

Определение 2.2. Величина y (tn ) yn называется полной погреш ностью метода (2.5) в точке t = tn.

Величина ошибок округлений в сильной степени зависит от кон кретного вычислительного устройства. Особое внимание ошибкам ок руглений следует уделять при наличии в алгоритме интегрирования итерационных процессов. Обычно о сходимости судят по близости двух соседних итераций, разность между которыми невозможно свести к нулю из-за наличия ошибок округлений. В этом случае необходимо 2.1. Основные определения выбрать константу, ограничивающую количество итераций. С другой стороны, ошибки округлений можно трактовать как возмущение пра вой части задачи (2.2). Если исходная задача обладает малой чувстви тельностью к малым возмущениям, то и от метода ее решения естест венно требовать того же свойства. Такие методы называют устойчивыми. В дальнейшем будут рассматриваться только устойчи вые методы, в которых не предполагается использование итерацион ных процессов. Более того, ниже будут изучаться в основном методы относительно невысокого порядка точности, для которых ошибки за счет неточности приближенной формулы существенно превышают ошибки округлений, возникающие при реализации (2.5). Поэтому вме сто полной погрешности метода будет рассматриваться только ошибка за счет неточности (2.5), называемая глобальной ошибкой.

Определение 2.3. Метод из класса (2.5) сходится, если для каж дой задачи (2.2) имеет место max 0n M || y (tn ) yn || 0 при h 0, и сxодится с порядком p, если max 0n M || y (tn ) yn ||= O(h p ), где Mh = tk t0, || || – некоторая норма в R N.

Очевидно, что для решения (2.2) имеет смысл рассматривать толь ко сходящиеся методы (2.5). При исследовании сходимости (2.5) тре буется определить ошибку n = y (tn ) yn, что в ряде случаев является достаточно сложной задачей. Поэтому в качестве предварительной оценки точности формулы (2.5) обычно используют погрешность ап проксимации.

Определение 2.4. Погрешностью аппроксимации n+1 схемы (2.5) в точке tn+1 [t0, tk ] называется величина, определяемая по формуле n +1 = y (tn +1 ) y (tn ) h f (tn, y (tn ), h), (2.6) где y (t ) есть точное решение задачи (2.2).

Если yn = y (tn ), то погрешность аппроксимации n+1 = = y (tn +1 ) yn +1 в точке tn+1 совпадает с величиной глобальной ошиб ки, т. е. погрешность аппроксимации – ошибка, которая получается за один шаг интегрирования, причем она возникает за счет отбрасыва ния членов при конечной аппроксимации производной дифференци ального уравнения. Поэтому в ряде работ погрешность аппроксимации Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ называют локальной ошибкой. Ниже нам также будет удобно пользо ваться этим термином.

Определение 2.5. Говорят, что метод (2.5) аппроксимирует систему (2.2), если h 1 max 0n M || n || 0 при h 0, и имеет max 0n M || n ||= O(h p +1 ), порядок аппроксимации p, если где N Mh = tk t0, || || – некоторая норма в R.

Метод, удовлетворяющий определению 2.5, еще называют согласо ванным с порядком p. Отметим, что в ряде публикаций условия схо димости и аппроксимации формулируются не при h 0, а при n, что одно и то же. При мельчении шага число узлов сетки возрастает, а при увеличении числа узлов величина шага стремится к нулю.

Теперь перейдем к изучению устойчивости численных формул (2.5). Для этого рассмотрим вторую пару аппроксимаций y n и y n +1, которые удовлетворяют соотношению y n +1 = y n + h f (tn, y n, h).

Определение 2.6. Метод (2.5) называется С-устойчивым на клас се задач, если существуют такие вещественные числа h0 0 и C, что для всех задач из данного класса выполняется неравенство || yn+1 y n +1 || || yn y n ||, = 1 + Ch, h (0, h0 ).

Устойчивость в смысле определения 2.6 называют еще устойчиво стью сходимости. Если 1, то метод называется контрактивным. Из вестно, что любой согласованный метод вида (2.5) является С-устойчивым. Отсюда следует теорема о том, что одношаговый метод из класса (2.5) сходится тогда и только тогда, когда он является согла сованным. Поэтому для доказательства сходимости методов вида (2.5) достаточно показать, что они аппроксимируют задачу (2.2) с порядком p, p 1. Из определения 2.6 следует, что при C 0 и для одного и того же шага интегрирования для С-устойчивых методов разность ме жду двумя любыми численными решениями, определенными по фор муле (2.5), остается ограниченной величиной. Очевидно, что при прак тической реализации методов на ЭВМ этого свойства может оказаться недостаточно, потому что соответствующие алгоритмы интегрирова ния могут быть малоэффективными. Поэтому были рассмотрены дру гие свойства устойчивости, такие как абсолютная устойчивость, 2.1. Основные определения A -устойчивость, L -устойчивость и др. Данные свойства называют вычислительной устойчивостью.

Рассмотрим понятие абсолютной устойчивости. Оно вводится на линейном скалярном уравнении y = y, y (0) = y0, t 0, (2.7) с комплексным, Re() 0. Ниже будем предполагать, что в форму ле (2.5) функция f линейна по y, если исходная система (2.1) ли нейная. Тогда, применяя (2.5) для решения (2.7), получим yn+1 = Q ( z ) yn, z = h, (2.8) где Q ( z ) называют функцией роста или функцией устойчивости мето да (2.5). Для рассматриваемых ниже методов Q ( z ) является либо мно гочленом относительно z, либо рациональной функцией.

Определение 2.7. Метод (2.5) называется абсолютно устойчивым для данного значения z = h, если | Q( z ) | 1. Область R комплексной плоскости z называется областью абсолютной устойчивости мето да (2.5), если он устойчив при всех z R. Пересечение области устой чивости с вещественной осью называется интервалом устойчивости.

Абсолютная устойчивость является естественным требованием, ес ли выполняется неравенство Re( z ) 0, поскольку модуль точного ре шения y (t ) = exp(t ) y0 задачи (2.7) есть невозрастающая функция. Не смотря на простой вид уравнения (2.7), оно успешно и достаточно долго используется в качестве модельного для предсказания устойчи вости практических численных методов решения нелинейных систем общего вида. Аргументация в пользу уравнения (2.7) заключается в следующем. Пусть в некоторый момент времени t = t1 решение y (t ) задачи (2.1) возмущено до некоторого значения y (t ). Тогда при t t имеет место x = A(t ) x, где x(t ) = y (t ) y (t ), а матрица A(t ) вычислена в соответствии с теоремой о среднем значении для векторных функций A(t ) = f (t, y (t ) + (1 ) y (t ))d, f (t, y (t )) – матрица Якоби системы (2.1). Теперь жесткость задачи (2.1) может быть полностью описана в терминах жесткости получен Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ ной линейной системы с переменными коэффициентами. Этот прием называется глобальной линеаризацией, поскольку линеаризация осу ществляется вдоль всей кривой y (t ) в целом. Следующий шаг – ло кальная линеаризация, которая сводится к замораживанию матрицы A(t ) в точке t1, т. е. A = A(t1 ) = f (t1, y (t1 )). Переобозначая x через y, получим линейную систему с постоянными коэффициентами, т. е.

y = Ay, (2.9) и теперь жесткость задачи (2.1) описывается в терминах линейной сис темы (2.9) с постоянными коэффициентами. Недостатки использова ния локальной линеаризации связаны с тем, что спектр замороженной матрицы не всегда дает правильную информацию о качественном по ведении решения задачи (2.1) – есть соответствующие примеры. Одна ко для многих практических задач переход к системе (2.9) вполне оп равдан.

Применяя (2.5) для решения задачи (2.9), получим yn+1 = Q(hA) yn, где, в отличие от (2.8), Q (hA) есть матрица. Матричная функция Q(hA) существует, если комплекснозначная функция Q ( z ) определена на спектре матрицы hA, т. е. Q(h1 ), …, Q(h N ) – собственные числа матрицы Q (hA). Если в задаче (2.7) интерпретировать как произ вольное собственное число матрицы A, то становится понятным ис пользование (2.7) для прогноза устойчивости методов (2.5). Ниже бу дем в основном пользоваться понятием абсолютной устойчивости, которую для сокращения назовем просто устойчивостью.

Основная проблема при решении жестких задач явными методами – проблема численной устойчивости. Поясним этот факт на примере яв ного метода Эйлера yn+1 = yn + hf (tn, yn ). Применяя его для решения задачи (2.7), получим yn+1 = (1 + h) yn, откуда условие устойчивости имеет вид | 1 + h | 1 или h 2/ | |. При больших значениях Re() и достаточно большом промежутке интегрирования это условие есть весьма обременительное ограничение на размер шага интегрирования.

Поэтому были рассмотрены A -устойчивые методы.

Определение 2.8. Численный метод называется A -устойчивым, если его область устойчивости включает всю полуплоскость Re(h) 0.

2.2. Контроль точности вычислений Свойство A -устойчивости гарантирует, что | Q(h) | 1 для всех Re(h) 0. Однако для многих одношаговых A -устойчивых методов Q(h) таково, что | Q(h) | 1 при Re(h). В результате чис ленные приближения к быстрозатухающим фундаментальным реше ниям с малыми временными постоянными могут затухать очень мед ленно. Следствием этого является понижение эффективности алгоритма интегрирования. Поэтому было введено понятие L устойчивости.

Определение 2.9. Одношаговый метод называется L -устой чивым, если он A -устойчив и если | Q(h) | 0 при Re(h).

Современные эффективные алгоритмы интегрирования жестких задач основаны на L -устойчивых численных схемах. Приведем в каче стве примера неявный метод Эйлера, который имеет вид yn+1 = yn + hf (tn +1, yn+1 ). Применяя его для решения задачи (2.7), получим yn +1 = Q(h) yn, Q(h) = 1 / (1 h), откуда следует L -устой чивость неявного метода Эйлера.



Pages:   || 2 | 3 | 4 | 5 |   ...   | 9 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.