автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Жесткие и плохо обусловленные нелинейные модели и методы их расчета

На правах рукописи

Пошивайло Илья Давлович

Жесткие и плохо обусловленные нелинейные модели и методы их расчета.

Специальность 05.13.18 — Математическое моделирование, численные методы и комплексы программ

Автореферат диссертации на соискание учёной степени кандидата физико-математических наук

2 9 ДПР 2015

005567947

Москва — 2014

005567947

Работа выполнена в федеральном государственном бюджетном учреждении науки Институт прикладной математики им. М.В. Келдыша Российской академии наук.

Научный руководитель: Калиткнн Николай Николаевич,

член-корреспондент РАН, доктор физико-математических наук, профессор.

Официальные оппоненты: Гончаров Виктор Анатольевич,

доктор физико-математических наук, доцент, профессор кафедры «Высшая математика №1», Федеральное государственное автономное образовательное учреждение высшего образования «Национальный исследовательский университет «Московский институт электронной техники». Скворцов Леонид Маркович, кандидат технических наук, зав. лабораторией САПР факультета «СМ», Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования «Московский государственный технический университет имени Н.Э. Баумана».

Ведущая организация: Кафедра математики физического факультета

Московского государственного университета им. М.В. Ломоносова

Защита состоится «21» мая 2015 г. в Ц часов на заседании диссертационного совета Д 002.024.03 на базе института прикладной математики им. М.В. Келдыша РАН по адресу: 125047, Москва, Миусская пл., д.4.

С диссертацией можно ознакомиться в библиотеке института прикладной математики им. М.В. Келдыша РАН.

Автореферат разослан «20» апреля 2015 г.

Ученый секретарь диссертационного совета Д 002.024.03, д.ф.-м.н.

Змитренко Н.В.

Общая характеристика работы

Актуальность темы. Настоящая диссертационная работа посвящена численному решению моделей, возникающих в задачах физики плазмы, химической кинетики и ряде других прикладных областей. Такие модели часто описываются жесткими или плохо обусловленными системами обыкновенных дифференциальных уравнений.Моделированию и расчету жестких систем посвящено множество работ. Основные методы и подходы были предложены в работах Далквиста, Гиршфельдера, Кертисса, Ваннера, Розенброка и их последователей. В настоящее время в данной области ведутся исследования по нескольким направлениям. Во-первых, исследуются семейства явно-неявных схем Розенброка. Эти схемы обладают высокой устойчивостью и являются экономичными, так как не требуют итераций. Во-вторых, исследуются возможности применения явных схем с автоматикой выбора шага интегрирования для решения жестких задач определенного типа. В-третьих, строятся эффективные реализации диагонально-неявных методов Рунге-Кутты, позволяющие экономить вычисления за счет алгоритмов автоматики выбора шага, использования модифицированного метода Ньютона для решения нелинейных систем и нетривиального предиктора для начального приближения на новом временном слое. Обилие направлений исследований говорит о востребованности эффективных и надежных методов решения жестких и плохо обусловленных систем ОДУ.

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

1. Разработать полностью неявные Ьр-устойчивые схемы Рунге-Кутты порядка точности р от 1 до 4, допускающие экономичную реализацию по сравнению с классическими полностью неявными схемами Рунге-Кутты.

2. Провести исследование возможности автономизации исходной задачи методом перехода к длине дуги. Для преобразованной новой задачи найти оценку точности получаемого решения по методу Ричардсона.

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

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

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

Научная новизна:

1. В рамках данной работы выведен новый подкласс разностных схем из семейства полностью неявных схем Рунге-Кутты порядка точности р от 1 до 4. Эти схемы обладают Ьр-устойчивостью, а по трудоемкости вычислений эквивалентны лишь одной стадии диагонально-неявных схем Рунге-Кутты, обладающих лишь 1,1-устойчивостью.

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

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

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

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

Методы исследования. При работе над диссертацией использовались теория дифференциальных уравнений, разностных схем, численных методов, математического анализа. Численное моделирование проводилось с использованием программ на языках С++ и МАТЬАВ/Ос1ауе.

Основные положения, выносимые на защиту:

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

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

3. Проведены расчеты моделей, описывающих прикладные задачи из плазменной теплопроводности, химической кинетики и нелинейной радиотехники.

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

Личный вклад. Автор под руководством H.H. Калиткина построил оптимальные обратные схемы Рунге-Куггы, разработал численный алгоритм их реализации, получил коэффициенты в форме матрицы Бутчера, выполнил адаптацию к дифференциально-алгебраическим задачам и трансформировал эти схемы для выбора длины дуги в качестве аргумента.

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

