РЕДАКЦИОННАЯ КОЛЛЕГИЯ СЕРИИ «МОНОГРАФИИ НГТУ»



д-р техн, наук, проф. (председатель) Н В. Пустовой
д-р техн. наук, проф. (зам. председателя) А.Г Вострецов д-р техн. наук, проф. (отв. секретарь) В.Н. Васюков

д-р техн. наук, проф. А.А. Батаев
д-р техн. наук, проф. А. С. Востриков
д-р техн. наук, проф. А.А. Воевода
д-р техн. наук, проф. В.В. Губарев
д-р техн. наук, проф. В.И Денисов
д-р физ.-мат. наук, проф. А.К Дмитриев
д-р физ.-мат. наук, проф. В.Г Дубровский
д-р филос. наук, проф. В.И. Игнатьев
д-р филос. наук, проф. В.В. Крюков
д-р техн. наук, проф. ГИ Расторгуев
д-р физ.-мат. наук, проф. В.А. Селезнев
д-р техн. наук, проф. Ю.Г Соловейчик
д-р техн. наук, проф. А.А. Спектор
д-р экон. наук, проф. В. А. Титова
д-р юр. наук, доц. В.Л Толстых
д-р техн. наук, проф. А.Г. Фшиов
д-р техн. наук, проф. А. Ф. Шевченко

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







КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ЖЕСТКИХ ГИБРИДНЫХ СИСТЕМ













НОВОСИБИРСК
2013

УДК 004.9
     Н731



Рецензенты:
Заслуженный деятель науки РФ, д-р техн. наук, профессор В.И. Денисов; д-р физ.-мат. наук, гл. науч. сотр. ИВТ СО РАН Л.Б. Чубаров





      Новиков Е.А.
Н731 Компьютерное моделирование жестких гибридных систем: монография / Е.А. Новиков, Ю.В. Шорников. - Новосибирск: Изд-во НГТУ, 2013. - 451 с. (Серия «Монографии НГТУ»).

         ISBN 978-5-7782-2023-2

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




УДК 004.9





ISBN 978-5-7782-2023-2

             © Новиков Е.А., Шорников Ю.В., 2012, 2013
             © Новосибирский государственный технический университет, 2012, 2013

Ministry ofEducation and Science of the Russian Federation
NOVOSIBIRSK STATE TECHNICAL UNIVERSITY





E.A. Novikov, Yu.V. Shornikov







COMPUTER SIMULATION OF STIFF HYBRID SYSTEMS

Monograph











NOVOSIBIRSK
2013

UDC 004.9
     N731



Reviewers:
Prof. V.I. Denisov, D. Sc. (Eng.), RF Honoured Scientists; L.B. Chubarov, D. Sc. (Phys. & Math.), ChiefScientist, SB RAS Computing Technology Institute







      Novikov E.A.
N731 Computer Simulation of Stiff Hybrid Systems: monograph / E.A. Novikov, Yu.V. Shornikov. - Novosibirsk : NSTU Publishers, 2013.- 451 p. («NSTU Monographs» series).


         ISBN 978-57782-2023-2


             The monograph is concerned with designing unconventional numerical methods to solve the Cauchy problem for stiff systems of ordinary differential equations. Special 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 algorithms. The methodology of hybrid systems is described in detail and their classification is given in the monograph. The capabilities of instrumental environment of hybrid 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 physical, chemical, biological and other processes.


UDC 004.9









                                        © Novikov E.A., Shornikov Yu.V., 2012, 2013
                                        © Novosibirsk State Technical

                                                University, 2012, 2013


ISBN 978-57782-2023-2

        ВВЕДЕНИЕ




        о многих приложениях возникает проблема численного решения задачи Коши для жестких систем обыкновенных диф-

ференциальных уравнений. Основные тенденции при построении чис

