Коллокационные методы со старшими производными и неявная экстраполяция численного решения
на правах рукописи
Хрусталева Екатерина Юрьевна КОЛЛОКАЦИОННЫЕ МЕТОДЫ СО СТАРШИМИ ПРОИЗВОДНЫМИ И НЕЯВНАЯ ЭКСТРАПОЛЯЦИЯ ЧИСЛЕННОГО РЕШЕНИЯ Специальность: 05.13.18 математическое моделирование, численные методы и комплексы программ
Автореферат диссертации на соискание учной степени е кандидата физико-математических наук
Ульяновск – 2010
Работа выполнена на кафедре Информационных технологий в Госу дарственном образовательном учреждении высшего профессионального образования Ульяновский государственный университет.
Научный консультант: доктор физико-математических наук, доцент Куликов Геннадий Юрьевич
Официальные оппоненты: доктор физико-математических наук, профессор Горбунов Владимир Константинович доктор физико-математических наук, профессор Кузнецов Евгений Борисович
Ведущая организация: ГОУ ВПО Ульяновский государственный технический университет
Защита состоится 24 ноября 2010 года в 1130 часов на заседании диссертационного совета Д 212.278.02 при Ульяновском государственном университете по адресу: г. Ульяновск, ул. Набережная р. Свияги, 106, корп. 1, ауд. 703.
С диссертацией можно ознакомиться в научной библиотеке Ульянов ского государственного университета, с авторефератом – на сайте ВУЗа http://www.uni.ulsu.ru.
Автореферат разослан 2010 года.
Просьба направлять отзыв на автореферат в одном экземпляре, заверенном печатью организации по адресу: 432000, г. Ульяновск, ул. Л.Толстого, д. 42, УлГУ, Управление научных исследований.
Ученый секретарь диссертационного совета Волков М.А.
Общая характеристика работы
Актуальность темы. Системы обыкновенных дифференциальных уравнений (ОДУ) представляют собой неотъемлемую часть математи ческого моделирования динамических объектов. Однако размер и слож ность современных моделей делает аппарат аналитических вычислений практически неприемлемым в данной ситуации. Поэтому проблеме чис ленного интегрирования задач вида x (t) = g(t, x(t)), t [0, T ], (1а) x(0) = x0, (1б) где правая часть g : Rm+1 Rm является достаточно гладкой функци ей, отводится существенное место в вычислительной математике. Такие задачи возникают как непосредственно при математическом моделирова нии в области биологии, медицины, механики, электротехники, химии и т. д., так и при решении более сложных систем уравнений (например, при дискретизации задач математической физики методом прямых).
В настоящее время к наиболее перспективным алгоритмам числен ного решения задачи Коши для ОДУ можно отнести экстраполяци онные методы. Общая теория экстраполяционных методов, использую щих неявные одношаговые схемы, была впервые предложена Кулико вым Г.Ю.1 Наряду с хорошо изученными стандартными одношаговыми методами Рунге-Кутты в качестве опорных для построения экстрапо ляционных схем, можно предложить сравнительно новые одношаговые методы Рунге-Кутты, использующие старшие производные (см., напри мер, Аульченко С.М., Латыпов А.Ф., Никуличев Ю.В.2, Куликов Г.Ю., Kulikov G.Yu. On implicit extrapolation methods for ordinary dierential equations// Russian Journal of Numerical Analysis and Mathematical Modelling. 2002. V. 17, № 1. P. 41–69.
Аульченко С.М., Латыпов А.Ф., Никуличев Ю.В. Метод численного интегрирования систем обыкновенных дифференциальных уравнений с использованием интерполяци онных полиномов Эрмита// Журнал вычислительной математики и математической физики. 1998. Т. 38, № 10. С. 1665–1670.
Меркулов А.И.3, Fehlberg E. 4, Kastlunger K.H.5, Nrsett S.P.6 ). Однако, исследования в этой области пока достаточно ограничены, а определение условий симметричности таких методов и использование их в качестве базовых для построения экстраполяции не изучалось вовсе.
Далее следует отметить, что важнейшим этапом в эффективной ре ализации любого численного метода для интегрирования обыкновенных дифференциальных уравнений является разработка механизма автомати ческого управления размером шага с целью обеспечения требуемой точно сти вычисления за минимальное (или приемлемое) время. Дополнитель ным преимуществом экстраполяционных методов является возможность автоматического подбора оптимальных размера шага и порядка в каждой точке сетки. Вопрос о надежной и достоверной оценке точности числен ного решения ОДУ неразрывно связан с проблемой оценки глобальной ошибки, которой посвящено достаточно большое число работ. В то вре мя как обсуждению алгоритмов автоматического контроля глобальной ошибки уделено значительно меньше внимания. Можно указать лишь от носительно небольшое количество статей, так или иначе затрагивающих эту тематику: Новиков Е.А.7, Dahlquist G.8, Lindberg B.9, Skeel R.D.10, Stetter H.J.11. А важное преимущество экстраполяционных методов, поз воляющее управлять не только размером шага, но и порядком самого Куликов Г.Ю., Меркулов А.И. Об одношаговых коллокационных методах со стар шими производными для решения обыкновенных дифференциальных уравнений// Журнал вычислительной математики и математической физики. 2004. Т. 44, № 10. С. 1782–1807.
Fehlberg E. New high-order Runge-Kutta formulas with step size control for systems of rst and second order dierential equations// ZAMM. 1964. V. 44, Sonderheft P.T17–T19.
Kastlunger K.H., Wanner G. Runge-Kutta processes with multiple nodes// Computing. 1972. V. 9. P. 9–24.
Nrsett S.P. One-step methods of Hermite type for numerical integration of sti systems// BIT. 1974. V. 14. P. 63–77.
Новиков Е.А. Оценка глобальной ошибки A-устойчивых методов решения жестких систем// Доклады академии наук. 1995. Т. 343, № 4. С. 452–455.
Dahlquist G. On the control of the global error in sti initial value problems// Lecture Notes in Mathematics. Berlin-New York: Springer, 1982. V. 912. P. 38–49.
Lindberg B. Characterization of optimal stepsize sequences for methods for sti dierential equations// SIAM Journal on Numerical Analysis. 1977. V. 14, № 5.
P. 859–887.
Skeel 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.
метода, практически не используется.
Таким образом, объектом исследования являются неявные одноша говые экстраполяционные схемы, а предметом исследования – алгорит мы комбинированного управления шагом и порядком экстраполяционных процессов, позволяющие находить численное решение, удовлетворяющее заданным требованиям к точности, в автоматическом режиме.
Цель и задачи работы. Целью работы является построение и тести рование новых численных методов, использующих старшие производные и экстраполяционную технику, для решения обыкновенных дифференци альных уравнений. А также разработка теории комбинированного управ ления размером шага численного интегрирования и порядком одноша говых методов (в том числе и коллокационных методов со старшими производными) с использованием экстраполяций, позволяющей достичь наперед заданную точность вычислений в автоматическом режиме. В со ответствии поставленной целью решались следующие задачи:
• разработка алгоритмов неявного локально-глобального управления порядком и шагом экстраполяционных методов;
• исследование свойств одношаговых методов Рунге-Кутты, использу ющих старшие производные;
• эффективная реализацию разработанных в диссертации алгоритмов с использованием новых упрощенных итераций Ньютона;
• разработка комплекса программ для численного интегрирования си стем обыкновенных дифференциальных уравнений с заданной точ ностью в автоматическом режиме.
Методы исследования. В диссертационной работе использованы ме тоды вычислительной математики, математического и функционального анализа, линейной алгебры и программирования.
Научная новизна. Все основные результаты настоящей диссертации являются новыми и состоят в следующем:
1. Разработаны алгоритмы контроля локальной и глобальной оши бок неявных экстраполяционных методов для систем обыкновенных дифференциальных уравнений.
2. Построены экстраполяционные схемы, использующие коллокацион ные методы со старшими производными, для нахождения численного решения систем обыкновенных дифференциальных уравнений и изу чены их фундаментальные свойства.
3. Разработаны алгоритмы контроля локальной и глобальной ошибок экстраполяционных методов, использующих коллокационные мето ды со старшими производными, для систем обыкновенных диффе ренциальных уравнений.
4. Предложена эффективная схема упрощенных итераций Ньютона, оп тимизирующих реализацию разработанных в диссертации алгорит мов.
Основные положения, выносимые на защиту:
1. Методика построения присоединенных и симметричных коллокаци онных методов со старшими производными, формулы для вычисле ния коэффициентов таких методов.
2. Способы автоматического контроля локальной и глобальной ошибок экстраполяционных методов, использующих неявные одношаговые методы Рунге-Кутты и коллокационные методы со старшими произ водными, для систем обыкновенных дифференциальных уравнений.
3. Методика построения упрощенных итераций Ньютона, доказатель ство сходимости комбинированных схем неявных экстраполяцион ных методов с использованием нового итерационного процесса, вы вод асимптотической оценки погрешности.
4. Комплекс программ на языке C++, позволяющий интегрировать си стемы обыкновенных дифференциальных уравнений с заданной точ ностью в автоматическом режиме.
Достоверность результатов. Достоверность результатов обеспечи вается строгостью постановок задач и доказательств теорем, использова нием аналитических и численных методов расчета, методов математиче ского моделирования и применением современных методик эксперимен тального исследования.
Практическая и теоретическая значимость. Полученные в дис сертации новые результаты представляют интерес для численного реше ния систем ОДУ экстраполяционными методами с контролем локальной и глобальной ошибок. В прикладном аспекте такие методы могут с успе хом быть использованы для разработки автоматизированных систем ма тематического моделирования в различных областях науки и техники, поскольку способность обеспечивать заданную пользователем точность вычислений в автоматическом режиме существенно повышает достовер ность математического моделирования на практике.
Апробация работы. Результаты диссертации докладывались и об суждались: на пятой международной научно-технической конференции Математическое моделирование физических, экономических, техниче ских, социальных систем и процессов (Ульяновск, 2003), на XXVI кон ференции молодых ученых механико-математического факультета МГУ им. М.В. Ломоносова (Москва, 2004), на 4th International Conference on Computer Science (Krakow, Poland, 2004).
Личный вклад автора. Постановка задач осуществлялась совмест но с научным руководителем, д.ф.-м.н. Г.Ю. Куликовым, теоретические положения, проведение расчетов и анализ полученных результатов вы полнены автором самостоятельно.
Публикации. По теме диссертации опубликовано 10 работ, в том чис ле 3 в рецензируемых научных изданиях, рекомендованных ВАК, их спи сок помещен в конце автореферата.
Структура диссертации. Диссертация состоит из введения, четы рех глав, заключения, списка литературы (162 наименования) и прило жения. Общий объем диссертации составляет 128 страниц, из них страниц основного текста, 10 страниц приложения. Работа включает таблиц и 20 рисунков.
Содержание работы Первая глава, состоящая из 3 параграфов, носит вводный вспомога тельный характер и содержит обзор литературы, посвященной как раз личным аспектам теории экстраполяционных методов для решения ОДУ, так и вопросам получения и практического использования оценок локаль ной и глобальной ошибок с целью управления размером шага и порядком численного интегрирования.
Вторая глава, состоящая из трех параграфов, посвящена неявным экстраполяционным методам и проблеме локального контроля точности.
В параграфе 2.1 рассматривается модификация известного алгоритма управления порядком и шагом в явных экстраполяционных методах12.
Предлагается улучшить производительность процедуры управления за счет более точного учета трудоемкости экстраполяции. Новая концепция работы на шаге базируется не на числе вызовов функции для вычис ления правой части исходной задачи, а на подсчете реальных затрат про цессорного времени, что позволяет проводить численное интегрирование максимально эффективно. Идея оптимального выбора порядка состоит в минимизации затрат машинного времени, необходимых на вычисление каждого шага экстраполяционного метода.
Данные вычислительных экспериментов, приведенные в параграфе 2.1, говорят о том, что модифицированное (т.е. с новой идей для вычис лении работы на шаге) автоматическое управление длиной шага и поряд ком ГБШ-алгоритма действительно работает на практике и позволяет находить численное решение за очень короткое время и с приемлемой точностью. Причем, при сравнении результатов работы стандартного и модифицированного алгоритмов, мы наблюдаем значительное увеличе ние эффективности последнего за счет более точного учета трудоемкости экстраполяции. Дополнительным преимуществом предложенной модифи кации является тот факт, что алгоритм Хайрера с соавт. с новой концеп цией вычисления работы может быть успешно адаптирован к неявным экстраполяционным методам, которые более устойчивы, чем явные про цессы.
Далее, в параграфе 2.2, производится адаптация описанной выше про цедуры управления для неявных экстраполяционных методов и показы вается, что модифицированный алгоритм достаточно работоспособен и в этом случае. Неявные экстраполяционные методы представляют осо бый интерес с практической точки зрения, поскольку именно они обла дают свойствами, необходимыми для квадратичной экстраполяции, кото рая значительно эффективней обычной. Однако при реализации следует учесть невозможность нахождения точного численного решения, опреде ляемого одношаговым методом. В этом случае применение стандартного алгоритма Хайрера с соавт. попросту неэффективно, так как число вы зовов функции для вычисления правой части уже не будет достоверно отражать трудоемкость алгоритма, поскольку для выполнения каждого шага метода необходимо использовать итерационные приближения. Одна Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных урав нений. Нежесткие задачи. М.: Мир, 1990.
ко идея для вычисления работы на шаге, предложенная в параграфе 2.1, остается работоспособной и в этом случае, так как мы используем такой универсальный критерий трудоемкости как затраты машинного времени.
Результаты практической реализации указанной процедуры контроля, приведенные в параграфе 2.3, однозначно подтверждают работоспособ ность рассмотренного способа автоматического управления длиной шага и порядком неявных экстраполяционных методов. Здесь же отмечается, что, учитывая только величину локальной ошибки, невозможно получить численное решение с наперед заданной точностью в автоматическим ре жиме, что является наиболее важным при решении задачи.
Третья глава, состоящая из двух параграфов, занимает важное ме сто в диссертационной работе, она посвящена проблеме контроля точно сти численного решения систем дифференциальных уравнений на основе автоматического управления размером шага интегрирования и порядком в неявных одношаговых экстраполяционных методах.
В параграфе 3.1 предлагается новая концепция управления размером шага и порядком экстраполяционного метода, в основе которой лежит контроль как локальной, так и глобальной ошибки численного решения.
Конечно, управление глобальной ошибкой требует дополнительных за трат машинного времени, но зато позволяет находить численное решение, удовлетворяющее наперед заданным требованиям к точности, в автома тическом режиме.
В параграфе 3.1 представлена улучшенная версия алгоритма устойчи вого неявного локально-глобального контроля точности, адаптирован ная к неявным экстраполяционным методам и основанная на контроле большего числа членов в асимптотическом разложении глобальной ошиб ки. В основе этого алгоритма лежит методика Куликова13 оценки гло бальной ошибки одношаговых численных методов.
Итак, введем на отрезке [0, T ] произвольную неравномерную сетку w = {tk+1 = tk + k, k = 0, 1,..., K 1, t0 = 0, tK = T } и обозначим через диаметр сетки w, т.е. = max {k }.
0kK Тогда одношаговый метод вида Kulikov G.Yu. A local-global version of a stepsize control for Runge-Kutta methods// Korean Journal of Computational and Applied Mathematics. 2000. V. 7, № 2.
P. 289–318.
xk k, Niter c E Находим xk+1, s+q2 (tk ), s+q1 (tk ) E Вычисляем k 3c 33 33 T T T T Нет 3 s+q2 (tk ) k s+q 3 loc Да c Находим k+1,k+2,s+q2 (tk+1),s+q1 (tk+1) x c Нет Вычисляем s+q3 (tk+1 ) 33 3 33 s+q3 ' 3 s+q3 (tk+1 ) k lg 333 c Вычисляем s+q3 (tk+1 ) Да c E Находим k Вычисляем s+q2 (tk+1 ) 33 Нет 3 33 s+q2 ' 3 s+q2 (tk+1 ) k lg c Да Вычисляем s+q2 (tk+1 ) c Находим k cr r r k = k = rr 33 k rr 33 rr r Нет 3 33 Нет loc ' s+q2 s+q2 (tk+1 ) k+1 33 Да 33 c Принимаем xk+1, k+1 ;
Да Bычисляем Niter c xk+1 k+1, Niter q + +;
Вычисляем Niter c Рис. 1: Схема неявного локально-глобального контроля точности переменного порядка tki = tk + k ci, (2а) l xki = xk + k aij g(tkj, xkj ), i = 1, 2,..., l, (2б) j= l xk+1 = xk + k bi g(tki, xki ), k = 0, 1,..., K 1, (2в) i= где x0 = x, а aij, bi и ci – вещественные коэффициенты, называется l-стадийным методом Рунге-Кутты (РК-методом).
Известно, что если интегрирование исходной задачи (1) велось с по стоянным шагом k, то глобальная ошибка метода (2) имеет следующее асимптотическое представление в виде ряда по степеням размера этого шага:
s s+ x(tk+1 ) xk+1 = s (tk+1 )k + s+1 (tk+1 )k +... + (3) S + S (tk+1 )k + k, k = 0, 1,..., K 1, S+ где k = O(k ) и i (t), i = s, s + 1,..., S, суть решения задач Коши i (t) = x g(t, x(t))i (t) + i+1 (t), i (0) = 0, (4) в которых i+1 (t) означает коэффициент главного члена локальной ошиб ки некоторого одношагового метода. При i = s этот член, т.е. s+1 (t), совпадает с коэффициентом главного члена локальной ошибки исходного РК-метода (2).
Суть представленного на рисунке 1 контроля точности заключается в следующем. Допустим, что мы уже вычислили два главных члена гло бальной ошибки с необходимой точностью и выяснили, что размер шага k является слишком большим и не позволяет найти численное решение с оценкой глобальной ошибки, не превосходящей установленного предела glob. Тогда новая концепция оптимального порядка состоит в том, чтобы поднять точность интегрирования за счет увеличения порядка экстрапо ляционного метода для фиксированного размера шага k. Поэтому, если возможно, мы поднимаем порядок экстраполяционного метода на едини цу (или на два в случае квадратичной экстраполяции) и вычисляем коэф фициенты старших членов локальной ошибки s+q2 (tk+1 ) и s+q1 (tk+1 ) (после численного интегрирования исходной задачи экстраполяционным методом порядка s + q 1 на локальном интервале [tk, tk+1 ]). Здесь необ ходимо отметить, что члены асимптотического разложения глобальной ошибки базового РК-метода (2), могут интерпретироваться как компо ненты, содержащие главный член локальной ошибки экстраполяционного метода при соответствующем выборе числа экстраполяций, при условии, конечно, что глобальная ошибка в предыдущей точке сетки tk равна нулю (в пределах установленной точности). В итоге, величины s+q2 (tk+1 ) и s+q1 (tk+1 ) означают приближенные значения для неоднородных членов в уравнении (4) при i = s + q 3 и i = s + q 2, т.е. для s+q2 (tk+1 ) и s+q1 (tk+1 ), найденные с точностью по крайней мере O(k ). Далее, если найденные величины удовлетворяют локальному контролю точности (не превосходят установленного предела loc ), численно интегрируем диффе ренциальное уравнение (4) неявным методом Эйлера с целью нахожде ния достаточно точных оценок для коэффициентов двух первых членов в разложении глобальной ошибки (3). Проверив погрешность этого инте грирования, мы осуществляем контроль точности приближенного реше ния с помощью двух последних членов в разложении глобальной ошибки, исключенных в результате применения экстраполяции. Контроль точно сти численного решения задачи (4), особенно необходимый при решении жестких дифференциальных уравнений вида (1), осуществляется за счет повторного вычисления оценок для коэффициентов двух главных членов в разложении глобальной ошибки (3), но только методом второго поряд ка (рекомендуется применять либо неявное правило средней точки, либо правило трапеций). Разности этих коэффициентов, найденных числен ными методами первого и второго порядков, дают нам оценки локальных ошибок для решений уравнения (4) при i = s+q3 и i = s+q2, получен ных неявным методом Эйлера. Эти оценки обозначены как s+q3 (tk+1 ) и s+q2 (tk+1 ) на рисунке 1. Затем, если вычисление двух главных чле нов проведено достаточно точно, мы проверяем, что полученная оценка глобальной ошибки удовлетворяет заданному критерию точности. Под этим понимается, что как главный член в разложении глобальной ошиб ки, так и сумма двух первых членов не превосходят в некоторой нор ме наперед заданной точности вычислений (glob ). Итак, если найденная оценка глобальной ошибки удовлетворяет установленному критерию ка чества численного решения, то мы переходим к интегрированию задачи (1) в следующей точке сетки. Если нет, то увеличиваем порядок экстра поляционного метода еще раз. При этом, порядок экстраполяционного метода в каждой точке сетки фиксируется как минимальный обеспечи вающий заданную точность вычислений. Если повышение порядка далее невозможно, то мы достигаем необходимую точность за счет уменьшения размера шага.
Подчеркнем, что мы контролируем глобальную ошибку численного ме тода порядка s + q 3, а интегрирование ведем на самом деле численным методом порядка s + q 1. Поэтому мы имеем некоторый запас точно сти для ошибок округления, ошибок итерационных методов, вовлечен ных в неявную экстраполяцию, и ошибок, возникающих из-за усечения разложения (3) (мы пренебрегаем членами порядка s + q 1 и выше). В квадратичной экстраполяции такой запас прочности составляет четыре порядка, что плодотворно сказывается на точности вычислений. Обратим внимание, что текущий порядок экстраполяционного метода, а значит и значение параметра q, устанавливающего количество численных реше ний исходной задачи (1) на локальном интервале [tk, tk+1 ], определяется в каждой точке сетки автоматически. Естественным ограничением сверху на порядок неявного экстраполяционного метода является число итера ций, которое мы выполняем при интегрировании неявным РК-методом (2). Поэтому, в случае если поднять точность интегрирования экстраполя ционным методом уже нельзя из-за погрешности итерационного процесса, вовлеченного в реализацию базового РК-метода (2), а обеспечить задан ную точность вычислений для размера шага k еще не удалось, то мы проводим повторное интегрирование исходной задачи (1) на локальном интервале [tk, tk+1 ], увеличив число итераций Niter на единицу.
Заметим, что приведенный выше алгоритм автоматического управле ния размером шага и порядком экстраполяции будет работать и в явных методах, рассмотренных в параграфе 2.1. Однако контроль глобальной ошибки в таким методах будет приводить к значительным дополнитель ным затратам машинного времени. Для неявных же экстраполяционных методов дополнительные затраты компьютерного времени будут относи тельно небольшими, так как реализация неявной экстраполяции и вы числение оценки глобальной ошибки требуют решения линейных систем с матрицами коэффициентов, задаваемыми матрицей Якоби исходной задачи (1). Таким образом, информация, необходимая для вычисления глобальной ошибки, доступна (и не требует дополнительных затрат) при использовании итераций ньютоновского типа.
Итак, улучшенный локально-глобальный контроль точности численно го решения требует от пользователя задания максимального количества экстраполяций (qmax ), а также трех параметров: loc допуск на локаль ную ошибку численного метода, glob допуск на глобальную ошибку численного метода и lg допуск на локальную ошибку неявного метода Эйлера, примененного для приближенного решения уравнения (4) отно сительно коэффициентов двух главных членов разложения глобальной ошибки (3).
Далее, в параграфе 3.2 представлено всестороннее тестирование раз работанной процедуры управления размером шага и порядком на диф ференциальных уравнениях с известным решением и анализ результа тов экспериментов. Во всех экспериментах установлена следующую связь между вышеуказанными допусками: loc = glob и lg = glob /10. Эти со отношения используются для того, чтобы показать какую пользу можно извлечь из дополнительного контроля глобальной ошибки при условии, что ограничения наложенные на локальную и глобальную погрешности численного решения идентичны. Таким образом, достаточно задать один из этих допусков, чтобы определить их все.
Четвертая глава, состоящая из 4 параграфов, содержит основные результаты диссертационной работы, полученные в классе одношаговых коллокационных методов со старшими производными, т.е. в ней теория комбинированного управления порядком и шагом расширена на упомя нутые выше вычислительные схемы.
В параграфе 4.1 обсуждаются свойства присоединенности и симметрии для РК-методов со старшими производными общего вида.
Определение 1 Одношаговый метод вида tki = tk + ci k, (5а) pj l (r) k aij g (r) (tkj, xkj ), r xki = xk + k i = 1, 2,..., l, (5б) j=1 r= pj l (r) k bj g (r) (tkj, xkj ), r xk+1 = xk + k k = 0, 1,..., K 1, (5в) j=1 r= где x0 = x0, g (r) (tkj, xkj ) обозначает r-ую производную правой части за дачи (1) по независимой переменной t, вычисленную в точке (tkj, xkj ), (r) (r) k размер k-го шага, а коэффициенты aij, bj и cj заданы таблицей:
(0) (1) (p ) (0) (1) (p ) c1 a11 a11... a111... a1l a1l... a1l l (0) (1) (p ) (0) (1) (p ) c2 a21 a21... a211... a2l a2l... a2l l................
.......,....... (6) (0) (1) (p1 ) (0) (1) (pl ) cl al1 al1... al1... all all... all (0) (1) (p ) (0) (1) (pl ) b1 b1... b1 1... bl bl... bl в которой A вещественная матрица размера l (p1 +p2 +...+pl +l), а b и c вещественные вектора размера (p1 + p2 +... + pl + l) и l соответ ственно, i, j = 1, 2,..., l, r = 0, 1,..., pj, причем (cj [0, 1]), называется l-стадийным методом Рунге-Кутты со старшими производными.
Известно, что симметричные одношаговые методы определяются с по мощью присоединенных. Таким образом, необходимо определить коэффи циенты присоединенного метода для РК-формулы (5), содержащей про изводные высоких порядков.
Теорема 1 Пусть метод (5) представляет собой l-стадийный РК метод со старшими производными, коэффициенты которого заданы таблицей (6). Тогда присоединенным для метода (5) будет также l стадийный РК-метод со старшими производными и коэффициентами (r) (r) aij, bj, c вида i c = 1 cl+1i, (7а) i (r) (r) (r) = (1)r (bl+1j al+1i,l+1j ), aij (7б) (r) (r) = (1)r bl+1j, bj r = 0, 1,..., pj, i, j = 1, 2,..., l. (7в) По определению, симметричным называется такой одношаговый ме тод, который совпадает со своим присоединенным. Применяя эту идею к РК-методам со старшими производными (5), мы заключаем, что для сим метрии достаточно потребовать, чтобы коэффициенты метода удовлетво ряли следующим равенствам:
ci = 1 cl+1i, (8а) (r) (r) (r) aij = (1)r (bl+1j al+1i,l+1j ), (8б) (r) (r) bj = (1)r bl+1j, r = 0, 1,..., pj, i, j = 1, 2,..., l. (8в) Подчеркнем, что формулы (8) необходимы для симметрии несократимого РК-метода со старшими производными. Кроме того, условия (8) подразу мевают выполнение следующего равенства:
pj = pl+1j, j = 1, 2,..., l. (9) Таким образом, условие (9) становится естественным требованием для симметричных РК-формул со старшими производными.
Симметричные методы являются базой для квадратичной экстраполя ции, которая весьма эффективна на практике. Поэтому в дальнейшем нам потребуются примеры симметричных РК-методов со старшими производ ными. Перейдем к семейству РК-методов со старшими производными, по строенному Куликовым и Меркуловым в процитированной выше статье.
Эти методы получили название Е-методов со старшими производными:
p (r) (r) (r) (r) (0) (0) r xk+1/2 = xk + k k (a1 gk + a3 gk+1 ) + k a2 gk+1/2, (10а) r= p (r) (r) (r) (r) (0) (0) r xk+1 = xk + k k (b1 gk + b3 gk+1 ) + k b2 gk+1/2, (10б) r= (r) k = 0, 1,..., K 1, где x0 = x0, gk+i = g (r) (tk+i, xk+i ), i = 0, 1/2, 1 и tk+i = tk + ik.
Чтобы Е-метод сходился с максимальным порядком 2p + 4, коэффици енты метода (10) должны удовлетворять следующим соотношениям:
(1)l (i + r)!
p + 1 pr i+r p+ (r) a1 = r!2p+r+2 i=0 l=0 j=0 l!(i + r l)!j!(p + 1 j)!(l + j + 2) (11а) i (p + q)!
, r = 0, 1,..., p, q!2q q= (1)l (p + 1)! p+ (0) (11б) a2 =, 2 l=0 l!(p + 1 l)!(2l + 1) (1)r+1 (p + 1) pr i+r p+1 (1)j (i + r)!
(r) a3 = r!2p+r+2 i=0 l=0 j=0 l!(i + r l)!j!(p + 1 j)!(l + j + 2) i (p + q)!
, r = 0, 1,..., p, q!2q q= (11в) (r) (r) (r) (0) (0) (r) (r) (r) b1 = a1 + (1)r a3, b3 = (1)r a1 + a3, b2 = 2a2, (11г) В настоящей работе предлагается некоторое упрощение выражений для вычисления коэффициентов Е-методов вида (10).
(r) (r) Теорема 2 Коэффициенты a1 и a3 коллокационного метода (10) со старшими производными вычисляются по формулам p + 1 pr p+1 i (p + q)!
(i + r)!(l + 1) (r) a1 =, (12а) r!2p+r+2 i=0 l=0 (p + 1 l)!(i + r + l + 2)! q=0 q!2q (1)r+1 (p + 1) pr i+r i (p + q)!
(i + r)!(l + 1) (r) a3 =, (12б) r!2p+r+2 (i + r l)!(p + l + 3)! q=0 q!2q i=0 l= где r = 0, 1,..., p.
Важным следствием параграфа 4.1 является доказательство того, что все Е-методы со старшими производными симметричны, а значит при годны для квадратичной экстраполяции.
Теорема 3 Е-метод (10) с коэффициентами (11)–(12) является сим метричным для любого целого p 0.
Заметим, что, как и любой неявный РК-метод, Е-методы со старши ми производными требуют привлечения дополнительных итерационных процессов для практического нахождения численного решения, а число и тип итераций существенным образом влияют как на порядок сходимости, так и на эффективность вычислений.
Параграф 4.2 посвящен сокращению дополнительных вычислитель ных затрат, возникающих при использовании неявной экстраполяции. В настоящей работе предлагается новая наиболее эффективная схема упро щенного метода Ньютона для систем нелинейных уравнений, возникаю щих при дискретизации методами вида (10) (в смысле затрат машинного времени).
Итак, введем следующие обозначения: x(tk+1 ) значение точного ре шения задачи (1) в точке tk+1, xk+1 значение точного решения нелиней ной задачи (10), а xk+1 = xk+1 (N ) значение приближенного решения задачи (10) в точке tk+1, полученное после N итераций некоторого ите рационного процесса. Очевидно, что применение Е-методов со старшими производными приводит на практике к следующей дискретной задаче:
k Xk+1 = G Xk+1, k = 0, 1,..., K 1 (13) относительно вектора неизвестных T R2m, Xk+1 = (xk+1/2 )T, (xk+1 )T а отображение G : R2m R2m задано формулой k p r (r) (r) (r) (r) (0) (0) xk + k k (a1 gk + a3 gk+1 ) + k a2 gk+1/ k G Xk+1 = r=, p (r) (r) (r) (r) (0) (0) b3 gk+1 ) + k b2 gk+1/ r xk + k k (b1 gk + r= где p означает наперед заданное неотрицательное целое число, а коэффи (r) (r) циенты ai, bj, i, j = 1, 2, 3, r = 0, 1,..., p, удовлетворяют соотношени ям (11)–(12). Упрощенные итерации Ньютона для нелинейной задачи (13) запишем следующим образом:
Xk+1 = Xk+1 A(Xk+1 )1 Fk Xk+1, 1 = 1, 2,..., N, (14а) Xk+1 = (T, xT )T R2m, k = 0, 1,..., K 1, xk k (14б) где Fk = E2m G (E2m единичный оператор в R2m ) и k (0) (0) Em a2 x g(tk+1/2, x0 a3 x g(tk+1, x0 ) k+1/2 ) A(Xk+1 ) = k+.
(0) (0) Em b3 x g(tk+1, x0 ) b2 x g(tk+1/2, xk+1/2 ) k+ Итак, построенные итерации (14) позволяют вычислять матрицу A и ее LU -разложение только один раз в каждой точке сетки, что являет ся наиболее эффективным по затратам машинного времени. При этом, использование упрощенных итераций Ньютона вида (14) не приводит к ухудшению теоретических оценок достаточного числа итераций в узле сетки. Это вытекает из следующего результата.
Теорема 4 Пусть правая часть задачи (1) достаточное число раз диф ференцируема в окрестности решения x(t) на отрезке [0, T ]. Тогда при диаметре сетки 0 приближенное решение этой задачи xk (N ), k = 1, 2,..., K, полученное с помощью комбинированного алгоритма (14), построенного на основе коллокационного метода (10) с коэффици ентами (11) (12), существует и сходится к точному решению x(t).
Кроме того, для сеток с достаточно малым диаметром справедлива оценка ошибки x(tk ) xk (N ) Cµ µ, k = 1, 2,..., K, (15) где µ = min{2N, 2p + 4}, = max {k }, а Cµ константа.
0kK Формула (15) имеет важное практическое значение, поскольку позво ляет оценить минимальное число упрощенных итераций Ньютона, кото рые необходимо выполнить для достижения максимального порядка схо димости 2p + 4. Итак, из (15) получаем, что число упрощенных итераций Ньютона (14) в каждой точке сетки следует подчинить условию N p + 2. (16) Данные вычислительных экспериментов, приведенные в конце пара графа, подтверждают все полученные выше теоретические результаты.
В параграфе 4.3 процедура автоматического управления длиной шага и порядком, разработанная для одношаговых экстраполяционных мето дов Рунге-Кутты в главе 3, обобщена на Е-методы со старшими произ водными. При построении этого управления, необходимо учитывать неко торые нюансы, связанные с использованием Е-методов со старшими про изводными. Обсуждаемый локально-глобальный контроль точности ба зируется на экстраполяционной технике, а значит порядок численного решения, вычисляемого в каждой точки сетки tk+1, будет меняться в за висимости от числа q интегрирований на локальном интервале [tk, tk+1 ].
Очевидно, что минимальное число упрощенных итераций (14), которое ограничивалось снизу условием (16), должно теперь соответствовать точ ности численного решения и тоже зависеть от параметра q. Принимая во внимание, что Е-методы (10) с коэффициентами (11)–(12) являются сим метричными, а значит, могут быть использованы для построения квад ратичной экстраполяции, приходим к выводу, что порядок экстраполя ционной схемы, основанной на Е-методах со старшими производными, вычисляется по формуле 2(p + q 1) + 4. Тогда из оценки ошибки (15) можно получить условие для достаточного числа упрощенных итераций Ньютона. А затем, чтобы обеспечить максимальную точность вычисле ний, полностью исключив влияние ошибки итерационного процесса на ошибку комбинированного метода, и сделать контроль точности асимп тотически корректным, мы увеличиваем количество итераций в каждой точке сетки на единицу, т.е. в дальнейшем будем вычислять достаточное число итераций (14) по формуле N = p + q + 2. (17) Отметим, что условие (17) является ключевым для надежного функцио нирования всей процедуры локально-глобального контроля точности пе ременного порядка. Как и в главе 3, работоспособность обсуждаемой про цедуры контроля локальной и глобальной ошибок Е-методов со старшими производными зависит от трех параметров, устанавливаемых пользовате лем, а именно от loc допуска на локальную ошибку комбинированного Е-метода, glob допуска на глобальную ошибку комбинированного Е метода и lg допуска на локальную ошибку неявного метода Эйлера, примененного для приближенного решения уравнения (4) относительно коэффициентов двух главных членов разложения глобальной ошибки (3).
Для вычислительных экспериментов мы выбираем следующее соотноше ние между ними:
loc = glob = lg.
Таким образом, достаточно задать один из этих допусков, чтобы опреде лить их все.
Далее, в параграфе 4.4 представлено всестороннее тестирование раз работанной вычислительной техники на дифференциальных уравнениях с известным решением и анализ результатов экспериментов. Результаты вычислительных экспериментов подтверждают теоретические положения главы 4.
В Заключении кратко описаны основные направления диссертацион ного исследования и перечислены следующие основные результаты:
1. Построены и протестированы алгоритмы автоматического контро ля локальной и глобальной ошибок экстраполяционных методов, ис пользующих неявные одношаговые методы Рунге-Кутты и коллока ционные методы со старшими производными, для систем обыкновен ных дифференциальных уравнений.
2. Разработана методика построения присоединенных и симметрич ных коллокационных методов со старшими производными, выведены формулы для вычисления коэффициентов таких методов.
3. Разработана методика построения упрощенных итераций Ньютона, доказана сходимость комбинированных схем неявных экстраполяци онных методов с использованием нового итерационного процесса, вы ведена асимптотическая оценка погрешности.
4. Разработан комплекс программ, написанных на языке C++ и ориен тированных на решение систем обыкновенных дифференциальных уравнений с заданной точностью.
Приложение, состоящее из одного параграфа, содержит описание библиотеки программ, использованной для проведения всех практиче ских расчетов в рамках работы над диссертацией. Библиотека является набором процедур, написанных на языке C++ и ориентированных на ре шение систем обыкновенных дифференциальных уравнений с заданной точностью.
Автор выражает глубокую благодарность своему научному руководи телю доктору физико-математических наук, доценту Куликову Геннадию Юрьевичу за постоянное внимание и помощь в работе.
Список публикаций по теме диссертации.
Публикации в изданиях, входящих в перечень ВАК:
1. Хрусталева, Е.Ю. Об автоматическом управлении размером ша га и порядком в явных одношаговых экстраполяционных методах / Г.Ю. Куликов, Е.Ю. Хрусталева// Журнал вычислительной матема тики и математической физики. 2008. Т. 48, № 8. C. 1392–1405.
2. Хрусталева, Е.Ю. Об автоматическом управлении размером шага и порядком в неявных одношаговых экстраполяционных методах / Г.Ю. Куликов, Е.Ю. Хрусталева// Журнал вычислительной матема тики и математической физики. 2008. Т. 48, № 9. C. 1580–1606.
3. Хрусталева, Е.Ю. Об автоматическом управлении длиной шага и порядком в одношаговых коллокационных методах со старшими производными / Г.Ю. Куликов, Е.Ю. Хрусталева// Журнал вычис лительной математики и математической физики. 2010. Т. 50, № 6. C. 1060–1077.
Публикации в прочих изданиях:
4. Хрусталева, Е.Ю. Реализация и тестирование алгоритма управле ния порядком и шагом явных экстраполяционных методов// Фунда ментальные проблемы математики и механики. Ульяновск: УлГУ, 2002. Вып. 11. С. 146–153.
5. Хрусталева, Е.Ю. О квадратичной экстраполяции, основанной коллокационных методах со старшими производными /Г.Ю. Ку ликов, А.И. Меркулов, Е.Ю. Хрусталева// Фундаментальные про блемы математики и механики. Ульяновск: УлГУ, 2003.
Вып. 1(13). С. 48–62.
6. Хрусталева, Е.Ю. О различных способах управления порядком и шагом явных экстраполяционных методов// Труды пятой междуна родной конференции Математическое моделирование физических, экономических, технических, социальных систем и процессов.
Ульяновск: УлГУ, 2003. С. 200–202.
7. Хрусталева, Е.Ю. Симметричные коллокационные одношаговые методы со старшими производными и квадратичная экстраполя ция /Г.Ю. Куликов, А.И. Меркулов, Е.Ю. Хрусталева// Материа лы XXVI конференции молодых ученых механико-математического факультета МГУ им. М.В. Ломоносова. Москва: Механико-мате матический факультет МГУ, 2004. С. 133–134.
8. Хрусталева, Е.Ю. Симметричные коллокационные одношаговые методы со старшими производными и квадратичная экстраполяция /Г.Ю. Куликов, А.И. Меркулов, Е.Ю. Хрусталева// Труды XXVI конференции молодых ученых механико-математического факульте та МГУ им. М.В. Ломоносова. Москва: Механико-математический факультет МГУ, 2004. Т. 2. С. 145–148.
9. Khrustaleva, E.Yu. On a family of A-stable collocation methods with high derivatives /G.Yu. Kulikov, A.I. Merkulov, E.Yu. Khrustaleva// In: Bubak M. et al (eds.): Computational Science ICCS 2004.
4th International Conference, Krakow, Poland, June 6–9 2004.
Proceedings, Part II. Lecture Notes in Computer Science, 2004.
V. 3037. P. 73–80.
10. Khrustaleva, E.Yu. Symmetric Runge-Kutta methods with higher derivatives and quadratic extrapolation / G.Yu. Kulikov, E.Yu. Khrustaleva, A.I. Merkulov// In: Alexandrov V.N. et al (eds.): Computational Science ICCS 2006. 6th International Conference, Reading, UK, May 28–31, 2006. Proceedings, Part I. Lecture Notes in Computer Science, 2006. V. 3991. P. 117–123.