Апробация работы. Результаты диссертации докладывались на конференции "Современные проблемы вычислительной математики и математической физики" памяти академика A.A. Самарского к 95-летию со дня рождения (Москва, 16-17 июня 2014). Также по материалам диссертации были сделаны доклады на семинарах кафедры вычислительной математики ВМК МГУ (декабрь 2013), Института прикладной математики им. М.В. Келдыша РАН и кафедры математики физического факультета МГУ (апрель 2014).

Публикации. По теме диссертации всего опубликовано 7 работ в журналах, входящих в список ВАК: «Доклады Российской академии наук» - 2, «Журнал вычислительной математики и математической физики» - 1, «Математическое моделирование» - 4.

Объем и структура работы. Диссертация состоит из введения, пяти глав и заключения. Полный объем диссертации 89 страниц текста с 23 рисунками и 7 таблицами. Список литературы содержит 62 наименования.

Содержание работы

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

$ = f(u,t), u,f€KM, u(0) = uo, 0 <t<T. (1)

at

Явные схемы Рунге-Кутты, Адамса и другие непригодны для решения систем, полученных в результате моделирования исследуемых задач. В мировой литературе принято называть такие задачи жесткими. Далее в работе дается классификация трудностей решения систем (1) в зависимости от спектра матрицы Якоби правой части fu, формулируются требования к численным методам, пригодным для решения жестких систем, и необходимость получения надежной оценки погрешности решения. Далее сформулированы цель работы, научная новизна, достигнутые результаты и их значимость. Приводится обзор научной литературы по теме работы.

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

Обозначим решение на исходном шаге через и, а на новом через й (и, й 6 ). В общем виде формулы схем Рунге-Кутты записываются следующим образом:

я L

й = u + wfc ~ f(u + r^afc|W/,i + тс*). (2)

к=1 1=1

Матрица коэффициентов {а^} называется матрицей Бутчера. Векторы коэффициентов {сЛ.} и {Ьд,.} можно рассматривать как столбцы, дополняющие квадратную матрицу Бутчера до прямоугольной.

Обычно выделяют следующие классы схем. Если L = к — 1, т.е. лишь находящиеся ниже главной диагонали элементы матрицы Бутчера отличны от нуля, схемы называют явными (ERK - explicit Runge-Kutta). В них переход на

новый слой происходит по явным формулам. Явные схемы просты в реализации,

б •

а трудоемкость вычислений мала. Если Ь = к,то схемы называют диагонально-неявными (DIRK - diagonal implicit Runge-Kutta). В таких схемах на каждой стадии для нахождения очередного wfc необходимо решить систему алгебраических уравнений порядка М. Если же L = s, схемы называют полностью неявными схемами PK (FIRK - fully implicit Runge-Kutta). В них для нахождения всех wfc требуется решить систему алгебраических уравнений порядка sM.

Явные схемы ERK с s < 4 стадиями хорошо изучены. Часть коэффициентов схем выбирается для выполнения условий порядка, т.е. чтобы порядок точности схемы р = s. Наложим дополнительные требования: интерполяционное™ схемы и минимального числа слагаемых в невязке схемы. Тогда для каждого s выделяется единственная оптимальная схема, у которой в матрице коэффициентов Бутчера ненулевыми оказываются только элементы нижней кодиагонали afcft_i, где 2 < к < s. Из условий аппроксимации следует, что ак.к-1 = ск, где ск - дополнительный столбец матрицы Бутчера (коэффициенты при г). Тогда, обозначив aktk-1 = ск = ак, оптимальные явные схемы PK можно записать следующим образом:

s

û = u + тЪкwfc, s <4; wfc = f(u + такwt_i, t + так).

(3)

fc=i

Обратные схемы PK построим следующим образом. Рассмотрим движение в обратном направлении. Напишем явную схему для перехода t t. Для этого в (3) нужно поменять местами t f) i, u н û и изменить знак т. Получим следующие оптимальные обратные схемы ("backward optimal Runge-Kutta",

BORK):

S

û= u + r^bfcWfc, Wfc = f (û — TdfcWfc-i, î — rak). (4)

k=1

Коэффициенты этих схем те же, что и для схем (3). Очевидно, схемы (4) также

имеют точность 0(rs), s < 4.

Получим рекурсивную запись обратных схем (4). Для простоты записи

ограничимся автономными системами: f(u,i) = f(u). Тогда сдвиг w^ равен

f(û). w2 выражается уже как рекурсивная функция: w2 = f(û - ra2f(u)). То

