На правах рукописи
АЛЬШИНА Елена Александровна
МЕТОДЫ ВЫЧИСЛЕНИЙ С КОНТРОЛЕМ ТОЧНОСТИ
НА КВАЗИРАВНОМЕРНЫХ СЕТКАХ
01.01.07 – вычислительная математика
Автореферат
диссертации на соискание ученой степени
доктора физико-математических наук
МОСКВА
2006
Работа выполнена в Институте математического моделирования Российской Академии Наук
Научный консультант:
|
член-корр. РАН
Николай Николаевич Калиткин
|
Официальные оппоненты:
|
член-корр. РАН,
доктор физико-математических наук, профессор
Сергей Тимофеевич Суржиков
|
|
доктор физико-математических наук, профессор
Александр Николаевич Боголюбов
|
|
доктор физико-математических наук, профессор
Александр Петрович Михайлов
|
Ведущая организация:
|
Институт прикладной математики
им. М.В.Келдыша РАН
|
Защита состоится « 21 » ноября 2006 г. в ___ ч. ____ мин. на заседании Диссертационного совета Д 002.058.01 в Институте математического моделирования РАН по адресу 125047, г. Москва, Миусская пл. 4а
С диссертацией можно ознакомиться в библиотеке ИММ РАН.
Автореферат разослан « 16 » октября 2006 г.
Ученый секретарь
диссертационного совета,
д.ф.-м.н.
|
|
Н.В. Змитренко
|
ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ
Актуальность темы. Проверяя адекватность математических моделей путем сравнения с экспериментом нужно быть уверенным, что расчеты проведены с точностью, по крайней мере, на порядок лучше, чем погрешность измеряемых входных данных. Для практического применения результатов расчетов крайне важно не только получить численное решение моделируемой задачи, но и знать, с какой математической точностью этот результат получен. Со времен работ Лагранжа и особенно Коши всякий, установленный числено, результат принято сопровождать надежной оценкой погрешности.
Если говорить, например, о численном решении ОДУ, то наиболее распространенными во всем мире сейчас являются программы расчетов с автоматикой выбора шага. Наши многочисленные тесты таких программ на задачах с известным точным решением показывают, что реальная точность расчета может отличаться от декларируемой в 10, 100 и даже несколько тысяч раз, а при расчетах на длительные времена численное решение вообще может не иметь ничего общего с точным. Сложившаяся ситуация не может не вызывать опасений.
Тем ни менее надежный способ получения апостериорной асимптотически точной оценки погрешности путем расчета на вложенных сгущающихся сетках предложен еще Ричардсоном уже почти 100 лет назад. В логически законченном виде этот подход изложен в работе Ричардсона [L.F. Richardson, 1927] и известен в учебной литературе как метод Ричардсона, экстраполяция по Ричардсону, правило Рунге. Монография [Г.И. Марчук, В.В. Шайдуров, 1979] посвящена практическим аспектам применения этого метода. Основополагающей является книга [Рябенький В.С., Филиппов А.Ф., 1956], где дано общее достаточное условие применимости метода Ричардсона к сеточным задачам.
Как писал сам Ричардсон «развиваемая теория сложна и наводит на мысль, что практика может быть столь же сложной, тогда как в действительности она проста». Возможно, упомянутая Ричардсоном сложность препятствует широкому применению этого метода. Но метод сгущения сеток, являющийся мощным средством для численного решения широкого круга задач с гарантированной точностью, явно недооценен исследователями. Даже математики прикладники не всегда знают возможности этого метода и как им грамотно пользоваться. Большинство же инженеров, рассчитывающих свои задачи на компьютерах, почти не слышали о нём.
Возможность многократного рекуррентного уточнения по методу Ричардсона дает впечатляющий выигрыш в точности и экономичности. В западной литературе этому вопросу посвящены работы Грэга [W.B. Gragg, 1965], Штеттера [B. Gragg, H.J. Stetter, 1964], [H.J. Stetter, 1970], Бауэра, Рутисхаузера, Штифеля [F.L. Bauer, H. Rutishauserand E. Stiefel, 1963] и многих других авторов. Частным случаем этого подхода является метод рациональной экстраполяции Булирша и Штера [R. Bulirsh and J. Stoer 1964, 1966]. Обоснование экстраполяционных методов можно найти в монографии [Штеттер Х., 1978], и в книге [Э. Хайрер и др, 1990].
Другая причина недостаточно активного использования метода сгущающихся сеток в том, что классическая оценка погрешности по Ричардсону была написана только для равномерных сеток. Равномерные сетки невыгодны для мало-мальски сложных задач. Функция и ее производные могут сильно меняться на рассматриваемом промежутке. На тех участках, где производные функции малы вполне допустимо выбирать крупные шаги сетки; там, где производные велики, целесообразно выбирать более густую сетку. Сетки, приспособленные к поведению функции, называют адаптивными. Кроме того, в расчетах используют неструктурированные сетки: с треугольными и пирамидальными ячейками. Этому направлению посвящено огромное число работ, например Азаренка Б.Н., Вабищевича П.Н., Забродина А.В., С.А. Иваненко, Лисейкина В.Д., Петренко В.Е., Г.П. Прокопова, А.Ф. Сидорова, О.В. Ушаковой, Тишкина В.Ф. (а также многих других авторов) и многочисленные конференции. Однако на таких сетках не предложено способов построения аппроксимации высокого порядка точности, не удается воспользоваться классическим методом Ричардсона, получать асимптотически точную оценку точности и проводить рекуррентное уточнение.
Среди неравномерных сеток есть один важный класс квазиравномерные. Квазиравномерные сетки предложил А.А.Самарский около 1952 г. (в закрытых отчетах). Опубликована эта идея была в 1960-е г. А.Ф. Сидоровым. Строгое определение таких сеток и возможность обобщения на бесконечную область были даны научным консультантом диссертационной данной работы Н.Н.Калиткиным и опубликованы в 1978 г.
В данной работе показано, что такие сетки легко адаптируются ко многим задачам, и при этом позволяют использовать все возможности метода сгущения сеток, получать гарантированную асимптотически точную оценку для погрешности расчета, проводить рекуррентное уточнение. Это позволяет добиваться высокой точности при умеренном объёме вычислений даже для весьма сложных задач.
В практике нередко возникают задачи с сингулярными решениями, однако, имеется ли эта сингулярности, и, в какой именно точке, – заранее не всегда известно. Примерами являются задачи на разрушение материалов, электрический пробой (в том числе возникновение молнии) и другие. Один из разделов этой работы посвящен построенному в начале 2005 г. в развитие идей расчетов с контролем точности методу диагностики сингулярности, позволяющему чисто расчетно определить заранее неизвестные местоположение и тип особенности точного решения.
Важным классом практических задач являются жесткие и дифференциально-алгебраические системы.
Разработкой численных методов решении таких задач в разные годы занимались Абрамов А.А., Артемьев С.С., Арушанян О.Б., Бобков В.В., Боднарчук П.И., Герасимов Б.П., В.Н. Гридин, Демидов Г.В, Заворин А.Н., Захаров А.Ю., Калиткин Н.Н., Кузнецов Е.Б., Лебедев В.И., Е.А. Новиков, Новиков В.А., Макаров В.Л., Мейгауз М.Г., В.Б. Михайлов, Павлов Б.В., Ракитский Ю.В., Устинов С.М., Федоренко Р.П., Филиппов С.С., Черноруцкий И.Г., Шалашилин В.И., Ширков П.Д., Raimondas Čiegis, Bill Gear, Ernst Hairer, Alan Hindmarsh, Olavi Nevanlinna, Kees Dekker, Jan Verwer, Gerhard Wanner, Michel Roche.
Жесткие системы возникают в задачах химической кинетики, радиотехники, электроники, при решении задач теплопроводности и диффузии методом прямых и многих других. Они традиционно трудны для численного решения.
Трудности решения сверхжестких и дифференциально-алгебраических систем схожи между собой. Обзор лучших существующих подходов дан, например, в монографии [Хайрер Э., Ваннер Г., 1999]. Методы, описанные в этой монографии, реализованы в стандартных программах, они доступны, например, в сети Internet. Но наши тесты таких программ неоднократно показывали, что их точность может сильно отличаться от заявленной даже на задачах умеренной жесткости.
Кроме того, в литературе, как правило, описаны только схемы с действительными коэффициентами, тогда как использование комплексных коэффициентов повышает число степеней свободы и позволяет строить гораздо лучшие по точности и устойчивости схемы с небольшим числом стадий. Применению схемы Розенброка с комплексным коэффициентом для решения сверхжестких и дифференциально-алгебраических задач с гарантированной точностью посвящена одна из глав диссертации.
Цель работы состоит в построении численных методов для широкого класса задач математической физики и их реализации в алгоритмах расчетов с гарантированной точностью. Особое внимание уделено надежности предложенных методов даже для таких сложных вычислительных задач, как жесткие и дифференциально-алгебраические системы или вычисление несобственных интегралов от функций со слабым степенным убыванием.
Научная новизна. Алгоритм расчета с контролем точности методом вложенных сгущающихся сеток, с использованием рекуррентного повышения точности обобщен и обоснован для случая квазиравномерных сеток, в том числе в неограниченных областях. Разработан ряд новых алгоритмов, пригодных для неограниченных областей. Показаны преимущества применения схемы Розенброка с комплексным коэффициентом для интегрирования сверхжестких и дифференциально-алгебраических систем с гарантированной асимптотически точной оценкой погрешности. Построен метод чисто расчетной диагностики местоположения и типа особенностей заранее неизвестного точного решения.
Достоверность. В важнейших случаях строго доказаны теоремы о сходимости предложенных методов. Надежность метода подтверждена большим числом тестовых расчетов на задачах с известным точным решением.
Практическая значимость. Построены и обоснованы эффективные методы решения ряда задач, которые ранее не поддавались усилиям математиков-прикладников или решались с большим трудом. В частности, показано, что ряд известных алгоритмов и стандартных программ на самом деле не обеспечивает той точности, которую они обещают пользователю. В то время как развитые в диссертации методы действительно позволяют проводить расчеты с гарантированной точностью. Написаны программы решения ряда прикладных задач, которые могут стать основой стандартных программ.
Апробация работы. Основные результаты работы докладывались на многих Всероссийских и зарубежных конференциях. В том числе, на Всемирных Конгрессах Математиков (Пекин, 2002, Мадрид, 2006), на Европейском Конгрессе Математиков (Стокгольм, 2004), на Международных конференциях EquaDiff 2003, 13-th European Conference on Mathematics for Industry (Эйндховен, 2004), 10th International Conference Mathematical Modelling and Analysis и 2nd International Conference Computational Methods in Applied Mathematics (Тракай, 2005), ICTMA 12 Teaching of Mathematical Modelling and Applications, London, 2005, на Международной конференция «Ломоносовские чтения» (МГУ, 2004), на II Всероссийской конференции "Актуальные проблемы прикладной математики и механики", посвященной памяти академика А.Ф.Сидорова (2004), на Международной конференций студентов, аспирантов и молодых учёных по фундаментальным наукам «Ломоносов-2005» (Москва, МГУ), IX, X, XI Всероссийских школах-семинарах “Современные проблемы математического моделирования” (Абрау-Дюрсо, 2001, 2003, 2005) и др.
В 2004 г. был сделан доклад на Семинаре Научно-технического совета по фундаментальным исследованиям ГНЦ РФ ФЭИ им. А.И. Лейпунского (г. Обнинск), в 2005 г. результаты докладывались на семинаре по нелинейным уравнениям под руководством члена-корр. РАН И.А. Шишмарева (ВМиК МГУ). Результаты, составившие суть диссертации, докладывались в марте 2005 г. на семинаре математического сектора Федерального Ядерного Центра (г. Саров) и на юбилейной конференции, посвященной 85-летию А.А. Самарского. По результатам диссертации по приглашению Лондонского Математического Общества был прочтен цикл лекций весной 2006 г. в Университетах Кембриджа и Манчестера.
Публикации. Основные результаты опубликованы в 1 монографии (издательство ФИЗМАТЛИТ), 23 журнальных статьях (включая 4 в «Докладах Академии Наук») и 27 полнотекстовых публикациях докладов на зарубежных и всероссийских конференциях.
Структура и объем работы. Диссертация состоит из введения и 5 глав, списка основных результатов и перечня цитируемой литературы. Диссертация содержит 208 страниц, в общей сложности 77 рисунков и 19 таблиц. Список цитируемой литературы содержит 127 наименований.
СОДЕРЖАНИЕ РАБОТЫ
Введение содержит обзор существующих результатов по теме диссертации. Показана актуальность и практическая значимость работы, сформулированы цели диссертации.
Глава I состоит из 5 параграфов. В ней рассмотрены способы апостериорной асимптотически точной оценки погрешности. Показано, как можно повысить точность, если выполнить расчёты на нескольких равномерных сетках с разным числом узлов. Изложены практические приёмы контроля точности и диагностики ошибок, полезные для составления и отладки программ.
В первом параграфе с использованием формулы Тейлора и фундаментальных теорем В.С. Рябенького и А.Ф. Филиппова описана типичная структура погрешности сеточных методов в виде разложения в сумму обратных степеней числа узлов сетки
где погрешность сеточного метода, число узлов сетки, порядок точности метода, константы, зависящие от конкретного решения и его производных, но не зависящие от .
Определена связь между гладкостью решения и числом слагаемых в разложении . На простейшем примере квадратурных формул выписаны главные члены в разложении .
Второй параграф посвящен сгущению равномерных сеток. Дана формула Ричардсона
,
позволяющая по результатам расчетов , на двух соседних сетках с числами узлов N и rN выразить главный член в разложении и тем самым получить асимптотически точную оценку для погрешности сеточного расчета.
Выражение можно учесть в качестве поправки, исключив главный член в разложении погрешности , и повысить тем самым точность расчета
.
При наличии расчетов на большем числе сеток, процесс уточнения можно сделать рекуррентным. Число уточнений, как и число слагаемых в сумме ограниченно гладкостью решения.
Также обсуждается применение метода Ричардсона к многомерным задачам и даны оптимальные наборы сеток.
Третий параграф посвящен практическим аспектам написания и отладки работы программ расчетов с контролем точности. Описан алгоритм тестирования программ на задачах с известным точным решением. Можно определить эффективный порядок точности метода по углу наклона графика убывания погрешности с ростом числа узлов в двойном логарифмическом масштабе (рис. 1).
Формула для определения эффективного порядка точности -ого уточнения на -ой сетке следующая
.
Здесь соответствует расчету по базовому сеточному алгоритму (без уточнения).
Если эффективный порядок ниже теоретического, тогда как известно, что решение обладает необходимой гладкостью, это свидетельствует об ошибке в программе.
Если при расчете по отлаженной программе эффективный порядок точности близок к теоретическому, это свидетельствует о том, что все остальные слагаемые в разложении пренебрежимо малы по сравнению с главным и можно провести уточнение . Если эффективный порядок точности очередного уточнения заметно ниже теоретического значения, то гладкость решения недостаточна для дальнейших уточнений по Ричардсону. Повышение точности при этом возможно лишь за счет дальнейшего сгущения сетки.
Если при очередном сгущении сетки убывание погрешности не соответствует теоретическому порядку точности, тогда как сама величина погрешности мала, то это свидетельствует о выходе на ошибки округления. Стрелочки на рис. 1 указывают на заметное влияние ошибок округления.
|
|
Рис.1 а,б. Погрешности (а) и эффективные порядки точности (б) при рекуррентном сгущении.
|
Эффективный алгоритм расчетов с контролем точности состоит в том, что нужно провести серию расчетов, каждый раз сгущая число узлов сетки в одно и то же число раз, оценивая погрешность по формуле Ричардсона, реализовать возможное число уточнений. Если модуль очередной поправки не превосходит заданного уровня точности, то расчет можно прекратить, выдавая пользователю результат последнего уточнения в качестве ответа, а модуль последней поправки в качестве асимптотически точной оценки погрешности результата.
Включение в программу описанного алгоритма позволяет достигнуть заданной высокой точности при весьма скромных вычислительных затратах. Алгоритм не имеет насыщения в том смысле, что автоматически будет использовано столько непрерывных производных решения, сколько у него есть.
В четвертом параграфе идеи расчетов с контролем точности обобщаются на случай, когда гладкость решения недостаточна для реализации теоретического порядка точности выбранного численного метода. В большинстве таких случаев можно считать, что погрешность сеточного метода разлагается в сумму по дробным степеням числа узлов сетки [Марчук Г.И., Шайдуров В.В.]
,
Степени , как правило, заранее неизвестны. Проведем серию расчетов на сгущающихся сетка , каждый раз увеличивая число узлов в одно и то же число раз . Тогда по результатам расчетов на трех соседних сетках можно выразить главный член погрешности в разложении по формуле Эйткена:
, .
В работе показано, что формула дает асимптотически точную оценку погрешности расчета, что позволяет сделать процесс уточнения по Эйткену рекуррентным. Приведен пример применения алгоритма рекуррентного уточнения по Эйткену к вычислению несобственного интеграла . Стандартные пакеты математических программ (за исключением тех, которые «умеют» брать этот интеграл аналитически) не позволяют получить на этой задаче точности выше, чем 10%. Расчет по данному алгоритму позволил достичь точности при весьма скромном числе узлов сетки .
Если идеи рекуррентного сгущения по Ричардсону были известны ранее, то обоснование метода рекуррентного сгущения по Эйткену является новым, выносимым на защиту результатом.
В пятом параграфе описан контроль точности в стандартных пакетах математических программ. При численном решении ОДУ наиболее распространенными во всем мире сейчас являются программы расчетов с автоматикой выбора шага. Оценка погрешности в них производится априорно, она не является ни мажорантой, ни асимптотически точной. Она делается по «кухонным алгоритмам», отлаженным на большом числе тестовых примеров. Как правило, стандартная программа поставляется пользователю в виде «черного ящика», в котором можно задать требуемую точность (так называемый параметр tolerance) и программа выполняет расчет «якобы с этой точностью». Наши многочисленные тесты таких программ на задачах с известным точным решением показывают, что реальная точность расчета может сильно отличаться от декларируемой.
В качестве примера приведен расчет орбиты Аренсторфа по рекомендованной в [Хайрер Э. и др., 1990] программе DOPRI5, где реальная точность по скоростной компоненте решения оказалась в раз хуже заявленной.
Возможной причиной описанного сбоя стандартных программ является неинтерполяционность используемых в них схем. В данном параграфе исследован выбор оптимальных параметров явных схем Рунге-Кутты невысокого порядка точности, минимизирующий погрешность и обеспечивающий интерполяционность. На конкретном примере проиллюстрирована опасность применения неинтерполяционных схем.
|