ленных методов связаны с расширением их возможностей для решения жестких систем все более высокой размерности [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]. Он заключается в том, что одношаговой формуле р-го порядка точности в соответствие ставится многошаговая схема (р + 1)-го порядка. Затем данная схема преобразуется таким образом, чтобы после подстановки в нее приближений получилась оценка локальной ошибки 8 „ _ одношагового метода. , р
   Недостатком данного способа является многошаговость оценки, что приводит ко всем недостаткам многошаговых методов.
   В последнее время все чаще оценку локальной ошибки вычисляют с помощью вложенных методов. Приближение к решению в каждой точке вычисляется двумя методами вида р-го и (р + 1)-го порядков точности. Затем локальная ошибка метода р-го порядка оценивается через разность полученных приближений к решению. Обычно такой способ используется тогда, когда для расчетов по методу р-го порядка не требуется дополнительных вычислений правой части и матрицы Якоби. Следует отметить оперативность и относительную дешевизну оценки локальной ошибки с помощью вложенных методов, по затра

ВВЕДЕНИЕ

7

там на шаг лежащей между оценкой ошибки с помощью правила Рунге и многошаговой оценкой. В то же время по отношению к многошаговой оценке в ней при вычислении ошибки используется информация только с данного шага, что повышает ее надежность. Данный способ успешно применялся в работах [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], в которой было введено понятие А -устойчивости. Однако требование оказалось слишком обременительным для линейных многошаговых методов, и поэтому были введены менее ограничительные определения устойчивости [163, 176, 177, 272]. Среди многошаговых методов наибольшее распространение получили формулы дифференцирования назад [271, 272, 274, 275, 286, 287], обладающие свойством жесткой устойчивости [163].
   Понятие А -устойчивости привело к рассмотрению неявных методов типа Рунге-Кутты. Наиболее полное исследование этих методов освещалось в работах [229-230, 232-237, 301], а позднее - в монографиях [176, 177, 207, 237]. В [234, 235] доказана теорема о том, что для каждого т существует неявная т -стадийная схема, порядок точности которой равен 2т . Заметим, что аналогичная теорема для явных методов типа Рунге-Кутты отсутствует. Согласно [207] функция роста или функция устойчивости схемы максимального порядка является диагональной аппроксимацией Падэ для функции ехр(х). Если отказаться от максимального порядка, то можно построить схемы с лучшими свойствами устойчивости [176, 177].
   Несмотря на хорошие свойства точности и устойчивости, практическое использование неявных методов типа Рунге-Кутты ограничено

ВВЕДЕНИЕ

9

вследствие больших вычислительных затрат на шаг интегрирования. При реализации требуется итерационный процесс. Метод простой итерации не эффективен при решении жестких задач, потому что он приводит фактически к такому же ограничению на размер шага, что и явный метод. Поэтому возникает необходимость использования метода Ньютона или каких-либо его модификаций. Это приводит к необходимости обращения матрицы размерности (т ■ N х т ■ N), где т есть число стадий, a N - размерность системы. Значительно сократить вычислительные затраты при реализации можно за счет использования одной и той же матрицы на нескольких шагах интегрирования. Замораживание матрицы становится возможным благодаря тому, что это не влияет на порядок точности численной схемы, а определяет только скорость сходимости итерационного процесса. Поэтому необходимость в ее пересчете возникает при значительном замедлении сходимости итерационного процесса.
    Трудности с реализацией неявных методов типа Рунге-Кутты, привели к поискам более простых их модификаций. В работах [78, 163, 176, 177, 237] был рассмотрен класс полуявных формул типа Рунге-Кутты, т. е. таких методов, для которых имеет место Ргу = 0 при i < j. В этом случае итерационная матрица является блочно-диагональной, причем число блоков совпадает с числом стадий, а размерность каждого блока - с размерностью вектора решения. В результате вместо обращения матрицы размерности (т ■ N х т ■ N) надо обратить т матриц размерности N каждая. Исследование полуявных методов содержится в работах [163, 176, 177]. Заметим, что если матрица, составленная из коэффициентов Р, 1 < i, j < т , имеет одно т-кратное собственное число, то неявный метод можно реализовать с теми же затратами, что и полуявный метод [163]. Вопросам построения диагонально-неявных методов посвящена работа [78]. Дальнейшего сокращения вычислительных затрат можно добиться, если положить равными все Рц, 1 < i < т , и аппроксимировать все диагональные матрицы одной. Тогда на шаге требуется обратить только одну матрицу размерности N. Для этого случая в [207] доказана теорема о том, что порядок точности (т + 2) не может быть достигнут ни для какого т-стадийного полуяв-ного метода при рп =... = Ртт.