же самое справедливо для w3 и w4. Тогда обратная схема 4-го порядка примет

следующий рекурсивный вид:

û = u + rF(û), F(û) = bif(û) + 62f(û - ra2f(û))+ (5)

b3f(û - Ta3f(û - ra2f(û))) + biî(û - ra4f(û - ra3f(û - ro2f(û)))).

Схема (5) представляет собой нелинейную алгебраическую систему порядка M для определения Û, записанную через рекурсивные функции. Трудоемкость ее решения соответствует одной стадии диагонально-неявных схем. Тем

самым построенная схема действительно имеет невысокую по сравнению с другими неявными схемами трудоемкость.

Среди неоптимальных обратных схем можно рекомендовать одну схему второго порядка точности (s = р = 2), которая является обратной к модифицированному явному методу Эйлера (в англоязычной литературе его обычно называют "explicit midpoint method"). Этот метод также входит в семейство явных двухстадийных схем Рунге-Кутты (2). Формула для явного метода выглядит следующим образом:

u = u + rf(u+irf(u,0,i + ^)). (6)

Соответствующая обратная схема ("backward midpoint method", BMP) получается аналогично оптимальным обратным схемам:

u = u + rf(u-irf(u,i-|-r),i+ir)). (7)

Оптимальные обратные схемы (4) и схема (7) обладают Ls-устойчивостью.

Алгоритм. При решении задачи (1) обратными схемами (4) или (7) для перехода на новый слой необходимо решить систему нелинейных уравнений:

й - tF(u) - и = 0, (8)

где F(fi) записано в (5). Для решения нелинейной системы обычно применяют модифицированный метод Ньютона. Классический ньютоновский итерационный процесс выглядит следующим образом:

[Е — rFuiu"7)] (й,+1 — и4) = — [й® — rF(u'J) — и], й° = и. (9)

Каждая итерация сводится к линейной системе относительно й?+1. Метод Ньютона очень чувствителен к выбору начального приближения: вблизи корня сходимость метода квадратичная, однако если начальное приближение выбрано неудачно, сходимости может не быть вовсе. Преставляется целесообразным применить следующий простой подход. Рассмотрим обобщение классического ньютоновского итерационнного процесса:

D(ufc)Ajt = -F(ujfc), Ufc+1 = u, + 9k Дь О<0*<1. (10)

При вь = 1 оно переходит в классический метод Ньютона. Выбирать значение параметра на каждой итерации нужно таким образом, чтобы

|P(ufc+i)|<|F(ufc)|, (П)

т.е. модуль правой части |F(ujt) j монотонно уменьшался от итерации к итерации. Сначала нужно положить дк = 1. Если соотношение (11) не выполняется, то нужно уменьшить значение 9к вдвое, снова произвести проверку соотношения и т.д.

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

при производных стремятся к нулю.

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

G— = f(u,t), G = const. (12)

dt

Здесь G есть постоянная матрица, которая в общем случае является особенной: О < rankG < М. Величина rankG есть число дифференциальных компонент. Если rank G = М, система является чисто дифференциальной, а при rank G = О

- чисто алгебраической.

Применим к задаче (12) обратные схемы РК (4) с s < 4:

5

й = u + T^i>fcW,;, Gwfc = f(u-rafcWfc_i,i-таО- (13)

t=i

Обратные схемы Рунге-Кутты (4) с s < 4 и схема (7) являются жестко точными: у их матрицы Бутчера первую строку составляют коэффициенты bk. Значит эти схемы обеспечивает порядок точности р = s на дифференциально-алгебраических системах индекса 1.

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

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

Пусть исходная задача (1) решается на отрезке 0 < t < Т разностной схемой порядка аппроксимации р. Введем на [0;Г] равномерную или квазиравномерную сетку nN с числом узлов N. Выполним на ней расчет и найдем сеточное решение u(t; N), где t € П„. Затем сгустим сетку вдвое и на ней также найдем численное решение u(t-2N). По двум полученным расчетам можно находить асимптотически точную оценку погрешности на сетке n2N

A(t;2N) = [u(t;2iV) — u(t; N)]/(2P — 1). (14)

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

||A2W||c = maxn(|A(xn;2iV)|),

1/2 (15)

||Д2*||/а = - —

Затем строится график зависимости lg || ДдгЦ от lg N. Если зависимость погрешности от шага сетки (числа узлов) носит степенной характер, то этот график является прямой линией. Поэтому если наблюдаемый график с хорошей точностью можно считать прямой линией с наклоном tg а = р, то фактический расчет сходится с порядком точности р, и применение метода Ричардсона правомерно.

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

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

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

Для того, чтобы описать пограничный слой, продолжим сгущать сетки, воспользовавшись при необходимости повышенной разрядностью вычислений. С некоторого сгущения узлы сетки начнут попадать в переходную область. В этот момент на графике log ||Д||(Л^) погрешность может начать возрастать. Это возрастание связано с повышенной локальной погрешностью в переходной зоне и пограничном слое.

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

10

10s 10' to"

Рис. 1: Погрешность при решении жесткой задачи ван дер Пола (32), вычисленная по методу Ричардсона на серии сгущающихся сеток.

Метод перехода к длине дуги. Метод введения нового аргумента интегрирования - длины дуги интегральной кривой в М + 1-мерном пространстве {u0 = t, ui,..., им} - впервые был предложен в 1972 году в работе E.Riks. С 1993 года метод'детально развили В.И.Шалашилин и Е.Б.Кузнецов. В частности, была доказана теорема о том, что введение длины дуги обеспечивает наилучшую

обусловленность задачи Коши.

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

Длина дуги определяется соотношением

м

(do2 = X) du° -dt (16)

m=0

Принимая l за новый аргумент, получим вместо (1) следующую систему:

du, ~dl

dUm Гш{ u) - /т 0 <т<М, /о si. (17)

Правые части (17) ^т(и) не зависят от нового аргумента I, поэтому система (17)

является автономной. Выполняется соотношение £ Р,2П = 1. Поэтому правые

т=0

части Я» невелики, и система (17) является мягкой даже в том случае, если

система (1) была плохо обусловленной или жесткой.

11

Систему (1) надо интегрировать до момента Т. Но до какого значения I = L нужно интегрировать систему (17) - не известно. Поскольку Fq > 0, то Uq = t есть монотонно возрастающая функция I. Следовательно, надо интегрировать (17) до тех пор, пока не выполнится щ(1) > Т.

Сохранение балансов в задачах химической кинетики. Модели химической кинетики образуют важный подкласс задач, сводящихся к жестким системам обыкновенных дифференциальных уравнений. Часто стандартных требований к численным методам решения жестких ОДУ, перечисленным во введении к данной работе, оказывается недостаточно для получения качественно верного результата. Например, одно из требований, специфичных для задач химической кинетики, связано с законом сохранения частиц. В химических реакциях образуются и разлагаются молекулы, но суммарное число атомов каждого химического элемента остается при этом неизменным. Для хорошего качественного поведения численного решения необходимо, чтобы в численном расчете также сохранялось число атомов "каждого сорта".

Каждое уравнение баланса можно записать следующим образом. Пусть каждая молекула Uj содержит aj атомов определенного элемента. Тогда полное

число атомов в системе 52 ajuj = const не зависит от времени. Легко видеть, j

что при этом 52 ajfj — 0: сколько атомов выходит из одних молекул, столько же j

попадает в другие. Эти соотношения удобно записать в матричном виде:

ccTf = 0, aTu = const, (18)

Различных столбцов а существует столько, сколько различных химических элементов входит в реакции; исходная система должна удовлетворять всем этим балансам.

В данной работе доказаны теоремы о том, что схемы Рунге-Кутты (явные и неявные) и одностадийные схемы Розенброка сохраняют химический баланс системы. Кроме того, показано, что химический баланс системы сохраняется при переходе к длине дуги.

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

Классический метод Ньютона. В случае единственного уравнения с одним неизвестным классический метод Ньютона имеет следующий вид:

г5+1 = г, - Пх,)//'(х,). (19)

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

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

Теорема 1. Пусть заданы одношаговый итерационный процесс а;5+1 = Ф(ха,1{х3), }'{х3)), конечное число итераций и нулевое приближение ж0- Тогда можно найти такой многочлен /(х), что за Б итераций указанный процесс не сойдется к корню из выбранного нулевого приближения

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

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

Определение кратности корня. Пусть у функции одного переменного f(x) существует непрерывная а х, есть р-кратный корень. Тогда в

малой окрестности корня

/(®) » аДр + ЬАр+\ А = х-Х4. (20)

Подставляя (20) в формулу для метода Ньютона (19) и вычитая из обеих частей равенства х„ получаем:

= —Д4. + 0(Д?). (21)

р

Отсюда видно, что метод Ньютона для простого корня (р = 1) сходится квадратично, а для кратного корня (р > 2) - линейно.

По скорости сходимости можно определить кратность корня. Из (21) следует, что

/Ч'"+0(Д5). (22)

р - 1 Д3+1

Р А*

Величины Д4+ь Д., нам неизвестны, поскольку в них входит неизвестное х,. Однако нетрудно показать, что таким же, с точностью до малых величин, является отношение (жа+1 - х8)/{хв - 2:^-1). Удобно ввести знаменатель линейной сходимости д., и аппроксимацию кратности корня р.,:

Яз = -~-, Р* = --• (23)

ха - Хц-\ р 1 - д4

При сделанных предположениях р5 —> р при в —> оо. Поэтому нужно ввести в ньютоновские итерации расчет р, по (23). Если итерации сходятся, а р., стремится к некоторому числу р, то мы находим значени корня х* и его кратность.

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

Для ускорения сходимости ньютоновских итераций, когда кратность корня заранее неизвестна, целесообразно применить следующий прием. Если значение знаменателя линейной сходимости д4 (23) установилось, то по трем значениям ха, можно просуммировать дальнейшую геометрическую прогрессию и найти экстраполяцию:

х, и х.+х + (г,+1 - ха). (24)

1 — ?а

Это усовершенствование не дает выигрыша на простых корнях, но для кратных корней приводит не только к существенному ускорению сходимости, но и

позволяет достичь более высокой точности.

14

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

Бегущая тепловая волна. Рассматривается квазилинейное уравнение теплопроводности на полупрямой 0 < х < +оо с заданными граничными и начальным условиями:

д и ди\

= (25)

и(0,г) = (шА/х0)1/т, и(+оо^) = 0, и(х, 0) = 0. Оно имеет частное непрерывное решение вида волны, бегущей со скоростью с:

1/т

' * £ <*' (26) х > сЬ.

cm, . . . , —tct — x) u(x,t) = П ^о

