Содержание к диссертации
Введение
1 Обзор литературы 9
1.1 Одношаговые методы 10
1.2 Экстраполяционные методы 16
1.3 Автоматический контроль точности вычислений 24
2 Неявные экстраполяционные методы и локальный контроль точности 36
2.1 Автоматическое управление размером шага и порядком в явных одно-шаговых экстраполяционных методах 36
2.2 Автоматическое управление размером шага и порядком в неявных од-ношаговых экстраполяционных методах 44
2.3 Вычислительные эксперименты 47
3 Неявные экстраполяционные методы и глобальный контроль точности 65
3.1 Локально-глобальное управление размером шага и порядком экстраполяционных методов 65
3.2 Вычислительные эксперименты 70
4 Коллокационные методы со старшими производными и глобальный контроль точности 81
4.1. Присоединенные и симметричные РК-методы со старшими производными 81
4.2 Эффективная реализация Е-методов и упрощенные итерации Ньютона 87
4.3 Локально-глобальное управление размером шага и порядком в одногаа-говых коллокационных методах со старшими производными 96
4.4 Вычислительные эксперименты 97
Заключение 105
Литература 108
- Автоматический контроль точности вычислений
- Автоматическое управление размером шага и порядком в неявных од-ношаговых экстраполяционных методах
- Эффективная реализация Е-методов и упрощенные итерации Ньютона
- Локально-глобальное управление размером шага и порядком в одногаа-говых коллокационных методах со старшими производными
Введение к работе
Актуальность темы. Системы обыкновенных дифференциальных уравнений (ОДУ) представляют собой неотъемлемую часть математического моделирования динамических объектов. Однако размер и сложность современных моделей делает аппарат аналитических вычислений практически неприемлемым в данной ситуации. Поэтому проблеме численного интегрирования задач вида
x\t) = g(t,x(t)), *Є[0,Т], (la)
ж(0) = х, (16)
где правая часть д : Rm+ —> Rm является достаточно гладкой функцией, отводится существенное место в вычислительной математике. Такие задачи возникают как непосредственно при математическом моделировании в области биологии, медицины, механики, электротехники, химии и т. д., так и при решении более сложных систем уравнений (например, при дискретизации задач математической физики методом прямых).
В настоящее время к наиболее перспективным алгоритмам численного решения задачи Коши для ОДУ можно отнести экстраполяци-онные методы. Общая теория экстраполяционных методов, использующих неявные одношаговые схемы, была впервые предложена Куликовым Г.Ю.1 Наряду с хорошо изученными стандартными одношаговыми методами Рунге-Кутты в качестве опорных для построения экстраполяционных схем, можно предложить сравнительно новые одношаговые методы Рунге-Кутты, использующие старшие производные (см., например, Аульченко СМ., Латыпов А.Ф., Никуличев Ю.В.2, Куликов Г.Ю.,
хКиыкоу G.Yu. On implicit extrapolation methods for ordinary differential equations// Russian Journal of Numerical Analysis and Mathematical Modelling. — 2002. — V. 17, № 1. - P. 41-69.
2Аульченко СМ., Латыпов А.Ф., Никуличев Ю.В. Метод численного интегрирования систем обыкновенных дифференциальных уравнений с использованием интерполяционных полиномов Эрмита// Журнал вычислительной математики и математической физики. - 1998. - Т. 38, № 10. - С. 1665-1670.
Меркулов А.И.3, Fehlberg Е. 4, Kasthmger К.Н.5, N0rsett S.P.6). Однако, исследования в этой области пока достаточно ограничены, а определение условий симметричности таких методов и использование их в качестве базовых для построения экстраполяции не изучалось вовсе.
Далее следует отметить, что важнейшим этапом в эффективной реализации любого численного метода для интегрирования обыкновенных дифференциальных уравнений является разработка механизма автоматического управления размером шага с целью обеспечения требуемой точности вычисления за минимальное (или приемлемое) время. Дополнительным преимуществом экстраполяционных методов является возможность автоматического подбора оптимальных размера шага и порядка в каждой точке сетки. Вопрос о надежной и достоверной оценке точности численного решения ОДУ неразрывно связан с проблемой оценки глобальной ошибки, которой посвящено достаточно большое число работ. В то время как обсуждению алгоритмов автоматического контроля глобальной ошибки уделено значительно меньше внимания. Можно указать лишь относительно небольшое количество статей, так или иначе затрагивающих эту тематику: Новиков Е.А.7, Dahlquist G.8, Lindberg В.9, Skeel R.D.10, Stetter H.J.11. А важное преимущество экстраполяционных методов, позволяющее управлять не только размером шага, но и порядком самого
3Куликов Г.Ю., Меркулов А.И. Об одношаговых коллокационных методах со старшими производными для решения обыкновенных дифференциальных уравнений// Журнал вычислительной математики и математической физики. — 2004. — Т. 44, № 10. - С. 1782-1807.
4Fehlberg Е. New high-order Runge-Kutta formulas with step size control for systems of first and second order differential equations// ZAMM. — 1964. — V. 44, Sonderheft — P.T17-T19.
5Kastlunger K.H., Wanner G. Runge-Kutta processes with multiple nodes// Computing. - 1972. - V. 9. - P. 9-24.
6Norsett S.P. One-step methods of Hermite type for numerical integration of stiff systems// BIT. - 1974. - V. 14. - P. 63-77.
7Новиков Е.А. Оценка глобальной ошибки А-устойчивых методов решения жестких систем// Доклады академии наук. — 1995. — Т. 343, № 4. — С. 452-455.
8Dahlquist G. On the control of the global error in stiff initial value problems// Lecture Notes in Mathematics. — Berlin-New York: Springer, 1982. — V. 912. — P. 38-49.
9Lindberg B. Characterization of optimal stepsize sequences for methods for stiff differential equations// SIAM Journal on Numerical Analysis. — 1977. — V. 14, № 5. — P. 859-887.
10Skeel R.D. Thirteen ways to estimate global error// Numerische Mathematik. — 1986. - V. 48. - P. 1-20.
"Stetter H.J. Local estimation of the global discretization error// SIAM Journal on Numerical Analysis. - 1971. - V. 8. - P. 512-523.
метода, практически не используется.
Таким образом, объектом исследования являются неявные одногла-говые экстраполяционные схемы, а предметом исследования - алгоритмы комбинированного управления шагом и порядком экстраполяционных процессов, позволяющие находить численное решение, удовлетворяющее заданным требованиям к точности, в автоматическом режиме.
Цель и задачи работы. Целью работы является построение и тестирование новых численных методов, использующих старшие производные и экстраполяционную технику, для решения обыкновенных дифференциальных уравнений. А также разработка теории комбинированного управления размером шага численного интегрирования и порядком одноша-говых методов (в том числе и коллокационных методов со старшими производными) с использованием экстраполяции, позволяющей достичь наперед заданную точность вычислений в автоматическом режиме. В соответствии поставленной целью решались следующие задачи:
разработка алгоритмов неявного локально-глобального управления порядком и шагом экстраполяционных методов;
исследование свойств одношаговых методов Рунге-Кутты, использующих старшие производные;
эффективная реализацию разработанных в диссертации алгоритмов с использованием новых упрощенных итераций Ньютона;
разработка комплекса программ для численного интегрирования систем обыкновенных дифференциальных уравнений с заданной точностью в автоматическом режиме.
Методы исследования. В диссертационной работе использованы методы вычислительной математики, математического и функционального анализа, линейной алгебры и программирования.
Научная новизна. Все основные результаты настоящей диссертации являются новыми и состоят в следующем:
Разработаны алгоритмы контроля локальной и глобальной ошибок неявных экстраполяционных методов для систем обыкновенных дифференциальных уравнений.
Построены экстраполяционные схемы, использующие коллокацион-ные методы со старшими производными, для нахождения численного
решения систем обыкновенных дифференциальных уравнений и изучены их фундаментальные свойства.
Разработаны алгоритмы контроля локальной и глобальной ошибок экстраполяционных методов, использующих коллокационные методы со старшими производными, для систем обыкновенных дифференциальных уравнений.
Предложена эффективная схема упрощенных итераций Ньютона, оптимизирующих реализацию разработанных в диссертации алгоритмов.
Основные положения, выносимые на защиту:
Методика построения присоединенных и симметричных коллокаци-онных методов со старшими производными, формулы для вычисления коэффициентов таких методов.
Способы автоматического контроля локальной и глобальной ошибок экстраполяционных методов, использующих неявные одношаговые методы Рунге-Кутты и коллокационные методы со старшими производными, для систем обыкновенных дифференциальных уравнений.
Методика построения упрощенных итераций Ньютона, доказательство сходимости комбинированных схем неявных экстраполяционных методов с использованием нового итерационного процесса, вывод асимптотической оценки погрешности.
Комплекс программ на языке C++, позволяющий интегрировать системы обыкновенных дифференциальных уравнений с заданной точностью в автоматическом режиме.
Достоверность результатов. Достоверность результатов обеспечивается строгостью постановок задач и доказательств теорем, использованием аналитических и численных методов расчета, методов математического моделирования и применением современных методик экспериментального исследования.
Практическая и теоретическая значимость. Полученные в диссертации новые результаты представляют интерес для численного решения систем ОДУ экстраполяционными методами с контролем локальной
и глобальной ошибок. В прикладном аспекте такие методы могут с успехом быть использованы для разработки автоматизированных систем математического моделирования в различных областях науки и техники, поскольку способность обеспечивать заданную пользователем точность вычислений в автоматическом режиме существенно повышает достоверность математического моделирования на практике.
Апробация работы. Результаты диссертации докладывались и обсуждались: на пятой международной научно-технической конференции «Математическое моделирование физических, экономических, технических, социальных систем и процессов» (Ульяновск, 2003), на XXVI конференции молодых ученых механико-математического факультета МГУ им. М.В. Ломоносова (Москва, 2004), на 4th International Conference on Computer Science (Krakow, Poland, 2004).
Личный вклад автора. Постановка задач осуществлялась совместно с научным руководителем, д.ф.-м.н. Г.Ю. Куликовым, теоретические положения, проведение расчетов и анализ полученных результатов выполнены автором самостоятельно.
Публикации. По теме диссертации опубликовано 10 работ, в том числе 3 в рецензируемых научных изданиях, рекомендованных ВАК, их список помещен в конце автореферата.
Структура диссертации. Диссертация состоит из введения, четырех глав, заключения, списка литературы (162 наименования) и приложения. Общий объем диссертации составляет 128 страниц, из них 118 страниц основного текста, 10 страниц приложения. Работа включает 49 таблиц и 20 рисунков.
Автоматический контроль точности вычислений
Итак, в настоящем параграфе исследуется вопрос получения оценок локальной и глобальной ошибок численных методов и эффективного использования таких оценок для контроля точности численного решения.
Идея, лежащая в основе всех методов с автоматическим выбором шага интегрирования, состоит в следующем: на каждом шаге определяется некоторая величина, характеризующая ошибку метода, затем на ее основе вычисляется оптимальная длина для шага интегрирования. При этом достигается определенный компромисс между точностью и временем интегрирования задачи (1.1.1). Традиционным и наиболее простым является использование оценки локальной ошибки (1.1.10), этот подход весьма популярен, хотя и имеет ряд серьезных недостатков. Локальная ошибка соответствует погрешности вносимой при численном решении задачи (1.1.1) на одном шаге. В настоящее время разработано большое число различных способов контроля локальной погрешности, однако наиболее часто численные методы с автоматическим выбором шага интегрирования основаны на вычислении главного члена асимптотического разложения локальной ошибки в ряд Тейлора (см. i s+i(tk) в формуле (1.1.10)) и последугощем выборе такого размера для очередного шага, который является максимальным для заданного предела локальной ошибки.
Решение первой задачи - нахождение главного члена локальной ошибки - в большинстве случаев осуществляется на практике с помощью правила Рунге или численных методов разного порядка (см., например, [4, 5, 33, 41, 43, 44, 51, 100, 101]). Оба способа опираются на следующее предположение об асимптотическом представлении локальной ошибки метода порядка s (см., например, [1, 4, 5, 44]):
При использовании правила Рунге оценивают только первый член в разложении локальной погрешности (1.3.1). Для оценки погрешности данным способом нужно выполнить три шага численного метода. Сначала из точки tk делают шаг размера тк, при этом справедливо соотношение Отсюда, принимая во внимание (1.3.1), нетрудно получить практическую формулу для оценки главного члена локальной ошибки при достаточно малом т . Итак, применение правила Рунге для вычисления главного члена локальной ошибки предполагает вычисление трех шагов численного метода, что может служить серьезным препятствием для его использования в том случае, когда накладные расходы на один шаг численного метода велики. Другой недостаток данного метода состоит в том, что при выводе формулы (1.3.2) существенно использовалась гладкость коэффициента фя(Ь), что не всегда справедливо на практике. Полное теоретическое обоснование этой методики оценки погрешности впервые было предложено Ричардсоном (см. [146, 147]). Другой подход к оценке локальной ошибки основан на использовании численных методов разных порядков. Здесь для получения оценки проводят два интегрирования из точки tk с шагом т методами, имеющими разные порядки. Для определенности будем считать, что для решения задачи применяются методы порядков s и р, причем s р. Предположим, что приближенное решение Xk+i вычислено методом порядка s, a Xk+i — методом порядка р. Допустим также, что правая часть задачи (1.1.1а) р + 1 непрерывно дифференцируема в некоторой окрестности точного решения этой задачи. Тогда, используя асимптотическое разложение локальной ошибки (1.3.1) для каждого из двух численных методов, легко видеть, что локальная ошибка метода порядка s удовлетворяет равенству На практике, чаще всего, ограничиваются оценкой только главного члена локальной ошибки, хотя изложенный выше способ позволяет делать и более точные оценки. Заметим, что способ, опирающийся на два метода разных порядков, требует только два интегрирования исходной задачи (1.1.1) на каждом шаге, что в некоторых случаях дает заметный выигрыш во времени по сравнению с правилом Рунге. Следовательно, были предприняты попытки повысить эффективность именно второго способа оценки локальной погрешности. К примеру, для РК-методов идея заключалась в построении таких РК-формул, которые сами содержали бы кроме численного приближенного значения &k+i некоторое выражение ж +і более низкого порядка. Это могло бы, в частности, удешевить процедуру выбраковки шагов. В результате появились вложенные РК-формулы [44, 100], где для нахождения численного решения методами порядков s и s + 1 применяют одни и те же стадийные величины, что значительно снижает затраты машинного времени, так как основные затраты машинного времени в РК-методах (особенно неявных) приходятся на вычисление стадийных величин. После получения асимптотически корректной оценки локальной погрешности можно переходить к решению второй задачи, а именно: к управлению шагом интегрирования, базирующемуся на полученной оценке. Основная цель всех существующих процедур управления размером шага — определение максимального размера шага интегрирования для заданной границы локальной ошибки, что в большинстве случаев решается аналитически (см., например, [1, 4, 5, 44, 51, 100, 101]). В силу предположения о достаточной гладкости правой части уравнения (1.1.1а), очевидно, что функции 4 i(t), і = s,...,S, будут ограничены на всем отрезке [0,Т]. Таким образом, уменьшая размер шага интегрирования, мы можем сделать локальную ошибку сколь угодно малой (в отсутствии ошибок округления). Это и является причиной значительной популярности численных методов для решения задачи (1.1.1) с переменным шагом интегрирования на основе контроля локальной ошибки. На практике обычно задают некоторый допуск є/ос и требуют, чтобы на каждом шаге интегрирования выполнялось неравенство Тогда выполнение условия (1.3.4) всегда можно обеспечить выбирая подходящий размер шага интегрирования т .. Типичную процедуру управления локальной ошибкой при численном интегрировании задачи (1.1.1) можно описать следующим образом. Пусть в k-oii точке сетки нам известны величины Жь = x(tk) и Tfc. Тогда а) вычисляем величины Xk+i и Дж +1; б) проверяем справедливость условия (1.3.4) для Ах +и в) если оценка локальной ошибки Дхк+і не удовлетворяет (1.3.4), то отбрасываем величину Xk+i, уменьшаем шаг интегрирования Тк и переходим к пункту а); г) если условие (1.3.4) выполнено, то принимаем величину Xk+i в качестве при ближенного решения задачи (1.1.1) в точке tk+i, вычисляем размер нового шага интегрирования Тк+і и переходим к пункту а) в следующей точке сетки wT. В качестве размера шага интегрирования Тк выбирают максимальную величину, при которой выполняется условие (1.3.4). В силу того, что на практике обычно ограничиваются оценкой только главного члена локальной погрешности, условие (1.3.4) принимает вид \шікж+і єІ0С. (1.3.5) Теперь, используя (1.3.5), можно получить простой и дешевый способ вычисления оптимального (в смысле локального контроля) размера шага интегрирования. Итак, пусть на к + 1-ом шаге соотношение (1.3.5) не выполнено. Тогда для вычисления оптимальной длины шага приходим к следующему простому выражению: Если же на fc+1-ом шаге неравенство (1.3.5) имеет место, то, предполагая, что главный член локальной ошибки меняется достаточно гладко, для вычисления нового шага интегрирования также используют формулу (1.3.6), но с заменой т на Тк+\. Следует заметить, что при практической реализации, учитывая погрешность оценки (1.3.6), а также ошибки округления, в формулу вычисления оптимального шага добавляют так называемый гарантийный фактор 7- Он используется для сокращения числа пересчетов численного решения и, как следствие, затрат машинного времени. Обычно это некоторая константа, принимающая значения в промежутке от 0.5 до 0,9. В этом случае формула (1.3.6) для практического вычисления размера шага интегрирования принимает вид
Автоматическое управление размером шага и порядком в неявных од-ношаговых экстраполяционных методах
Коротко, процедуру контроля можно изобразить в виде блок-схемы, представленной на рисунке 1.1. Необходимые детали практической реализации автоматического управления порядком и длиной шага в ГБШ-алгоритме описаны в [44, с. 244-247].
Этот алгоритм доказал свою высокую эффективность для явных экстраполяцион-ных методов. К недостаткам рассмотренной процедуры управления длиной шага и порядком следует отнести вычисление работы на шаге по формуле (1.3.15). Она проста и удобна, но вряд ли дает реальное представление о вычислительной работе при интегрировании дифференциальных уравнений. Более того, формулы (1.3.15) и (1.3.16), которые подсчитывают число вызовов функции для вычисления! правой части задачи (1.1.1а) за единичный шаг, в случае неявных методов работать недбудут. Точнее, если для явных экстраполяционных методов число вызовов функции для вычнеленнй правой части является приемлемой оценкой трудоемкости экстраполяции, то для неявных экстраполяционных методов это уже не так. Для последних, особенно в случае использования итерации ньютоновского типа при решении нелинейных задач, более важное значение приобретает решение возникающих линейных систем. Но с другой стороны, вычисление правой части задачи (1.1.1а) тоже может быть очень трудоемким и вносить существенных вклад в затраты машинного времени. Поэтому в следующей главе будет предложена другая концепция подсчета вычислительной работы за единичный шаг, которая не страдает от указанных недостатков.
Отметим также, что учет только главного члена асимптотического разложения локальной ошибки при выборе шага интегрирования (а на этой идее основано большинство известных процедур) не гарантирует малость всей локальной ошибки для выбранного шага и будет корректным только для разбиений с малым диаметром.
Другой проблемой почти всех известных методик автоматического управления порядком и шагом экстраполяционных методов (а также методов фиксированного порядка), включая описанный выше алгоритм Хайрера с соавт., является то, что все они осуществляют только контроль локальной ошибки, а значит — не позволяют контролировать глобальную погрешность, которая является более существенной для точности решения задачи. К тому же для явных методов (или других методов с ограниченной областью устойчивости) контроль локальной ошибки часто входит в противоречие с устойчивостью численного метода, что не только увеличивает затраты машинного времени, но и вообще может привести к непредсказуемым результатам.
Использование же неявных численных методов с автоматическим выбором шага интегрирования на практике существенно усложняется тем, что и реализация правила Рунге, и использование неявных численных методов разных порядков при решении нелинейных задач требуют значительных затрат машинного времени.
В итоге, локальный контроль не способен находить численное решение, удовлетворяющее наперед заданным требованиям к точности, в автоматическом режиме. Однако последнее свойство не только желательно, но и необходимо для целого ряда прикладных задач. В связи с этим, в вычислительной математике не раз предпринимались попытки оценивать и контролировать глобальную погрешность численных методов. Достаточно полный и хорошо систематизированный обзор всевозможных способов оценки глобальной ошибки численных методов, записанных в операторном виде, дал Скил [154]. Один из возможных подходов к оценке глобальной ошибки А-устойчивых методов представлен в [38]. РІнтересная идея минимизации глобальной ошибки на сетках с фиксированным количеством узлов описана в [4]. Но все эти результаты весьма далеки от практического использования, так как они отличаются либо большой трудоемкостью, либо малой эффективностью. Кроме того, не было разработано никаких более или менее реалистичных способов автоматического управления глобальной ошибкой. Так, тот же Скил в [155] отмечает, что наиболее популярный способ управления глобальной ошибкой состоит в предположении о пропорциональности локальной и глобальной ошибок с тем, чтобы на основе управления локальной ошибкой добиться необходимой малости глобальной. Однако там же указано, что даже если такая пропорциональность и имеет место, то коэффициент этой пропорциональности зависит от многих факторов, весьма специфичных для каждой конкретной задачи, учесть которые на практике не представляется возможным. Другой способ управления шагом интегрирования, предложенный в [38], тоже нельзя признать удовлетворительным, так как он не принимает во внимание противоречие между глобальным характером контролируемой величины (ошибки численного метода) и локальным характером управления (размера шага интегрирования). Практическое исследование работоспособности различных методик оценки глобальной ошибки описано в [53].
Дополнительные сложности возникают при использовании неявных методов для решения ОДУ. Как уже говорилось выше, в случае неявного метода мы не можем вычислить точно приближенное решение, даваемое этим методом, и вынуждены ограничиться некоторым итерационным приближением. Однако большинство алгоритмов оценки и контроля локальной и глобальной ошибок численных методов в основном базируются на их асимптотическом разложении. Наличие же неявностн предполагает, что в каждой точке сетки мы вынуждены применять итерационные алгоритмы, которые могут серьезно влиять на асимптотическое разложение погрешности численного решения, а значит, и на саму процедуру управления размером шага интегрирования, делая ее совершенно несостоятельной. Видимо, по этим причинам в современных и наиболее полных монографиях по численным методам для дифференциальных и дифференциально-алгебраических систем уравнений (см., например, [65], [100], [101]) вообще отсутствует какое-либо упоминание о способах автоматического контроля глобальной ошибки. Сравнительно недавно, Куликов разработал эффективный способ оценки глобальной ошибки неявных одношаговых численных методов и предложил практический алгоритм для использования таких оценок в процедуре автоматического управления точностью численного решения (см., например, [18, 21, 116, 118, 130]). Следуя идеям Грегга [93, 94], ему удалось получить асимптотическое разложение глобальной ошибки неявных одношаговых методов для решения дифференциальных и дифференциально-алгебраических уравнений индекса 1 при размере шага, стремящемся к нулю. Л также вывести уравнения для коэффициентов этого разложения (его главных членов), на основе которых построено несколько вариантов алгоритма локально-глобального контроля точности численного решения, которые наряду с локальной управляют и глобальной ошибкой. В результате, Куликов разработал эффективный метод для автоматического решения систем дифференциальных и дифференциально-алгебраических уравнений индекса 1 с любой наперед заданной точностью (в отсутствие ошибок округления). Позже ею методика вычисления глобальной ошибки и процедура управления этой оценкой были с успехом адаптированы к коллокационным одношаговым методам со старшими производными (см. [26]).
Рассмотрим здесь улучшенную версию алгоритма Куликова, т.е. устойчивый неявный локально-глобальный контроль точности, который основан на контроле двух главных членов в асимптотическом разложении глобальной ошибки (1.2.1) (см. [21, 121]). Компактно, этот алгоритм можно изобразить в виде блок-схемы, представленной на рисунке 1.2.
Эффективная реализация Е-методов и упрощенные итерации Ньютона
Во второй главе был рассмотрен механизм автоматического управления размером шага и порядком в экстраполяционных процессах. За основу мы взяли алгоритм Хайрера с соавт., предложенный для ГБШ-алгоритма, и показали как можно улучшить его производительность за счет более точного учета трудоемкости экстраполяции. Новая концепция работы базируется не на числе вызовов функции для вычисления правой части исходной задачи, а на учете реальных затрат процессорного времени, что позволяет проводить численное интегрирование максимально эффективно. Далее, мы адаптировали процедуру управления для неявных экстраполяционных методов и показали, что модифицированный алгоритм достаточно работоспособен. Однако, учитывая только величину локальной ошибки интегрирования, невозможно получить численное решение с наперед заданной точностью в автоматическим режиме, что является наиболее важным при решении задачи. Чтобы разрешить подобную проблему, в настоящей главе мы предложим другую методику управления размером шага и порядком, которая основана на контроле как локальной так и глобальной ошибок численного решения. Конечно, управление глобальной ошибкой требует дополнительных затрат машинного времени, но зато позволяет находить численное решение, удовлетворяющее наперед заданным требованиям к точности, в автоматическом режиме.
Для управления глобальной погрешностью неявных экстраполяционных методов воспользуемся методикой Куликова оценки глобальной ошибки одношаговых численных методов, а также его алгоритмом устойчивого неявного локально-глобального контроля точности, предложенными впервые в [18] и [118]. В параграфе 1.3 представлена улучшенная версия этого алгоритма, основанная на контроле двух главных членов в асимптотическом разложении глобальной ошибки (1.2.1) (см. блок-схему на рисунке 1.2).
Итак, целью настоящего параграфа является адаптация указанного алгоритма контроля локальной и глобальной ошибок к неявным экстраполяционным методам и повышение его эффективности за счет вычисления большего числа членов в разложении глобальной ошибки. Допустим, что мы уже вычислили два главных члена глобальной ошибки с необходимой точностью и выяснили, что размер шага Тк является слишком большим и не позволяет найти численное решение с оценкой глобальной ошибки, не превосходящей установленного предела єді0ь- Тогда исходная версия устойчивого локально-глобального контроля точности требует уменьшения размера шага тк, для увеличения точности численного интегрирования, и повторного вычисления всех найденных ранее величин. Для экстраполяционных методов это не является обязательным. Наша идея состоит в том, чтобы поднять точность интегрирования за счет увеличения порядка экстраполяционного метода для фиксированного размера шага г ;. Поэтому, если возможно, мы поднимаем порядок экстраполяционного метода на единицу (или на два в случае квадратичной экстраполяции) и вычисляем коэффициенты ips+q-2{tk+i) и ifts+q-i(tk+i) (после численного интегрирования исходной задачи экстраноляционным методом порядка s + q — 1 на локальном интервале [tk,ti+{\). Фактически, j-ые компоненты этих векторов есть ни что иное как последние две компоненты в векторе Y(j), деленные на размер шага тк в соответствующей степени. Конечно, на практике следует избегать указанного деления (а значит, и последующего умножения на то же самое число) через замену численной аппроксимации уравнения (1.2.2) относительно коэффициентов членов разложения глобальной ошибки, на разностное уравнение относительно самих членов асимптотического разложения глобальной ошибки (подробное разъяснение преимуществ такой замены можно найти в [21], [26] или [121]). Напомним, что векторы Y(j), j = l,2,...,m, получаются в результате применения экстраполяционного процесса как решения линейных систем (см. формулу (1.2.4), в случае обычной экстраполяции, или (1.2.26), в случае квадратичной экстраполяции, в параграфе 1.2). Здесь надо пояснить, что компоненты векторов Y(j), которые содержат члены асимптотического разложения глобальной ошибки базового РК-метода (1.1.2), могут интерпретироваться как компоненты, содержащие главный член локальной ошибки экстраполяционного метода при соответствующем выборе числа экстраполяции, при условии, конечно, что глобальная ошибка в предыдущей точке сетки tk равна нулю (в пределах установленной точности). В итоге, величины ips+q-2(tk+i) и s+q-i(tk+i) означают приближенные значения для неоднородных членов в уравнении (1.2.2) при i = s + q — 3ni = s + q- 2, т.е. для s+7_o(ifc+i) и Vs+f/-i(ifc+i), найденные с точностью по крайней мерс 0(тк). Затем мы интегрируем уравнение (1.2.2) при г = s+q — 3 и г = s + q —2 неявным методом Эйлера и получаем оценки s+(7_3(fc+1) и i/ s+(7_2(fc+i). Проверив погрешность этого интегрирования, мы осуществляем контроль точности приближенного решения с помощью двух последних членов в разложении глобальной ошибки, исключенных в результате применения экстраполяции, а именно — с помощью i/)s+7_3(fc+i) и s+a_2(t +i). Контроль точности численного решения задачи (1.2.2), особенно необходимый при решении жестких дифференциальных уравнений вида (1.1.1), осуществляется за счет повторного вычисления оценок для коэффициентов двух главных членов в разложении глобальной ошибки (1.2.1), но только методом второго порядка. Для этого рекомендуется применять либо неявное правило средней точки, либо правило трапеций. Разности этих коэффициентов, найденных численными методами первого и второго порядков, дают нам оценки локальных ошибок для решений уравнения (1.2.2) при г = s + q — 3 и г = s + q — 2, полученных неявным методом Эйлера. Эти оценки обозначены как Аф s+q_3(tk+x) и Ai/ s+9_2(ifc+1) на рисунке 3.1 Затем, если вычисление двух главных членов проведено достаточно точно, мы проверяем, что полученная оценка глобальной ошибки удовлетворяет заданному критерию точности. Под этим понимается, что как главный член в разложении глобальной ошибки, так и сумма двух первых членов не превосходят в некоторой норме наперед заданной точности вычислений. Напомним, что это организовано в виде проверки двойного равенства тк TV = тк, где тк означает максимальный размер постоянного шага для интегрирования на локальном отрезке [fc,fc+i] который обеспечивает выполнения условия \\ s+q-3{tk+\){Tk)s+q 2 Єдіоь. Размер шага т означает тоже самое, но для условия !І 3+г/-з( +і)(-г )8+7 3 + +(,_2( +і)(т )5+,3-2 е01оЪ. Итак, если найденная оценка глобальной ошибки удовлетворяет установленному критерию качества численного решения, то мы переходим к интегрированию задачи (1.1.1) в следующей точке сетки. Если нет, то увеличиваем порядок экстраполяционного метода еще раз. Процедура, описанного выше, улучшенного локально-глобального управления размером шага и уюрядком экстпраполяционных методов изображена на рисунке 3.1.
Подчеркнем, что мы контролируем глобальную ошибку численного метода порядка s+q—З, а интегрирование ведем на самом деле численным методом порядка s+q—1. Поэтому мы имеем некоторый запас точности для ошибок округления, ошибок итерационных методов, вовлеченных в неявную экстраполяцию, и ошибок, возникающих из-за усечения разложения (1.2.1) (мы пренебрегаем членами порядка s + q — 1 и выше). В квадратичной экстраполяции такой запас прочности составляет четыре порядка, что плодотворно сказывается на точности вычислений (см. ниже результаты численных экспериментов).
Локально-глобальное управление размером шага и порядком в одногаа-говых коллокационных методах со старшими производными
В параграфе 4.3 мы воспользуемся очень гибким алгоритмом одновременного контроля локальной и глобальной ошибок численного интегрирования, предложенным для неявных экстраполящгонных методов в параграфе 3.1 (см. блок-схему на рисунке 3.1).
Заметим, что теорема 1 (см. параграф 1.1), являющаяся теоретическим основанием для всей процедуры локально-глобального контроля, справедлива только для методов с постоянным шагом, поэтому ее адаптация к переменным сеткам требует определенных усилий. Так, в [26] можно найти методику вычисления и контроля главного члена глобальной ошибки применительно к Е-методам со старшими производными. Ранее, аналогичная процедура была использована для контроля глобальной ошибки в стандартных методах Рунге-Кутты (см. [21] или [18], [118]). В процитированных работах она носит название устойчивого нелепого локально-глобального контроля точности. Однако, свой окончательный вид локально-глобальный контроль точности численного решения приобрел в параграфе 3.1 с использованием экстраполяционной техники. Напомним, что его дополнительная гибкость состоит в том, что мы меняем как размер шага интегрирования, так и порядок численного метода. Компактно, этот механизм представлен в виде блок-схемы, изображенной на рисунке Прежде, чем приступить к построению процедуры управления, необходимо обсудить некоторые нюансы, связанные с использованием Е-методов со старшими производными. Во-первых, обсуждаемый локалыю-глобальный контроль точности базируется на экстраполяционной технике, а значит — порядок численного решения, вычисляемого в каждой точки сетки ifc+i, будет меняться в зависимости от числа q интегрирований на локальном интервале [fc,fc+i]- Очевидно, что минимальное число упрощенных итераций (4.2.2), которое ограничивалось снизу условием (4.2.16), должно теперь соответствовать точности численного решения и тоже зависеть от параметра q. Здесь следует вспомнить результаты параграфа 4.1, а именно, теорему 15, в которой доказано, что Е-методы (1.1.6) с коэффициентами (4.1.11) являются симметричными, а значит, могут быть использованы для построения квадратичной экстраполяции. Напомним, что квадратичная экстраполяция означает, что каждый новый столбец в экстраполяционной таблице увеличивает порядок экстраполяционного метода на два. Таким образом, порядок экстраполяционной схемы, основанной на Е-методах со старшими производными, вычисляется по формуле 2(p+q —1)+4. Тогда из оценки ошибки (4.2.11) можно получить условие для достаточного числа упрощенных
Во-вторых, численные эксперименты параграфа 3.2 подтвердили, что использование квадратичной экстраполяции в локально-глобальном контроле точности позволяет находить численное решение, соответствующее запросу пользователя, в автоматическом режиме и за короткое время. Следовательно, такие алгоритмы вычислений очень эффективны на практике, что является дополнительным аргументом для применения квадратичной экстраполяции в Е-методах со старшими производными. В-третьих, если мы возьмем в качестве достаточного числа упрощенных итераций (4.2.2) в каждой точке сетки правую часть неравенства (4.3.1), то это будет означать, что Е-метод со старшими производными и итерационный процесс (4.2.2) вносят соразмерный вклад в ошибку комбинированного метода, что не очень хорошо, так как алгоритм локально-глобального контроля точности учитывает только ошибку Е-метода и игнорирует ошибку итерационного процесса. Таким образом, чтобы сделать наш контроль точности асимптотически корректным, мы увеличиваем количество итераций в каждой точке сетки на единицу, т.е. в дальнейшем будем вычислять достаточное число итераций (4.2.2) по формуле
Отметим, что условие (4.3.2) является ключевым для надежного функционирования всей процедуры локально-глобального контроля точности переменного порядка. В-четвертых, работоспособность обсуждаемой процедуры контроля локальной и глобальной ошибок Е-методов со старшими производными зависит от трех параметров, устанавливаемых пользователем, а именно от Єіос — допуска на локальную ошибку комбинированного Е-метода, Єдш, — допуска на глобальную ошибку комбинированного Е-метода и Єід — допуска на локальную ошибку неявного метода Эйлера, примененного для приближенного решения уравнения (1.2.2) относительно коэффициентов двух главных членов разложения глобальной ошибки (1.2.1) (см. детали в параграфе 3.1). Для вычислительных экспериментов мы выбираем следующее соотношение между ними: Таким образом, достаточно задать один из этих допусков, чтобы определить их все. Ниже мы представим данные численных экспериментов в зависимости от параметра єді0ь, подразумевая, что равенство (4.3.3) всегда выполняется.
Далее, продемонстрируем на примерах работоспособность методов (1.1.6) с неявным локально-глобальным контролем точности численного решения переменного порядка.
Итак, для тестирования наших численных методов выберем те же самые задачи, что и в параграфах 2.3 и 3.2, т.е. (2.1.1), (2.1.4) и (2.1.5). Основное достоинство указанных задач состоит в том, что мы знаем особенности поведения их точного решения (см. параграф 2.1), а значит — можем вычислить точные ошибки с целью проверки надежности процедуры управления глобальной ошибкой Е-методов со старшими производными. Напомним также, что задача (2.1.5) представляет собой пример с негладким решением и позволяет проверить работоспособность Е-методов с автоматическим контролем глобальной ошибки в такой ситуации.
Для реализации процедуры локально-глобального контроля возьмем Е-методы 4-го, 6-го и 8-го порядков из параграфа 4.2, заданные таблицами коэффициентов (4.2.17)-(4.2.19).
Итак, дополняя комбинированные Е-методы (4.2.17)-(4.2.19) с упрощенными итерациями Ньютона (4.2.2) процедурой локально-глобального контроля точности переменного порядка из параграфа 3.1 и применяя их к задачам (2.1.1), (2.1.4) и (2.1.5), получаем следующие результаты. В таблицах 4.4-4.6 представлены данные численного интегрирования задачи (2.1.1) Е-методами (4.2.17)-(4.2.19), соответственно. Точнее, там показаны точность вычислений и затраты процессорного времени в зависимости от установленного допуска на глобальную ошибку численного метода, т.е. в зависимости от є дм. Напомним, что остальные параметры процедуры локально-глобального контроля точности (єіос п Єід) подчиняются условию (4.3.3). Дополнительно отметим, что, как и во всей работе, численные эксперименты проводились на персональном компьютере с процессором Intel Celeron., 2233MHz, а программы написаны на языке Borland C++, Release 5.02. Аналогично, в таблицах 4.7-4.9 даны результаты для решения задачи (2.1.4), а в таблицах 4.10-4.12 — для численного интегрирования задачи (2.1.5). Во всех случая удалось решить тестовые задачи с погрешностью, не превосходящей установленную границу є9[0ь что подтверждает работоспособность алгоритма управления глобальной ошибкой одношаговых методов из параграфа 3.1 применительно к Е-методам со старшими производными, обсуждаемыми в этой главе. Подчеркнем отдельно, что требуемая точность вычислений для всех задач достигается за весьма незначительное время, что тоже указывает на серьезный прикладной потенциал Е-методов с автоматическим контролем глобальной ошибки.
Напомним, что процедура локально-глобального контроля точности основана на комбинированном управлении размером шага и порядком базового метода. За счет этого удается повысить гибкость управления, а значит — сократить затраты машинного времени. Покажем, как эта методика работает на практике для достижения поставленной цели. Точнее, рисунки 4.1-4.3 демонстрируют изменение размера шага численного интегрирования и его порядка в динамике с целью контроля погрешности численного решения для задач (2.1.1), (2.1.4) и (2.1.5), соответственно. Верхние два графика на каждой фигуре дают информацию для Е-метода (4.2.17), средние графики — для Е-метода (4.2.18) и нижние графики — для Е-метода (4.2.19). На каждом рисунке изображено изменение размера шага и порядка в процессе численного решения
