Содержание к диссертации
Введение
1, Математическое моделирование и дифференциально-алгебраические уравнения 26
1.1. Модели, содержащие процессы с существенно различными характерными временами 26
1.1.1. Математические модели в биологии и медицине 27
1.1.2. Математические модели в электротехнике 31
1.1.3. Корректность представления математических моделей системами дифференциально-алгебраических уравнений 33
1.2. Системы дифференциально-алгебраических уравнений 36
1.2.1. Системы линейных дифференциально-алгебраических уравнений с постоянными коэффициентами 37
1.2.2. Системы линейных дифференциально-алгебраических уравнений с переменными коэффициентами 41
1.2.3. Дифференциально-алгебраические уравнения общего вида . 52
1.2.4. Задача Коши с алгебраической связью на фазовые неременные 54
2. Одношаговые численные методы 60
2.1. Асимптотические оценки погрешности итерационных методов . 60
2.1.1. Метод простых итераций 60
2.1.2. Модифицированный метод Ньютона 65
2.1.3. .Лемма о выпуклости 67
2.2. Методы Рунге-Кутты 68
2.2.1. Построение итеративных методов Рунге-Кутты 68
2.2.2. Сходимость итеративных методов Рунге-Кутты 71
2.2.3. Дифференциальные уравнения 89
2.3. Одношаговые методы общего вида 90
2.3.1. Построение итеративных одношаговых методов 90
2.3.2. Сходимость итеративных одношаговых методов 92
2.4. Численные примеры 96
2.4.1. Простые итерации 96
2.4.2. Ньютоновские итерации 101
3. Практическая реализация и эффективность итеративных методов Рунге-Кутты 104
3.1. Оптимизация простых итераций 105
3.1.1. Стандартные методы 105
3.1.2. Диагонально оптимальные методы 106
3.1.3. Эффективность диагонально оптимальных методов 111
3.1.4. Псевдодиагоналыю оптимальные методы 114
3.1.5. Практическая реализация 117
3.2. Итеративные методы Рунге-Кутты с нетривиальным предиктором . 122
3.2.1. Построение методов с нетривиальным предиктором 122
3.2.2. Численные примеры 124
3.2.3. Сходимость методов с нетривиальным предиктором 127
3.2.4. Практическая реализация 135
3.3. Оптимизация ньютоновских итерации 137
3.3.1. Способы реализации ньютоновских итераций 138
3.3.2. Оптимизация по времени 143
3.3.3. Оптимизация по памяти 148
3.3.4. Численные примеры 148
4. Автоматический контроль точности для одношаговых методов 151
4.1. Контроль локальной ошибки 151
4.1.1. Одношаговые методы с переменным шагом 152
4.1.2. Управление размером шага интегрирования 153
4.1.3. Численные примеры 156
4.2. Разложение глобальной ошибки 159
4.2.1. Дифференциально-алгебраические уравнения 159
4.2.2. Дифференциальные уравнения 164
4.3. Контроль глобальной ошибки 164
4.3.1. Неявное локально-глобальное управление шагом интегрирования 165
4.3.2. Неявные методы 168
4.3.3. Жесткие задачи t .,v 175
4.3.4. Полуявное локально-глобальное управление шагом интегрирования 181
4.3.5. Дифференциальные уравнения 188
5. Одношаговые экстраполяционные методы 190
5.1. Неявная экстраполяция 190
5.1.1. Экстраполяционные методы 191
5.1.2. Неявные экстраполяционные методы 194
5.2. Квадратичная экстраполяция 200
5.2.1. Симметричные одношаговые методы 201
5.2.2. Симметричные методы Рунге-Кутты 207
5.2.3. Неявная квадратичная экстраполяция 212
5.3. Минимально неявные методы 217
5.4. Дифференциально-алгебраические уравнения 226
5.4.1. Неявная экстраполяция 226
5.4.2. Симметричные одношаговые методы 228
5.4.3. Квадратичная экстраполяция 236
5.4.4. Практическая реализация 238
6. Многошаговые численные методы 242
6.1. Линейные многошаговые методы 242
6.1.1. Построение итеративных многошаговых методов 242
6.1.2. Сходимость итеративных многошаговых методов 244
6.1.3. Нетривиальный предиктор 253
6.2. Численные примеры 256
6.2.1. Формулы дифференцирования назад 256
6.2.2. Методы Адамса 260
7. Автоматический контроль точности для многошаговых методов 263
7.1. Управление локальной ошибкой 263
7.1.1. Многошаговые методы с переменным шагом 264
7.1.2. Вычисление главного члена локальной ошибки 267
7.2. Контроль глобальной ошибки 272
7.2.1. Локально-глобальное управление размером шага интегрирования272
7.2.2. Неявные многошаговые методы 281
7.2.3. Численные примеры 285
7.3. Многошаговые экстраполяционные методы 292
7.3.1. Многошаговая экстраполяция 293
7.3.2. Модифицированное локально-глобальное управление шагом интегрирования 297
Литература 301
- Корректность представления математических моделей системами дифференциально-алгебраических уравнений
- Сходимость итеративных методов Рунге-Кутты
- Сходимость методов с нетривиальным предиктором
- Неявное локально-глобальное управление шагом интегрирования
Корректность представления математических моделей системами дифференциально-алгебраических уравнений
В параграфе 1.2 дадим обзор литературы, посвященной как самим дифференциально-алгебраическим уравнениям, так и способам их решения. И прежде всего отметим, что данная тематика занимает особое место в науке. Во-первых, системы дифференциально-алгебраических уравнений, а значит, и численные методы дня них имеют большое прикладное значение, так как позволяют исследовать многие реальные процессы в различных естественно-научных областях знания (см. параграф 1.1). Во-вторых, несмотря на то, что такие задачи имеют много общего с системами обыкновенных дифференциальных уравнений, они обладают рядом существенных отличий, изучение которых приводит к появлению новых фундаментальных понятий и теоретических результатов.
Итак, как мы упоминали во введении, первые работы по этой тематике появились сравнительно недавно. Их можно датировать концом шести десятых-началом семидесятых годов (см., например, [9], [78], [88], [114]). Хотя отдельные результаты, относящиеся к линейным системам с постоянными коэффициентами, были получены несколько ранее.
В последние двадцать лет развитие теории и методов решения систем дифференциально-алгебраических уравнений происходило особенно интенсивно. Уже можно назвать множество монографий, посвященных этому вопросу: [10]—[12], [83], [91], [96], [97], [137], [139]. Однако, несмотря на несомненные достижения, в настоящее время еще не существует устоявшейся терминологии, связанной с дифференциально-алгебраическими уравнениями. Так, одни авторы называют системы такого рода дифференциально-алгебраическими (Gear C.W., Petzold L.R., Hairer Е., Wanner G., Brenan К.E., Engquist B.E., Reinboldt W.C., Кузнецов Е.Б., Шалашилин В.И), другие — сингулярными системами (Бояринцев Ю.Е., Данилов В.А., Логинов А.А., Чистяков В.Ф., Campbell S.L.), третьи обозначают их как нелинейные дифференциальные уравнения с алгебраическими связями (Lotstedt Р.) и, наконец, встречаются понятия "неявные дифференциальные уравнения77 (Campbell S.L.) и "дифференциальные уравнения на многообразиях" (Reinboldt W.C.). Поэтому введем некоторые определения, на которые будем опираться в дальнейшем.
В. случае обыкновенных дифференциальных уравнение если 5 = П Sj и z"(9j) ограничены, a Az0 — О, то из (1.32) следует сходимость приближенного решения к точному. Для дифференциально-алгебраических уравнений эти условия могут не выполняться. Во-первых, так как решение такой задачи связано с дифференцированием функции f(t) (см. раздел 1.2.1), то z (i) может иметь разрывы, следствием чего может явиться неограниченность z"(t). Вснвторых, Дт-о может не стремиться к нулю. В-третьих, матрица 5& может быть неограниченной, т.е. некоторые элементы этой матрицы могут неограниченно возростать по модулю.
Теорема 1.2.4. Пусть с421 d 1 и \\Qj\\, \\Qj\\, \\Qjl\\ ограничены для всех j. Тогда неявный метод Эйлера сходится для случая п = 2 при условии, что начальная огиибка имеет порядок малости о(т) и z" ограничена.
Теорема 1.2.4 показывает, что для задачи (1.25) при п = 2 еще удается получить сходимость. Но если рассмотреть задачи более высокой размерности (а следовательно, и более высокого индекса), то в матрице $& появятся ненулевые члены порядка 0(1 /г) и выше, что приведет к неограниченному росту $Р и, как следствие, к расходимости численного метода.
Таким образом, мы приходим tf выводу, что для интегрирования систем линейных дифференциально-алгебраических уравнений с переменными коэффициентами требуется их предварительное преобразование, приводящее к задачам, решаемым численными методами. Далее обсудим некоторые из таких преобразований.
Сходимость итеративных методов Рунге-Кутты
Во второй главе, состоящей из четырех параграфов, развита теория комбинированных одношаговых методов для решения систем дифференциально-алгебраических уравнений индекса 1. Такие методы строятся на основе одношаговых методов для дискретизации дифференциальных уравнений и итерационных процессов для решения нелинейных алгебраических. Полученные таким образом алгоритмы исследованы на сходимость. Для них доказаны оценки погрешности и выведены достаточные условия, гарантирующие сходимость максимального порядка. Теоретические результаты подтверждены численными примерами.
. Асимптотические оценки погрешности итерационных методов
В параграфе 2.1 рассматриваются наиболее известные и часто используемые на практике представители семейства итерационных методов для решения нелинейных алгебраических уравнений, а именно, — метод простых итераций и модифицированный метод Ньютона. Для этих методов получены асимптотические оценки погрешности, играющие важнейлгую роль в анализе сходимости комбинированных численных методов. Кроме того, здесь приводятся и другие дополнительные результаты необходимые для дальнейшего изложения. . В этом случае будет справедлива следующая теорема о методе простых итераций [59, с. 401]. Теорема 2.1.1. Пусть отображение G : D С Rm —» Rm F-дифференцируемо на выпуклом множестве DQ С D и производная этого отображения непрерывна по Липшицу, т. е. d&V) - dG{x")\\ 7s - х"\\ V х\х" Є Do. и допустим, что замкнутый шар1 S(x,p ) С D0. Тогда итерации (2.2) лежат в S(xiP ) и сходятся к единственному решению х системы (2.1) в S(x,p ) f] Do. При этом справедлива оценка ошибки s - N p"-PN, ЛГ-1,2,..., (2.3) г е PTV = І7Рлг_і + PN-I +17, Ро = 0. Итак, теорема 2.1.1 не только дает достаточные условия сходимости метода простых итераций, но и позволяет оценить погрешность этого метода через сходящуюся последовательность {рк}. Однако оценка (2.3) весьма неудобна для практического использования. Она не дает реального представления о скорости сходимости из-за правой части этой формулы, опредеіенной рекуррентно. Поэтому, учитывая неявный вид оценки (2.3), докажем теорему об асимптотическом разложении ошибки метода простых итераций при а —У 0. В результате мы не только доказали, что оценка (2.21) точнее оценки (2.22), но и выяснили, что последняя оценка остается справедливой и при более слабых предположениях относительно гладкости задачи (2.19). STRONG Теорема 2.1.5. Пусть выполнены условия теоремы 2.1.4. Тогда для. ошибки модифицированного метода Ньютона справедлива оценка (2.22). Следствием теоремы 2.L5 является тот факт, что при a — 0 погрешность N-ой рітерации модифицированного метода Ньютона ведет себя, по крайней мере, как 0(ог+). Однако с учетом (2.23) интересно выяснить асимптотическое поведение погрешности модифицированного метода Ньютона для более точной оценки (2.21). Теорема 2.1.6. Пусть выполнены условия теоремы 2.1-4- Тогда для ошибки модифицированного метода Ньютона при достаточно малом а справедлива оценка \\х - xN\\ -a v+1, ЛГ = 1,2, (2.25) где С — некоторая константа, независящая от у. Доказательство теоремы 2.1.(1 практически дословно повторяет доказательство теоремы 2.1.2 с учетом специфики производящей последовательности {рм}. Однако в последнем случае оно значительно проще технически (см. [30]).
Итак, мы показали, что хотя оценка погрешности модифицированного метода Ньютона (2.21) точнее (2.22), обе эти оценки являются асимптотически эквивалентными, так как при достаточно малом а дают сходимость модифицированного метода Ньютона со скоростью геометрической прогрессии.
В заключение этого раздела отметим, что оценка (2.25) хорошо соответствует оценке ошибки для полного метода Ньютона. Из теоремы Канторовича следует, что при одной итерации ошибка полного метода Ньютона ведет себя как 0(а2) при а - 0 (см. [59, с. 404]). С другой стороны, при одной итерации полный и модифицированный методы Ньютона совпадают, поэтому они должны иметь одинаковую погрешность, что следует из (2.25) при N = 1.
Сходимость методов с нетривиальным предиктором
В заключение отметим, что теоремы 3.2.2-3.2.5 сохраняют свою актуальность и в случае произвольных одиопгаговьгх методов, т. е. для методов (2.109)—(2.111) с нетривиальным предиктором. Надо только учесть, что в общем случае отсутствует понятие стадийного порядка и. Поэтому он должен быть исключен из оценок ошибок (3.69), (3.71), (3.73) и (3.75) для методов (2.109)—(2.111) с нетривиальным предиктором.
Обсудим теперь некоторые аспекты практической реализации численных методов с нетривиальным предиктором, а также дополнительные возможности для снижения числа итераций в точке сетки.
Во-первых, как показывают теоремы 3.2.2-3.2.5 для методов с нетривиальным предиктором существенное значение имеет точность задания стартовых величин 5 . к = ОЛ , / — 1, т. е. если мы хотим достичь сходимость порядка базовой РК формулы s. то такую же точность необходимо обеспечить при определении стартовых значений. Это подраз\ мевает выполнение первых УСЛОВИЙ В (3.70). (3.72) и (3.74). Однако с практической точки зрения эти условия не существенны, так как если начальный вектор (2.28в) известен точно (или, по-крайней мере, с точностью 0(rs)), то, применяя для вычисления гь к = 1,2,...,/ - 1, РКН или РКМН-метод с тривиальным предиктором, при некотором числе итераций можно добиться их выполнения (см. раздел 2.2.2). Затем число итераций в оставшихся точках сетки можно уменьшить в соответствии с (3.70), (3.72) или (3.74). Более того, нетрудно разработать процедуру разгона для методов с нетривиальным предиктором, при которой на каждом шаге порядок предиктора увеличивается на 1, а число итераций сокращается.
Во-вторых, при правильной организации вычислений число итераций для решения дифференциальных уравнений РКПИ, РКН и РКМН-методами может быть и меньше, чем предсказано формулами (3.76)-(3.78). Рассмотрим дискретную задачу (2.29). Для удобства изучения вопроса в теоретическом плане мы до сих пор считали, что тот или иной итерационный процесс применяется ко всем уравнениям этой задачи одновременно. Однако на практике поступают по-другому. В силу того, что уравнение (2.29в) разрешимо явно, итерации применяются только при вычислении стадийных величин (2.29а)5. Это позволяет сократить затраты машинного времени за счет уменьшения размерности решаемой задачи. Кроме того, следует отметить, что при подстановке найденных стадийных величии x t в (2.29в) ошибка итерационного метода будет домножаться на г. Последнее обстоятельство означает, что вклад итерационного процесса в ошибку комбинированного метода будет в данном случае на порядок меньше. Поэтому достаточное число итераций можно сократить в соответствии со следующими оценками:
Таким образом, кроме уменьшения порядка решаемой в итерационном процессе алгебраической системы мы дополнительно получили сокращение числа итераций, что еще больше увеличивает эффективность численного интегрирования дифференциальных уравнений неявными РК-методами. Очевидно, что аналогичное улучшение мы получаем и при правильном применении РКПИ, РКН и РКМН-методов с тривиальным предиктором к дифференциальным уравнениям. В этом случае достаточно в оценках для числа итераций (3.79)—(3.81) положить = 1 (сравнить с (2.100)-(2.102)).
К сожалению, при решении дифференциально-алгебраических уравнений, подобное сокращение числа итераций уже не наблюдается. Этому препятствует алгебраическое соотношение (2.29г). Однако даже здесь можно получить некоторый выигрыш во времени.
Напомним, что мы рассматриваем в ланнын момент только дифференциальные уравнения. поэтому алгебраические соотношения (2.296. г) и ц/-компоненты отсутствуют.
Нетрудно видеть, что решение задачи (2,29) можно разбить на три этапа. Сначала применяется один из вышеназванных итерационных методов для вычисления стадийных величин z . В этом случае мы решаем подзадачу (2.29а, 6) размерности 1(т + п). Затем х-компоненты приближенного решения задачи (2.28) находятся явно с помощью (2.29в). И наконец, опять применяя итерационный метод для решения (2.29г), вычисляем /-компоненты численного решения. Здесь мы решаем подзадачу размерности п. Такая реализация процесса численного интегрирования задачи (2.28) позволяет существенным образом сократить затраты машинного времени.
Так вот, при решении меньшей системы (последний этап) число итераций действительно сокращать нельзя. И исходя из теорем 3.2.2, 3.2.3, его нужно выбирать в соответствии с (2.70) или (2.72)6. Однако при решении большей системы (первый этап) итераций может быть меньше, так как на следующем этапе, т. е. при вычислении х-компонентов приближенного решения, ошибка итерационного метода будет уменьшена на порядок за счет умножения на т. Кроме того, здесь нужно иметь в виду, что точность вычисления стадийных величин не влияет на точность вычисления у-компонентов приближенного решения (см. (2.29г)). Поэтому при решении подзадачи (2.29а, б) число итераций следует выбирать, исходя из оценок: Таким образом, машинное время при решении дифференциально-алгебраических систем может быть дополнительно съэкономлено за счет уменьшения числа итераций при вычислении стадийных величин итеративных РК-мстодов (для методов с тривиальным предиктором в (3.82) и (3.83) достаточно положить = 1),
В параграфе 3.3 изучаются особенности применения ньютоновских итераций для решения систем дифференциально-алгебраических уравнений индекса 1 большой размерности неявными РК-методами. Численное интегрирование в этом случае на практике сводится к многократному решению систем линейных алгебраических уравнений с большими матрицами коэффициентов, что приводит к огромным затратам машинного времени и требует значительного объема оперативной памяти ЭВМ. Поэтому с учетом структуры таких матриц разработана специальная модификация гауссовского исключения неизвестных, позволяющая существенным образом сократить затраты машинного времени и памяти. Кроме того, здесь рассматривается один подход к реализации ньютоновских итерашш. основанный на негауссовском обращении матриц, что может быть весьма благоприятно для использования параллельных или векторных ЭВМ. Выше мы уже отмечали, что оценка (2.74) имеет в основном только теоретическое значение.
Неявное локально-глобальное управление шагом интегрирования
К сожалению, результаты последнего эксперимента не позволяют говорить о каком-либо преимуществе локально- глобально контроля точности численного решения над стандартным локальным для тестовой задачи (4.37) (см. табл. 4.18). Как мы видим, усложнение способа управления размером шага интегрирования не привело к ожидаемому увеличению точности численного решения.
Для объяснения этого феномена вспомним, что жесткие задачи характеризуются тем, что их якобиан обладает большими по модулю собственными числами no-крайней мере на некоторых интервалах отрезка интегрирования. Поэтому на таких интервалах скорость изменения глобальной ошибки численного метода значительно возрастает (см. (4.10)). Тогда шаг т , который обеспечивал необходимую малость локальный ошибки, может быть слишком большим для получения достоверной оценки глобальной ошибки с помощью (4.33) или (4.38). Более того, погрешность этих формул вообще может превышать глобальную ошибку численного метода, т. е. найденное значение глобальной ошибки не будет содержать НЕ одной значащей цифры.
Очевидно, что управление шагом интегрирования, которое базируется на такой информации не может давать приемлемые по точности результаты. Таким образом, становиться понятно, что при решении жестких задач с заданной точностью важное значение приобретает контроль погрешности численного решения системы дифференциально-алгебраических уравнений (4.10).
Для реализации этой идеи на практике воспользуемся формулами (4.33) и (4.38). Мы знаем, что они дают главный член глобальной ошибки метода (2.108) с разной точностью. Поэтому если взять их разность, то получим оценку компонентов локальной ошибки метода (4.33) (Д «( +і) и &Ф {Ік+і)) с точностью 0(т%). Окончательно, выбирая границу для этой погрешности елг, приходим к устойчивому неявному локально-глобальному контролю точности численного решения для од-ношаговых методов (см. рис. 4.5).
Проверим теперь работоспособность нового способа выбора шага интегрирования на тестовой задаче (4.37). Для этого решим ограниченную задачу трех тел КБН-методом с предиктором второго порядка и устойчивым неявным локально-глобальным контролем точности при тех же самых границах для локальной и глобальной ошибок, что и в таблицах 4.17, 4.18. Положим елг = ег/100.
Как мы видим, устойчивый вариант локально-глобального контроля точности численного решения внушает определенный оптимизм (см. табл. 4.19), так как в отличие от двух предыдущих экспериментов КБН-метод с новым управлением размером шага численного интегрирования везде достигает заданную точность вычислений, если число ньютоновских итераций в точке сетки не противоречит условию (4.6). Действительно, если мы подставим з = 6и = 3в (4.6) (так как стадийный и классическрш порядки КБ-метода равны соответственно трем и шести, а порядок ошибки интерполирования равен трем), то получим, что двух итераций Ньютона в точке сетки будет достаточно (см. вторую строку в табл. 4.19).
В заключение этого раздела продемонстрируем еще раз прекрасные потребительские свойства устойчивого локально-глобального контроля точности численного решения, но теперь для методов низкого порядка. В качестве базового РК-метода возьмем правило трапеций и построим ПТН-метод с предиктором второго порядка.
Протестируем этот метод па ограниченной задаче трех тел.
Как следует из данных таблицы 4.20, в этом случае нам также удается достичь любую наперед заданную точность для численного решения. Однако, чтобы убедиться в преимуществе нового способа управления шагом интегрирования более наглядно, приведем еще данные для ПТН-метода, но с локальным контролем точности, при тех же самым границах4для ошибки (см. табл. 4.21). Как и следовало ожидать, локальный контроль не только не обеспечивает требуемую точность вычислений, но и вообще дает весьма большую погрешность при одной или двух итерациях Ньютона в точке сетки. В отличие от него устойчивый локально-глобальный контроль достигает требуемую точность уже при одной итерации, что полностью соответствует оценке достаточного числа итераций Ньютона (4.6).
В предыдущих разделах мы предложили и исследовали некоторые способы контроля глобальной ошибки одношаговых методов для нежестких и жестких дифференциально-алгебраических уравнений индекса 1, Все эти способы базировались на применении неявных одношаговых методов к системе (4.10) для получения достоверной оценки главного члена глобальной ошибки и последующем управлении шагом интегрирования с целью его контроля. В результате было показано, что одношаговые методы с локально-глобальным управлением размером шага интегрирования действительно позволяют достигать любую разумную наперед заданную точность для численного решения в автоматическом режиме. В настоящем разделе сосредоточим свое внимание на эффективности новых алгоритмов.
Как и выше, для исследования эффективности различных способов управления шагом интегрирования рассмотрим средние шаги, которые дают соответственно локальный и устойчивый локально-глобальный контроль точности для КБН-метода, примененного к задаче (4.37). В таблице 4.22 приведены экспериментальные данные для локального контроля точности численного решения, а в таблице 4.23 — для локапьно-глобального.
Основываясь на этой информации, мы приходим к следующим выводам. Во-первых, в обоих случаях максимальный размер для среднего шага достигается только при двух итерациях в точке сетки, что подтверждает справедливость оценки (4.6). Во-вторых, средние шаги при локальном и устойчивом локально-глобальном контроле точности численного решения отличаются не более чем в два раза при одинаковых границах для локальной и глобальной ошибок. Причем, если требуется не слишком высокая точность вычислений, то разница становиться еще меньше (см. первые два столбца в табл. 4.22, 4/23). Однако если теперь мы обратимся к непосредственным затратам машинного времени, то окажется, что локально-глобальное управление шагом интегрирования требует в три раза больших затрат. Это легко объясняется высокой стоимостью вычисления локальной ошибки КБН-метода в точке fjfe+i Более того, мы должны пересчитывать г +2 и V s+i(2fc+i) в процессе численного интегрирования всякий раз, когда этого требует устойчивый локально-глобальный контроль точности. Таким образом, задачей настоящего раздела является упрощение устойчивого локально-глобальной контроля точности численного решения с целью повышения его эффективности в смысле сокращения затрат машинного времени.