Вводя равномерную сетку {х„} с шагом h и обозначая ип — и(хп), получим для (25) схему метода прямых точности 0(h2):

-¿jr = /»(u) = ■^[нп+1/2(ип+1 - ип) - xn_1/2(un - u„_i)]. (27)

Это жесткая система нелинейных ОДУ. Чем больше т, тем большая надежность требуется от численного метода. Здесь в тестовых расчетах использовались т = 5, с = 1, но = 4.

В системе (27) фигурируют значения коэффициента теплопроводности в серединах интервалов хп±1/2. Для его вычисления будем использовать аппроксимацию Xn+i/2 = (х(ип+1) + х(ип)) /2. Тогда система (27) принимает следующий вид:

Для расчета задачи (28) будем использовать следующий набор схем. Известная обратная схема Эйлера (OIRK1), (й — и) /т = f(u), входит в семейство оптимальных обратных схем Рунге-Кутты (4) при s = 1. Схема Ll-устойчива, имеет точность лишь 0(т), но обладает высокой надежностью. Схема Crank-Nicolson (CN), также известная как схема "с полусуммой", (й — и) /т = (f(u) + f(u)) /2, имеет точность 0(т2), но при этом она лишь Л-устойчива и немонотонна. Обратная схема (7) (BMP), как было показано выше, является Z/2-устойчивой и имеет точность 0(т2). Из семейства явно-неявных комплексных схем Розенброка возьмем следующую Ь2-устойчивую схему точности 0(т2) (CROS):

^Е — ~~~~~Tfu j w = f(u), u = u + rRew. (29)

2 J 15

Схемы порядка точности р > 2 не тестировались, поскольку пространственный оператор (28) имеет лишь 2-й порядок аппроксимации.

Расчеты показали, что безытерационные явно-неявные разностные схемы, такие, как схема CROS, или неявные схемы, в которых итерации выполняются не до сходимости (число итераций ограничено), не дают сходимости к решению задачи (28): счет становится немонотонным и "разваливается". Неявные схемы OIRK1 и CN, в которых выполняется только одна итерация, продвигаются с той же скоростью, что и безытерационная схема CROS, и все сильнее отстают

от точного решения.

Надежный результат дает только расчет с итерациями, выполненными до сходимости. Расчет задачи (28) при соотношении шагов т/h = 2 без ограничения количества итераций для схем BMP, OIRK1 и CN дает гладкий профиль, визуально близкий к точному решению.

Рис. 2: Профиль бегущей волны для момента времени Ь = 0.2. Жирная линия - точное решение. Маркерами показаны расчеты для г/Л. = 4: В - СТ<Г с 1 итерацией, □ - СИ с 3 итерациями, . - СЫ с 6 итерациями, о - СЫ с итерациями до сходимости.

