WWW.DISS.SELUK.RU

БЕСПЛАТНАЯ ЭЛЕКТРОННАЯ БИБЛИОТЕКА
(Авторефераты, диссертации, методички, учебные программы, монографии)

 


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

«КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ЖЕСТКИХ ГИБРИДНЫХ СИСТЕМ Е.А. Новиков, Ю.В. Шорников КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ЖЕСТКИХ ГИБРИДНЫХ СИСТЕМ НОВОСИБИРСК 2012 УДК 004.9 Н 731 Рецензенты: Заслуженный ...»

-- [ Страница 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. 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.

ISBN 978-57782-2023-5 © Novikov E.A., Shornikov Yu.V.,

ВВЕДЕНИЕ

В о многих приложениях возникает проблема численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений. Основные тенденции при построении численных методов связаны с расширением их возможностей для решения жестких систем все более высокой размерности [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] на основе схемы типа Розенброка второго порядка точности построен алгоритм интегрирования с замораживанием матрицы Якоби, уже при построении такого алгоритма на основе формулы третьего порядка возникают принципиальные трудности. В работах средственно в вычислительные формулы, приведены для решения дифференциально-алгебраических задач.

В [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].

Существует набор контекстов, в которых область определения правой части системы дифференциальных уравнений ограниченна. В одних задачах векторное поле определено не везде, поскольку модель спроектирована для применения только в некоторых областях фазового пространства. Например, нельзя применять сверхзвуковую аэродинамическую модель летательного аппарата на дозвуковых скоростях вием физической природы проблемы. Так, в задачах химической кинетики фазовые переменные соответствуют концентрациям различных компонентов. Их отрицательные значения не имеют физического смысла, но могут быть получены вследствие численных ошибок моделирования или неадекватности модели. Подобные ситуации возникают в проблемах с фазовыми переходами. Такие модели характеризуются тем, что имеют области неопределенности в фазовом пространстве. В этом случае управляющий компьютер должен переключиться на другой закон управления в этой области. Пример – контроллер, который вырабатывает управляющее решение в зависимости от результатов расчета инверсной кинематики для руки робота [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).

НЕКОТОРЫЕ СВЕДЕНИЯ ИЗ ТЕОРИИ

ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

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

Рассмотрим систему обыкновенных дифференциальных уравнений первого порядка вида где 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 ), что имеет место где y0 = ( y1, y2,…, y N )T – заданный вектор, верхний индекс T означает транспонирование. Вектор 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 ) удовлетворяет условию Липшица, т. е.

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

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

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

Замечание. Условие Липшица выполнено, если все компоненты вектор-функции f (t, y ) имеют в области E ограниченные частные | fi (t, y1, y2, …, y N ) / y j | C, 1 i, j N, где C – некоторое положительное число.

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

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

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

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

ОТ НАЧАЛЬНЫХ ДАННЫХ

Рассмотрим задачу Коши для системы дифференциальных уравнений, в которой правая часть и начальные данные зависят от параметров z, z = ( z1, z2,…, zm )T. Такая задача имеет вид Теорема о непрерывной зависимости и дифференцируемости решения по параметру. Пусть в некоторой области G пространства (t, y, z ) вектор-функция f имеет непрерывные производные по t и y порядка p, p 1. Тогда решение y (t, z ) задачи Коши (1.3) непрерывно и имеет p непрерывных производных по t и z.

Теорема о непрерывной зависимости и дифференцируемости решения по начальным данным. Пусть в задаче Коши

Г л а в а 1. НЕКОТОРЫЕ СВЕДЕНИЯ ИЗ ТЕОРИИ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

вектор-функция f имеет непрерывные производные порядка p, p по t и y. Тогда решение y (t, t0, y0 ) непрерывно и имеет p непрерывных производных по t, t0 и y0.

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

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

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

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

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

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

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

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

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

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

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