Более подробно результаты расчетов с ограниченным числом итераций для схемы СЫ представлены на рис.2. Как было упомянуто выше, счет с одной итерацией сильно отстает от точного решения. При ограничении в 2-5 итераций счет разваливается. Ограничение в б итераций дает ответ близкий к искомому решению, однако он является немонотонным. Заметим, что для отношения т/Л. = 4 шести итераций недостаточно для получения гладкого решения, хотя численное решение и не отстает так сильно от точного. Приведенные результаты иллюстрируют справедливость следующей теоремы.

Теорема 2. Пусть для уравнения (28) составлена любая разностная схема и фиксированы число итераций и отношение т/h. Всегда найдется такое точное решение, сходимости к которому не будет при г, h —)■ 0.

При увеличении соотношения шагов т/h счет с использованием ньютоновских итераций становится все менее устойчивым и требует все больше итераций до вхождения в малую окрестность решения и сходимости. Это связано с использованием профиля бегущей волны с предыдущего момента времени в качестве начального приближения для нового слоя. Чем больше соотношение т/h, тем худшим оказывается такое начальное приближение. Расчеты показали, что в таких условиях усечение (11) помогает значительно сократить количество итераций до сходимости во всех тестируемых схемах с итерациями.

Транзисторный усилитель. Расчет транзисторного усилителя является одной из традиционных задач для тестирования разностных схем решения дифференциально-алгебраических систем. Система, описывающая электрическую схему этого тринзисторного усилителя, представляет собой дифференциально-алгебраическую систему размерности M = 8, записанную в неявной форме (12). Ставится задача найти напряжение на выходе усилителя при заданных значениях параметров элементов цепи и периодическом входном напряжении.

Проведем расчет на последовательности сгущающихся сеток с помощью оптимальных обратных схем (13), одностадийной схемы Розенброка (29) и следующей двухстадийной схемы Розенброка с комплексными коэффициентами (CROS4):

(G-rcnfJu^k^fCu),

(G - та2{и (и + т Re a21k:)) k2 = f (u + т Re caifci)), (30)

û = u + т Re (biki + b2k2),

На рис. 3 приведены графики абсолютной погрешности решения в норме 1/2 в зависимости от числа узлов сетки для всех упомянутых схем. Видно, что для схем CRÛS и OIRK2 достигается второй порядок точности, и схема OIRK2 дает небольшой выигрыш в числе узлов. Схема Розенброка CROS4 дает разумное численное решение, близкое к точному, лишь на более мелкой сетке по сравнению с OIRX3 и другими схемами. Однако при дальнейшем сгущении она позволяет получить более точное решение.

Единственная схема 4-го порядка точности OIRK4 выходит на теоретический порядок точности лишь на последних четырех сетках. При этом, как было отмечено выше, время расчета по обратным схемам быстро увеличивается с ростом стадийности схемы и размерности задачи. При решении задачи о транзисторном усилителе время расчета по явно-неявным схемам Розенброка

17

Рис. 3: Погрешность численного решения задачи о транзисторном усилителе. Схемы: . -CROS, о - OIRK2, В - CROS4, □ - OIRK3, Д - OIRK4.

было в 10-50 раз меньше времени расчета по оптимальным обратным схемам Рунге-Кутты аналогичного порядка точности.

Уравнение ван дер Пола. Рассмотрим еще одну задачу, которая получается при моделировании электрической цепи с нелинейным элементом. Электрическая схема, впервые составленная Бальтазаром ван дер Полом, представляет собой RLC-контур, в котором пассивный резистор заменен на активный элемент, усиливающий ток в цепи, если его амплитуда падает ниже определенного уровня. Данная цепь описывается следующим обыкновенным дифференциальным уравнением второго порядка:

г" + o(z2 - l)z' + z = 0, (31)

где г пропорционально значению тока в контуре 1(f), а значение о зависит от параметров элементов цепи. Обозначив z = и, z1 = v, получим систему обыкновенных дифференциальных уравнений первого порядка:

du/dt = V, dv/dt = -u-o(u2 -l)v. (32)

Для численного интегрирования системы (32) использовались следующие разностные схемы. Из семейства явных схем Рунге-Кутгы были взяты оптимальные явные схемы (3) 2-го и 4-го порядка точности ERK2 и ERK4. Обе эти схемы не могут иметь А-устойчивости, и считаются непригодными для жестких задач. Из семейства явно-неявных схем Розенброка были взяты схемы CROS (29) и CROS4 (30). Из семейства полностью неявных схем Рунге-Кутты были взяты схемы OIRK1 и BMP (7). Эти схемы являются наиболее надежными для расчета

жестких задач благодаря ньютоновским итерациям до сходимости.

18

Расчеты проводились как при аргументе "время" ;£, так и при аргументе "длина дуги" I на последовательности равномерных сеток, рекуррентно сгущающихся вдвое. Сравнение разностных схем проводилось при сх = 100. Результаты расчета погрешности представлены на рис. 4. На нем по оси абсцисс как для аргумента £, так и для I нанесено количество узлов для одного цикла. Обсудим эти результаты.

Рис. 4: Задача (32) для а = 100. Маркеры черные для аргумента t и светлые для аргумента 1\ квадратики для ERK2 и ERK4, кружки для CROS и CROS4, треугольники для схемы BMP.

Верхняя линия есть слияние трех кривых для схем второго порядка при аргументе t. На первый взгляд это кажется странным: явная и неявные схемы дали одинаковую точность, хотя явные схемы принято считать плохими для жестких задач. Причина по всей видимости в том, что задача (32) не является чисто жесткой. В ней участки плохой обусловленности и комплексности спектра матрицы Якоби оказывают не меньшее влияние на расчет, чем участки собственно жесткости. Неявные же схемы эффективны на жестких участках решения. Третья кривая соответствует тем же схемам второго порядка точности, но при аргументе I. Опять эти схемы дают практически одинаковые результаты. Но теперь погрешность в ~ 106 раз меньше, чем при аргументе t. Для схем 4-го порядка выводы качественно аналогичны: разница в точности при одинаковых числах узлов по аргументам t или I составляет ~ 105.

Проиллюстрируем влияние жесткости задачи на результат расчета с помощью длины дуги. Возьмем для примера явную схему ERK4 (на остальных схемах получаются сходные результаты). В расчетах брались значения параметра от а = 1 до а = 1000. Первое соответствует мягкой задаче, последнее - очень жесткой. Результаты представлены на рис. 5.

10""'— 10

ю' ю'

N

10' Ю"

Рис. 5: Задача (32), схема ЕЛ/С4; В - расчеты по □ - расчеты по I. Около линий указаны

значения ст.

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

Химическая кинетика. Для уравнения ван дер Пола характерной константой жесткости является величина ст. Задача с а = 100...1000 считается достаточно трудной для тестирования методов решения жестких систем (задача ван дер Пола с параметром жесткости <7 = 104 входит в набор стандартных процедур системы МАТЬАВ). Однако имеются важные классы задач, в которых параметр жесткости может превышать Ю10. Например, такими являются задачи химической кинетики, поскольку скорости различных реакций в них могут отличаться на много порядков. Такие задачи можно назвать сверхжесткими (хотя правильнее назвать их сверхтрудными, поскольку в них есть как процессы сверхбыстрого затухания, так и процессы очень быстрого нарастания). Кроме того, практические задачи приводят к системам уравнений большой размерно-

сти.

Рассмотрим тестовую задачу, описывающую процесс горения природного газа в воздухе. Данная модель включает 25 химических уравнений и 20 компонент. Искомые концентрации этих компонент имеют сильно различающиеся масштабы, от ~ 10'2 до ~ 10"17. Система ОДУ, описывающая данную задачу, имеет канонический вид (1). В ней элементы правой части £ представляют собой линейные комбинации вспомогательных переменных г,, выражения для которых имеют вид г,- = Ь • и, ■ «т или г,- = Ь ■ щ. Здесь к} - константы реакций, щ,ит-

Ч - ^

компоненты вектора решения и. Начальные данные ио для каждой компоненты берутся по окончании основного горения и соответствуют догоранию. Интервал интегрирования по времени возьмем 0 < t < 1.2.

Явные схемы Рунге-Кутты оказались непригодными для решения таких задач. При интегрировании по времени счет разваливается на первых же шагах даже при очень маленьких значениях шага ht ~ Ю"8: значения концентраций по модулю возрастают и превосходят представимые на компьютере числа.

Переход к длине дуги меняет характер трудностей. Первый шаг после начального приближения дает правдоподобное решение. Однако на дальнейших шагах время t. практически не возрастает. Расчетные графики концентраций становятся пилообразными. При этом амплитуда пилы может быть на порядки больше, чем точные значения концентраций. Из-за этого концентрации могут принимать даже отрицательные значения, что химически бессмысленно. Большие концентрации приводят к большим значениям правых частей f(u). Они входят в знаменатель правой части F для функции t(l). В результате эта правая часть оказывается исчезающе малой, и увеличения расчетного t от шага к шагу практически не происходит.