Устойчивость системы (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 Определение 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 асимптотически устойчиво.

РАВНОВЕСИЯ ЛИНЕЙНОЙ СИСТЕМЫ

Рассмотрим линейную однородную систему с постоянными комплексными коэффициентами где 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, – собственные числа матрицы A, а число 0 можно выбрать сколь угодно малым.

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

1.7. Устойчивость положений равновесия нелинейных систем Критерий Рауса–Гурвица. Пусть все коэффициенты уравнения (1.8) вещественны и a0 0. Для того чтобы все корни уравнения (1.8) имели отрицательные вещественные части, необходимо и достаточно, чтобы все главные миноры определителя были положительны.

РАВНОВЕСИЯ НЕЛИНЕЙНЫХ СИСТЕМ

Исследование устойчивости нелинейной системы будем сводить к исследованию устойчивости линейной системы. Для этого применим метод линеаризации. Пусть a – положение равновесия автономной системы где y (t ) и f ( y ) – N -мерные вектор-функции, причем f ( y ) дважды непрерывно дифференцируема в некоторой окрестности U точки a.

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

где 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 = 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 (t ) и f (t, y ) – N -мерные вектор-функции, причем f (t, y ) дважды непрерывно дифференцируемая функция в некоторой области : || y || r, t 0, и при этом выполняется соотношение т. е. система (1.12) имеет решение y (t ) 0. В этом случае говорят об устойчивости по Ляпунову или об асимптотической устойчивости нулевого решения y (t ) 0. Эти понятия формулируются аналогично соответствующим понятиям для автономных систем.

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

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

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

Г л а в а 1. НЕКОТОРЫЕ СВЕДЕНИЯ ИЗ ТЕОРИИ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Определение 1.8. Функция V (t, x) называется функцией Ляпунова системы (1.12), если:

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

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

КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ

ОДНОШАГОВЫХ МЕТОДОВ

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

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

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

Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ

Ниже будет рассматриваться задача Коши для неавтономной системы обыкновенных дифференциальных уравнений где 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, получим Поэтому часть рассуждений проведем для автономной системы. Случаи, когда различие задач (2.1) и (2.2) принципиально, будут оговариваться отдельно. Ниже будем предполагать, что задача (2.2) является жесткой.

Определение 2.1. Задача Коши (2.2) называется жесткой на некотором интервале I [t0, tk ], если для t I имеет место где i есть собственные числа матрицы Якоби, вычисленной на решении y (t ).

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

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

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

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

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

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

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

Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ

Потому применение классических явных методов, для которых условие (2.3) необходимо для устойчивости, практически невозможно на современных ЭВМ.

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

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

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

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

и теперь жесткость задачи (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.

Свойство 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 -устойчивость неявного метода Эйлера.

2.2. КОНТРОЛЬ ТОЧНОСТИ ВЫЧИСЛЕНИЙ

Здесь и ниже будут рассматриваться оценки локальной и глобальной ошибок в смысле главного члена, т. е. далее будет оцениваться только первый член при разложении соответствующих ошибок в ряды Тейлора по степеням h. Итак, пусть для решения задачи (2.2) имеется метод вида (2.5) p -го порядка точности и пусть его локальная ошибка n, p представима в виде где – некоторая функция, не зависящая от размера шага интегрирования.

Для записи (2.10) требуется определенная гладкость правой части исходной задачи на промежутке [tn, tn +1 ], которая предполагается выполненной. Известно, что при условии (2.10) глобальную ошибку n, p можно вычислить по формуле

Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ

где x(t ) есть решение задачи Коши f = f / y – матрица Якоби системы дифференциальных уравнений.

Непосредственное использование оценки (2.11) в неравенстве для контроля точности вычислений может приводить к значительному увеличению вычислительных затрат. Дело в том, что для определения n, p требуется оценка локальной ошибки (2.10), что означает, как правило, дополнительные вычисления правой части исходной задачи.

Кроме того, решение задачи (2.12) связано с вычислением N 2 компонент матрицы Якоби системы (2.2) плюс затраты на интегрирование (2.12). Поэтому такой способ используется лишь иногда при решении жестких задач, когда матрица Якоби уже известна.

Обычно при реализации A -устойчивых или L -устойчивых методов вида (2.5) имеется матрица Dn вида где E – единичная матрица, a – известное число.

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

Теорема 2.1. Пусть для решения задачи (2.1) имеется метод p -го порядка точности и пусть его локальная ошибка n, p представима в виде где (t ) есть некоторая функция, не зависящая от размера шага интегрирования. Тогда для любого номера шага n имеет место соотношение где || || – некоторая норма в R N, x(t ) – решение задачи Коши (2.12), а xn вычисляется по следующей рекуррентной формуле:

Для доказательства подставим (t ) = f (t, y (t )) (t ) в (2.12), т. е. для определения x(t ) получим следующую задачу:

Численное решение (2.17) определим с помощью схемы типа Розенброка первого порядка точности, которая применительно к задаче (2.1) имеет вид где величина a при задании Dn по формуле (2.13) определена при описании исходного метода p -го порядка точности. Применяя (2.18) для решения (2.17), запишем a ( n n + xn xn ) и группируя соответствующим образом слагаемые, получим Dn ( xn+1 xn ) = a 1[( xn n ) + Dn ( n xn )]. Применяя слева к полученному соотношению матрицу Dn 1, получим формулу (2.16). Равенство (2.15) следует из сходимости схемы (2.18) с первым порядком. Теорема доказана.

Заметим, что схема (2.18) при a 0,5 является A -устойчивой, а при a = 1 L -устойчивой. Это означает, что при значениях a, 0 a 0,5, оценка глобальной ошибки, вычисленная с применением (2.16), может возрастать, в то время как истинная ошибка будет убывать. Однако для многих L -устойчивых методов (2.5) имеет место 0 a 0,5, и поэтому для них использование (2.16) при выборе величины шага интегрирования приведет к понижению эффективности

Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ

алгоритма интегрирования. С целью устранения этого недостатка рассмотрим (2, 1)-метод, который применительно к (2.1) имеет вид Нетрудно убедиться, что при 0,5 0, 25 2 a 0,5 + 0, 25 2 данная схема имеет первый порядок точности и является L -устойчивой.

Теорема 2.2. Пусть выполнены условия теоремы 2.1. Тогда имеет место соотношение (2.15), где Для доказательства применим численную схему (2.20) для решения задачи (2.17). Получим Добавим в правую часть соотношение ( n xn n + xn ). Учитывая вид матрицы Dn, сгруппируем определенным образом слагаемые и умножим слева полученное соотношение на Dn. В результате запишем Теперь добавим к выражению в квадратных скобках соотношение ( n xn n + xn ) и снова сгруппируем слагаемые, учитывая вид Dn.

Получим Умножив слева полученное уравнение на матрицу Dn 2, имеем соотношение (2.21). Равенство (2.15) следует из сходимости (2.20).

Теорема доказана.

Теперь рассмотрим вопрос о том, в каком случае применение (2.16) лительных затрат, чем непосредственное интегрирование (2.12). Прежде всего отметим, что при вычислении стадий A -устойчивых или L устойчивых методов (2.5) будет использоваться LU -разложение матрицы Dn с выбором главного элемента по строке или столбцу, а иногда по всей матрице, т. е. при вычислении приближенного решения по формуле (2.5) осуществляется декомпозиция матрицы Dn. Поэтому вычисление оценки ошибки по формуле (2.16) не приводит к значительному увеличению вычислительных затрат, которые практически определяются числом операций на выполнение обратного хода в методе Гаусса. Далее, как уже отмечалось выше, использование (2.12) связано с необходимостью оценки локальной ошибки, что, как правило, приводит к дополнительным вычислениям правой части задачи (2.2) и как следствие – увеличению вычислительных затрат. При использовании (2.16) требуется оценить функцию (t ), что в ряде случаев удается осуществить бесплатно. Ниже при оценке (t ) будем поступать следующим образом.

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

где c p 1 и c p – некоторые вычисляемые постоянные, а функция p (t ) не зависит от размера шага интегрирования. Вычисления предполагается проводить по методу p -го порядка, для определения глобальной ошибки которого требуется оценка функции (t ) = c p p (t ).

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

Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ

могут быть использованы те коэффициенты, которые остаются свободными после удовлетворения требованиям аппроксимации. Класс методов (2.5) с условиями (2.22) не пуст. Теперь докажем утверждение, с использованием которого далее будет определяться оценка функции p (t ).

Теорема 2.3. Пусть для решения задачи (2.2) имеется два метода вида (2.5) ( p 1) -го и p -го порядков точности и пусть их локальные ошибки удовлетворяют требованиям (2.22). Тогда выполняется соотношение Для доказательства предположим, что приближение к решению yn в точке tn вычислено методом p -го порядка, т. е. выполнено равенство Тогда приближенное решение yn+1 в точке tn +1, вычисленное методом p -го порядка, можно представить в виде Разлагая функцию f (tn, y (tn ) + n, h) в ряд Тейлора в окрестности точки y (tn ) и используя очевидное равенство получим Добавляя разность y (tn+1 ) y (tn +1 ) в правую часть полученного соотношения и учитывая определение погрешности аппроксимации, будем иметь Проводя аналогичные рассуждения для метода ( p 1) -го порядка, имеем Учитывая, что n, p = O(h p +1 ), n, p 1 = O (h p ), и вычитая (2.27) из (2.26), получим откуда, с учетом (2.22), следует (2.23). Теорема доказана.

При использовании описанного выше способа оценки глобальной ошибки экономия вычислительных затрат достигается за счет того, что в ряде случаев при построении метода ( p 1) -го порядка точности используются те же самые стадии, что и в методе p -го порядка, т. е. метод ( p 1) -го порядка является побочным продуктом метода p -го порядка.

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

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

Возможны и другие способы оценки функции p (t ), не связанные с построением вспомогательного метода. При построении конкретных методов этот вопрос будет обсуждаться более подробно. Вычисление оценки глобальной ошибки n, p метода вида (2.5) p -го порядка точности по формуле где xn определяется из соотношения (2.16), связано с вычислением матрицы Якоби. Во многих случаях при решении задач умеренной жесткости можно применить явные методы, в которых матрица Якоби не

Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ

используется и для которых оценка ошибки по формуле (2.16) означает существенное увеличение вычислительных затрат. В то же время при интегрировании умеренно жестких задач явными методами основные затраты приходятся на участок установления, где шаг интегрирования ограничен не требованием точности, а устойчивостью численной схемы и где точность расчетов, как правило, значительно завышена. Как показывают расчеты, в этом случае в качестве оценки ошибки метода можно использовать величину c p h p p (t ). Докажем следующее утверждение.

Теорема 2.4. Пусть для решения скалярной задачи (2.1) имеется метод p -го порядка точности, локальная ошибка которого имеет вид (2.22). Тогда выполняется неравенство || x(t ) || | c p | || p (t ) ||, где решение задачи Коши Для доказательства введем обозначение F = знакопостоянства f (t, y (t )) для скалярного уравнения вида (2.1) с использованием теоремы о среднем можно записать где z0 [t0, tk ]. Теорема доказана.

Основываясь на этом соображении при N 1, будем предполагать выполнение неравенства max1i N || xi (t ) || | c p | max1i N || pi (t ) ||, т. е. в качестве оценки глобальной ошибки n, p используем иногда следующую оценку: n, p c p h p p (tn ). Учитывая теорему 2.3, эту оценку можно записать в виде Теперь для контроля точности вычислений и при выборе величины шага интегрирования можно применить неравенство где – требуемая точность интегрирования, || || – некоторая норма в R N. Во всех приведенных ниже алгоритмах интегрирования левая часть неравенства (2.31) определяется по формуле || n, p || = max1i N | i, p | /(| yn | + r ), где i – номер компоненты, r – положительный параметр. Если по i -й компоненте решения выполняi ется неравенство | yn | r, то контролируется абсолютная ошибка r, в противном случае – относительная ошибка.

Отметим еще одну важную особенность построенных оценок. Для этого рассмотрим L -устойчивый метод (2.5), который применим для решения задачи (2.7). Тогда в (2.8) в силу L -устойчивости схемы выполняется Q ( z ) 0 при z. Так как для точного решения y (tn+1 ) = exp(h) y (tn ) задачи (2.7) выполняется аналогичное свойство, естественно требование стремления к нулю глобальной ошибки при условии z. В силу того что все приведенные выше рассуждения относились к оценке главного члена ошибки, это требование может не выполняться. Поэтому с целью исправления асимптотического поведения глобальной ошибки вместо построенных выше оценок n, p будем рассматривать оценку n, p ( jn ) вида n, p ( jn ) = D1 jn n, p, 1 jn J f, где параметр J f зависит от асимптотических свойств n, p. Теперь для контроля точности и при выборе величины шага интегрирования можно проверять следующее неравенство:

Г л а в а 2. КОНТРОЛЬ ТОЧНОСТИ И УСТОЙЧИВОСТИ ОДНОШАГОВЫХ МЕТОДОВ



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


Похожие работы:

«Министерство образования и науки Российской Федерации Государственное образовательное учреждение высшего профессионального образования Ивановский государственный энергетический университет имени В.И. Ленина А.И. Тихонов Живая планета или поиск нового подхода к миропониманию Иваново 2011 ББК 20 Т46 Тихонов А.И. Живая планета или поиск нового подхода к миропониманию / ГОУВПО Ивановский государственный энергетический университет имени В.И. Ленина. – Иваново, 2011. – 84 с. ISBN В данной монографии...»

«Н.Ф. ГЛАДЫШЕВ, Т.В. ГЛАДЫШЕВА, С.И. ДВОРЕЦКИЙ, С.Б. ПУТИН, М.А. УЛЬЯНОВА, Ю.А. ФЕРАПОНТОВ РЕГЕНЕРАТИВНЫЕ ПРОДУКТЫ НОВОГО ПОКОЛЕНИЯ: ТЕХНОЛОГИЯ И АППАРАТУРНОЕ ОФОРМЛЕНИЕ Монография Москва Издательство Машиностроение-1 2007 УДК 661.183:546.32-39+546.41-36 ББК Л113.2 Р177 Рецензенты: Доктор химических наук, профессор Воронежского государственного университета Г.В. Семенова Доктор технических наук, профессор Санкт-Петербургского государственного технологического института (технического...»

«В.Ю. ДАВЫДОВ, В.Б.АВДИЕНКО История спортивного плавания Сталинград - Волгоград Монография 2 Волгоград 2010 ББК. 75.717.5 Д. 138 Рецензенты: доктор педагогических наук, профессор А.А.Сучилин; доктор педагогических наук, профессор А.А.Кудинов Д. 138. Давыдов В.Ю., Авдиенко В.Б. История спортивного плавания Сталинград – Волгоград: Монография // В.Ю. Давыдов, А.Б.Авдиенко. – Волгоград: ФГОУ ВГАФК, 2011. – 212 с. Монография адресована в первую очередь профессорскопреподавательскому составу,...»

«Министерство образования и науки Республики Казахстан Институт зоологии П.А. Есенбекова ПОЛУЖЕСТКОКРЫЛЫЕ (HETEROPTERA) КАЗАХСТАНА Алматы – 2013 УДК 592/595/07/ ББК 28.6Я7 Е 79 Е 79 Есенбекова Перизат Абдыкаировна Полужесткокрылые (Heteroptera) Казахстана. Есенбекова П.А. – Алматы: Нур-Принт, 2013. – 349 с. ISBN 978-601-80265-5-3 Монография посвящена описанию таксономического состава, распространения, экологических и биологических особенностей полужесткокрылых Казахстана. Является справочным...»

«Южный научный центр РАН Институт социально-экономических и гуманитарных исследований В.В. Кондратьева, М.Ч. Ларионова ХУДОЖЕСТВЕННОЕ ПРОСТРАНСТВО В ПЬЕСАХ А.П. ЧЕХОВА 1890-х – 1900-х гг.: МИФОПОЭТИЧЕСКИЕ МОДЕЛИ Ростов-на-Дону 2012 УДК 821.161.1.0 ББК 83.3(2Рос–Рус)1 Программа фундаментальных исследований Президиума РАН Традиции и инновации в истории и культуре Проект Художественная литература как способ сохранения, трансляции и трансформации традиционной культуры Кондратьева В.В., Ларионова...»

«Государственное образовательное учреждение высшего профессионального образования ДАЛЬНЕВОСТОЧНАЯ АКАДЕМИЯ ГОСУДАРСТВЕННОЙ СЛУЖБЫ МИНИСТЕРСТВО ОБРАЗОВАНИЯ ХАБАРОВСКОГО КРАЯ ИНСТИТУТ ПРИЕМНОЙ СЕМЬИ: ОПЫТ КОМПЛЕКСНОГО ИССЛЕДОВАНИЯ Хабаровск-2008 4 ББК 60.542.4 И-712 Авторский коллектив: Байков Н.М., доктор социологических наук, профессор; Каширина Л.В., доктор психологических наук, профессор; Березутский Ю.В., кандидат социологических наук, доцент; Иванцева И.А., кандидат социологических наук;...»

«Е.А. Урецкий Ресурсосберегающие технологии в водном хозяйстве промышленных предприятий 1 г. Брест ББК 38.761.2 В 62 УДК.628.3(075.5). Р е ц е н з е н т ы:. Директор ЦИИКИВР д.т.н. М.Ю. Калинин., Директор РУП Брестский центр научно-технической информации и инноваций Государственного комитета по науке и технологиям РБ Мартынюк В.Н Под редакцией Зам. директора по научной работе Полесского аграрно-экологического института НАН Беларуси д.г.н. Волчека А.А Ресурсосберегающие технологии в водном...»

«РОССИЙСКИЙ УНИВЕРСИТЕТ ДРУЖБЫ НАРОДОВ ФИЛОЛОГИЧЕСКИЙ ФАКУЛЬТЕТ КАФЕДРА ПСИХОЛОГИИ И ПЕДАГОГИКИ Гагарин А.В. ЭКОЛОГИЧЕСКАЯ КОМПЕТЕНТНОСТЬ ЛИЧНОСТИ: ПСИХОЛОГО-АКМЕОЛОГИЧЕСКОЕ ИССЛЕДОВАНИЕ Монография Москва, 2011 1 Утверждено ББК 74.58 РИС Ученого совета Г 12 Российского университета дружбы народов Работа выполнена при финансовой поддержке РГНФ (проект № 10-06-0938а) Научный редактор: академик РАО, доктор психологических наук, профессор А.А. Деркач Р е ц е н з е н т ы: член-корр. РАО, доктор...»

«Библиотека адвоката Федеральная палата адвокатов Российской Федерации Центр правовых исследований, адвокатуры и дополнительного профессионального образования Федеральной палаты адвокатов Российской Федерации Л. А. Скабелина ПСИХОЛОГИЧЕСКИЕ АСПЕКТЫ АДВОКАТСКОЙ ДЕЯТЕЛЬНОСТИ Монография Москва 2012 УДК 347.965 ББК 67.75 С42 Автор: Л. А. Скабелина, канд. психолог. наук, доцент кафедры адвокатуры и нотариата МГЮА им. О. Е. Кутафина Рецензенты: Е. В. Семеняко, д-р юрид....»

«с? Ч ^ Q 1 X Эскиз-реконструкция Южного берега древней Таврики i Художник Л. Н. Тимофеев Суровый край каменных обвалов. Отсутствуют характерные кипарисы, завезенные архипелажскими греками лишь в XVIII в. Нт города Алупки. Выше его места, на подъеме к -Петри, господствующему над панорамой, видно таврское городище - святилище на горе Крестовой. Такие же городища-святилища видны вдали, слева, на горе Кошка, и справа, на мысе Ай-Тодор, с огнями, тоже, видно, сигнальными. К берегу правят галеры,...»

«П.П.Гаряев ЛИНГВИСТИКОВолновой геном Теория и практика Институт Квантовой Генетики ББК 28.04 Г21 Гаряев, Петр. Г21 Лингвистико-волновой геном: теория и практика П.П.Гаряев; Институт квантовой генетики. — Киев, 2009 — 218 с. : ил. — Библиогр. ББК 28.04 Г21 © П. П. Гаряев, 2009 ISBN © В. Мерки, иллюстрация Отзывы на монографию П.П. Гаряева Лингвистико-волновой геном. Теория и практика Знаю П.П.Гаряева со студенческих времен, когда мы вместе учились на биофаке МГУ — он на кафедре молекулярной...»

«Министерство образования Республики Беларусь УЧРЕЖДЕНИЕ ОБРАЗОВАНИЯ ГРОДНЕНСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ИМЕНИ ЯНКИ КУПАЛЫ В.Н. Черепица ИСТОРИЯ И ПОВСЕДНЕВНОСТЬ В ЖИЗНИ АГЕНТА ПЯТИ РАЗВЕДОК ЭДУАРДА РОЗЕНБАУМА Монография Гродно 2005 УДК 355.124.6 ББК 68.54 Ч46 Рецензенты: кандидат исторических наук, доцент А.Г.Устюгова; кандидат исторических наук, доцент Э.С.Ярмусик. Рекомендовано советом исторического факультета ГрГУ им. Я.Купалы Черепица, В.Н. История и повседневность в жизни агента пяти...»

«В.Ю. Кудрявцев, Б.И. Герасимов ЭКОНОМИЧЕСКИЙ АНАЛИЗ ТОПЛИВНО-ЭНЕРГЕТИЧЕСКОГО КОМПЛЕКСА (НА ПРИМЕРЕ ТАМБОВСКОЙ ОБЛАСТИ) Научное издание КУДРЯВЦЕВ Вадим Юрьевич, ГЕРАСИМОВ Борис Иванович ЭКОНОМИЧЕСКИЙ АНАЛИЗ ТОПЛИВНО-ЭНЕРГЕТИЧЕСКОГО КОМПЛЕКСА (НА ПРИМЕРЕ ТАМБОВСКОЙ ОБЛАСТИ) Монография Редактор З.Г. Ч ер нов а Компьютерное макетирование З.Г. Черново й Подписано в печать 07.07.2005. Формат 60 84 / 16. Бумага офсетная. Печать офсетная. Гарнитура Тimes New Roman. Объем: 5,22 усл. печ. л.; 5,2...»

«РОССИЙСКАЯ АКАДЕМИЯ НАУК ИНСТИТУТ СОЦИАЛЬНОЭКОНОМИЧЕСКОГО РАЗВИТИЯ ТЕРРИТОРИЙ РАН А.А. Шабунова, К.А. Гулин, М.А. Ласточкина, Т.С. Соловьева МОДЕРНИЗАЦИЯ ЭКОНОМИКИ РЕГИОНА: СОЦИОКУЛЬТУРНЫЕ АСПЕКТЫ Вологда 2012 УДК 316.4(470.12) ББК 60.524(2Рос–4Вол) Публикуется по решению М74 Ученого совета ИСЭРТ РАН Работа выполнена при поддержке гранта Российского гуманитарного научного фонда №11-32-03001а Социально-гуманитарный потенциал модернизации России Модернизация экономики региона: социокультурные...»

«Серия Historia Militaris исследования по военному делу Древности и Средневековья Р е д а к ц и о н н ы й с о в е т: Ю. А. Виноградов (Санкт-Петербург, Россия); В. А. Горончаровский (Санкт-Петербург, Россия); Н. Ди Космо (Принстон, США); Б. В. Ерохин (Санкт-Петербург, Россия); А. Н. Кирпичников (Санкт-Петербург, Россия); Б. А. Литвинский (Москва, Россия); А. В. Махлаюк (Нижний Новгород, Россия); М. Мельчарек (Торунь, Польша); В. П. Никоноров (Санкт-Петербург, Россия); В. Свентославский (Гданьск,...»

«УДК 882-1 ББК 84(2Рос-Рус)5 в 93 Редакционная коллегия : Н. В. Высоцкий,С. В. /'Кильцов,А. В. Максимов,В. Б. Назаров, Е. А. Трофимов Составление и комментарии П. Е. Фокина Подготовка текста. научное консультирование и текстологические комментарии С. В. /'Кильцова При составлении комментариев учтены воспоминания современников В. С. Высоцкого и наблюдения исследователей его творчества, зафиксированные в монографиях и научных публикациях, в частности в книгах: /'Кивая 1Кизнь. Штрихи к...»

«В.Н. Иванов, Л.С. Трофимова МОДЕЛИРОВАНИЕ ФОРМИРОВАНИЯ И РАЗВИТИЯ ПАРКОВ МАШИН ДОРОЖНЫХ ОРГАНИЗАЦИЙ Омск 2012 Министерство образования и науки РФ Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования Сибирская государственная автомобильно-дорожная академия (СибАДИ) В.Н. Иванов, Л.С. Трофимова МОДЕЛИРОВАНИЕ ФОРМИРОВАНИЯ И РАЗВИТИЯ ПАРКОВ МАШИН ДОРОЖНЫХ ОРГАНИЗАЦИЙ Монография Омск СибАДИ УДК 625.76. ББК 39.311.-06- И Рецензенты: д-р техн. наук,...»

«С.В.Бухаров, Н.А. Мукменева, Г.Н. Нугуманова ФЕНОЛЬНЫЕ СТАБИЛИЗАТОРЫ НА ОСНОВЕ 3,5-ДИ-ТРЕТ-БУТИЛ-4-ГИДРОКСИБЕНЗИЛАЦЕТАТА 2006 Федеральное агенство по образованию Государственное образовательное учреждение высшего профессионального образования Казанский государственный технологический университет С.В.Бухаров, Н.А. Мукменева, Г.Н. Нугуманова Фенольные стабилизаторы на основе 3,5-ди-трет-бутил-4-гидроксибензилацетата Монография Казань КГТУ 2006 УДК 678.048 Бухаров, С.В. Фенольные стабилизаторы на...»

«ЮРИДИЧЕСКИЙ ФАКУЛЬТЕТ КАЗАНСКОГО УНИВЕРСИТЕТА: два века образования и науки УДК 34 ББК 67Г Ю70 Печатается по решению Юбилейной комиссии по издательской деятельности Казанского университета Научный редактор доктор юрид. наук, профессор И.А.Тарханов Редакционная коллегия: профессор Р.М.Валеев, профессор Ф.Р.Сундуров, профессор М.В.Талан, фотоснимки И.Ф.Сафина Ю70 Юридический факультет Казанского университета: Два века образования и науки. – Казань: Изд-во Казанск. ун-та, 2004. – 180 с. ISBN...»

«Китай: угрозы, риски, вызовы развитию Под редакцией Василия Михеева МОСКОВСКИЙ ЦЕНТР КАРНЕГИ Москва 2005 УДК 327(510) ББК 66.2(5Кит) К45 Рецензенты: доктор экономических наук, профессор С. С. Суслина, доктор исторических наук С. Г. Лузянин China: Threats, Risks and Challenges to Development Электронная версия: http://www.carnegie.ru/ru/pubs/books Издание осуществляется на средства некоммерческой неправительст венной исследовательской организации — Фонда Карнеги за Междуна родный Мир при...»







 
© 2013 www.diss.seluk.ru - «Бесплатная электронная библиотека - Авторефераты, Диссертации, Монографии, Методички, учебные программы»

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