Эти соображения справедливы при вычислениях с бесконечной разрядностью. При конечной разрядности вычислений возникает одна тонкость. Расчеты с сильной разномасштабностью констант реакций невозможно выполнять при 32-разрядных числах: ошибки округления становятся сопоставимыми с самим решением. При 64-разрядных числах можно производить хорошие расчеты, если матрица Якоби системы записывается аналитически. Но если матрица Якоби вычисляется разностно, то для нее необходимо использовать 128-разрядные числа (решать же получающуюся линеаризованную систему можно и с 64-разрядными числами). Последняя тонкость особенно существенна, если аргументом является длина дуги.

Расчеты тестовой задачи проводились как для аргумента t, так и для аргумента I. Хорошие результаты в обоих случаях дала схема CROS, где матрица Якоби вычислялась разностно с использованием 128-разрядных чисел. Остальные вычисления проводились с 64-разрядными числами, что при решении вспомогательной линейной системы давало существенную экономию времени.

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

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

В пятой главе приведено описание и часть исходного кода пакета GEABORK ("Guaranteed Error for Arc Backward Optimal Runge-Kutta methods"), который был создан в процессе разработки и отладки численных методов, представленных в данной работе. Он написан на языке программирования MATLAB, полностью совместимом со свободно распространяемой средой для математических вычислений Octave. Комплекс включает в себя широкий выбор численных методов для интегрирования жестких и плохо обусловленных систем ОДУ и дифференциально-алгебраических систем, а также необходимые вспомогательные подпрограммы (обобщенный метод Ньютона для решения нелинейных алгебраических систем, определение погрешности численного решения по методу Ричардсона, разностное вычисление матрицы Якоби, методы автономизации задачи и т.д.). Все компоненты комплекса GEABORK сопровождаются автоматизированными тестами для быстрой проверки правильности реализации.

Основные результаты

1. Построен новый класс неявных схем - оптимальные обратные схемы Рунге-Кутгы с числом стадий s < 4. Доказано, что эти схемы Ls-устойчивы и сходятся с 5-тым порядком точности как на обыкновенных дифференциальных, так и на дифференциально-алгебраических системах индекса 1. По сочетанию устойчивости, точности и экономичности они являются наилучшими среди всех схем типа Рунге-Кутты. Эти схемы эффективны для решения жестких и сверхжестких задач.

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

22

дуги улучшает расчет пограничных слоев и также позволяет вычислять апостериорную оценку погрешности.

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

4. С помощью пакета GEABORK проведено моделирование бегущей тепловой волны, транзисторного усилителя и процесса горения метана в воздухе.

Публикации автора по теме диссертации

1. H.H. Калиткин, И.П. Поишвайло, Обратные Ls-устойчивые схемы Рунге-Кутгы // ДАН, 2012 г., т.442, №2, с.175-180.

2. H.H. Калиткин, И.П. Пошивайло, Гарантированная точность при решении задачи Коши методом длины дуги // ДАН, 2013 г., т.452, №5, с.499-502.

3. И.П. Поишвайло, Усеченный многомерный метод Ньютона // Математическое моделирование, 2012 г., т.24, №1, с.103-108.

4. H.H. Калиткин, И.П. Пошивайло, Вычисления с использованием обратных схем Рунге-Кутты // Математическое моделирование, 2013 г., т.25, №10, с.79-96.

5. H.H. Калиткин, И.П. Пошивайло, Решение задачи Коши для жестких систем с гарантированной точностью методом длины дуги // Математическое моделирование, 2014 г., т.26, №7, с.3-18.

6. H.H. Калиткин, И.П. Поишвайло, Определение кратности корня нелинейного алгебраического уравнения // ЖВМиМФ, 2008 г., т.48, №7, с.1181-1186.

7. H.H. Калиткин, И.П. Пошивайло, О вычислении простых и кратных корней нелинейного уравнения // Математическое моделирование, 2008 г., т.20, №7, с.57-64.

8. A.A. Белов, H.H. Калиткин, И.П. Пошивайло, Численные методы решения сверхжестких задач Коши // Международная конференция «Современные проблемы вычислительной математики и математической физики» памяти академика A.A. Самарского к 95-летию со дня рождения. Москва, 16-17 июня 2014 года. Тезисы докладов. С. 27-29.

23

Подписано в печать: Формат 60x84 1/16. Уч.-изд.л. Тираж 80 экз. Заказ №36 Отпечатано в типографии ИПК МИЭТ.

124498, г. Москва, г. Зеленоград, площадь Шокина, дом 1, МИ