CC BY
119
УДК 519.622
МАКСИМАЛЬНЫЙ ПОРЯДОК ТОЧНОСТИ (ш, 1)-МЕТ0Д0В РЕШЕНИЯ ЖЁСТКИХ ЗАДАЧ
Е. А. Новиков1,2
1 Институт вычислительного моделирования СО РАН,
660036, Красноярск, Академгородок.
2 Сибирский федеральный университет,
660041, Красноярск, пр. Свободный, 79.
E-mail: novikov@icm.krasn.ru
Исследованы (то, 1)-методы решения жёстких задач, в которых на каждом шаге один раз вычисляется правая часть системы дифференциальных уравнений. Показано, что максимальный порядок точности L-устойчивого (rn, 1)-метода равен двум, и построен метод максимального порядка.
Ключевые слова: жёсткие задачи, схемы Розенброка, (то, к)-методы, А-устой-чивость, Ь-устойчивость.
Введение. При решении задачи Коши для жёстких систем обыкновенных дифференциальных уравнений широкое распространение получили методы типа Розенброка [1] благодаря простоте реализации и достаточно хорошим свойствам точности и устойчивости. Данные численные схемы получены из полуявных методов типа Рунге—Кутта, в которых для решения нелинейной системы алгебраических уравнений, возникающей при вычислении каждой стадии, используется одна итерация метода Ньютона. Все остальные проблемы решаются выбором величины шага интегрирования.
Наибольшее распространение получили методы типа Розенброка, в которых при вычислении каждой стадии применяется одна и та же матрица Якоби. Известно [2], что в этом случае для m-стадийного метода Розенброка максимальный порядок точности равен т + 1, причем схема максимального порядка может быть только А-устойчивой. Если отказаться от максимального порядка, то можно построить L-устойчивую численную формулу т-того порядка точности. В практических расчётах, как правило, отказываются от максимального порядка в пользу L-устойчивости. Заметим, что на основе методов типа Розенброка нельзя построить схему с замораживанием матрицы Якоби выше второго порядка точности [3], что ограничивает применение данных методов расчётами с небольшой точностью или задачами небольшой размерности.
В [4] предложен класс (т, к)-методов, в которых нахождение стадий не связывается с обязательным вычислением правой части системы дифференциальных уравнений. Числа тик означают соответственно число стадий и количество вычислений правой части системы дифференциальных уравнений на шаг интегрирования. Реализация (т, к)-методов так же проста, как и методов Розенброка, однако (т, к)-схемы имеют лучшие свойства точности и устойчивости. В рамках (т, к)-методов значительно проще решаются
Евгений Александрович Новиков (д.ф.-м.н., профессор), главный научный сотрудник, отд. вычислительной математики1; зав. кафедрой, каф. математического обеспечения дискретных устройств и систем2.
проблемы замораживания матрицы Якоби и ее численной аппроксимации.
Здесь исследуются (т, 1)-методы решения жёстких задач, в которых на каждом шаге один раз вычисляется правая часть системы дифференциальных уравнений. Показано, что максимальный порядок точности Ь-устойчи-вого (т, 1)-метода равен двум, и построен метод максимального порядка.
1. Схемы типа Розенброка. Далее будет рассматриваться задача Коши для автономной системы обыкновенных дифференциальных уравнений
У1 = /(У), У(Ро) =Уо, *о < * < (1)
где у и / — вещественные Ж-мерные вектор-функции, £ — независимая переменная. Рассмотрение автономной задачи не снижает общности, потому что введением дополнительной переменной у'м+1 = 1, УN+1(^0) = £о неавтономную задачу можно привести к автономному виду. Для решения задачи (1) будем применять методы типа Розенброка вида
гп / г—1 \
Уп+1 =Уп + ^РгЬ, АЛ = Л,/ ( уп + ^ (Зф-) ) , (2)
1=1 ' 3=1 '
где А. = Е — аН/'п\ Е — единичная матрица; /!п = 9/(уп)/ду — матрица Яко-
би системы (1); кг— стадии метода; а, Рг, /3^ — коэффициенты, определяющие свойства точности и устойчивости (2); 1 ^ г ^ т, 1 ^ ^ г — 1. В настоящее
время методы типа Розенброка трактуются более широко. Под ними понимаются все численные схемы, в которых матрица Якоби или ее аппроксимация вводятся непосредственно в формулу интегрирования.
Рассмотрим в качестве примера одностадийный метод типа Розенброка
Уп+1 = Уп + Р1к1, АЛ = Ь/(уп), (3)
где матрица А определена в (2). Разлагая приближенное решение уп+1 в ряд Тейлора по степеням Л, до членов с Л2 включительно, получим
Уп+1 = Уп+рФ!п + арф2!п!п + ОО3),
где /п = /(уп)• Представление точного решения у{Ьп+1) в виде ряда Тейлора в окрестности точки Ьп имеет вид
У^п+1) = у(1п) + Л./ + ^2/7 + оо3),
где элементарные дифференциалы / и /;/ вычислены на точном решении у(гп)- Сравнивая ряды для точного и приближенного решений при условии Уп = У^п), видим, что схема (3) будет иметь второй порядок точности, если Р\ = 1 и ар\ = 0,5, то есть р\ = 1 и а = 0,5. Теперь исследуем устойчивость метода (3). Для этого применим его для решения скалярного тестового уравнения у1 = Ху, где Л есть произвольное комплексное число, 11е(А) < 0. Смысл Л — некоторое собственное число матрицы Якоби задачи (1). Обозначая X = 1г\, получим уп+\ = Я(х)уп, где функция устойчивости С}(х) имеет следующий вид:
1 + (Р1 ~ а)х 1 — ах
Подставляя сюда значения коэффициентов р\ = 1 и а = 0,5, имеем С}(х) = = (1 + 0,5ж)/(1 — 0,5ж), то есть схема (3) второго порядка является А-устой-чивой. Из вида функции С,)(х) следует, что схема (3) будет Ь-устойчивой, если Р\ = а = 1, что противоречит второму порядку точности. Обычно отказываются от второго порядка в пользу Ь-устойчивости, что приводит к более эффективному методу, хотя и первого порядка.
В случае большой размерности задачи (1) основные вычислительные затраты связаны с обращением матрицы Аг- Обычно вместо обращения решается линейная система алгебраических уравнений АЛ = Ь/(уп) с применением Ьи-разложение матрицы А. с выбором главного элемента по строке или столбцу, а иногда и по всей матрице, то есть при вычислении приближенного решения по формуле (3) осуществляется декомпозиция матрицы А. (порядка -/V3 арифметических операций). Обратный ход метода Гаусса стоит порядка -/V2 операций. Таким образом, при большой размерности исходной задачи общие вычислительные затраты определяются временем декомпозиции матрицы Аг- Возникает естественный вопрос: нельзя ли без значительного увеличения вычислительных затрат исправить схему (3) таким образом, чтобы она была Ь-устойчивой и имела второй порядок точности? Эта проблема решается в рамках (т, /с)-методов.
2. Класс (т, /с)-методов. Пусть заданы т, к € Ъ, к ^ т. Обозначим
через Мт множество чисел ! е Z таких, что 1 ^ г ^ т, а через Мд. и ^
подмножества из Мт вида
Мк = {гПг € Мт | 1 = Ш1 < Ш2 < • • • < ГПк ^ т},
^ = {mj-l € Мт \ j > 1, mj € Мк, mj ^ г}, 1 < г ^ т.
Рассмотрим следующие численные схемы:
т
Уп+1 =Уп + У~]рЛ, г=1
(г—1 \ г—1
Уп Н- ^ ^ | Н- ^ ^ + Ь/п ^ ^ ^ Мк 1 (4)
з=± ' З^г з=1
г—1
АЛ = кг-1 + ^ агуку + Ь/'п ^ с^ку, г € Мт \ Мк,
3= 1
где А = Е — аЛ,/^; а, Рг, /3^, а^, — постоянные коэффициенты; Л, —шаг
интегрирования; к и т— соответственно количество вычислений функции / и число обратных ходов в методе Гаусса (число стадий). На каждом шаге интегрирования осуществляются одно вычисление матрицы Якоби и одна декомпозиция матрицы Аг- Допускается аппроксимация матрицы Якоби матрицей Ап, удовлетворяющей условию
Ап = /га + Ь-Вп + 0(Л,2),
где матрица Вп не зависит от шага интегрирования. Данное условие позволяет применять методы (4) с замораживанием как численной, так и аналитической матрицы Якоби. Так как кит полностью определяют затраты на
шаг, а набор чисел т\, тг, • •., тк из множества Мд. только распределяет их внутри шага, то методы типа (4) называются (т, /с)-методами.
Отметим, что при к = т и = 0 методы (4) совпадают со схемами
типа Розенброка (2), а при к = т и ац = 0 —с 1Ю\^-методами [2]. Заметим также, что при рассмотрении методов такого типа все авторы изучали случай к = т, то есть когда число стадий и количество вычислений правой части системы дифференциальных уравнений (1) совпадают. В этом случае /с-стадийную схему (4) можно поставить в соответствие /с-стадийной полуяв-ной формуле типа Рунге—Кутта, при реализации которой на каждом шаге используется одна матрица размерности N. Относительно таких численных формул известно, что нельзя построить /с-стадийную схему выше (к + 1)-го порядка точности, причем схема максимального порядка является А-устой-чивой. Очевидно, этот факт распространяется и на (т, /с)-методы при т = к.
3. Численные схемы с одним вычислением правой части. Выберем Мт = = {1,2..., т} и Мк = {1} при к = 1, тогда ■]{ есть пустое множество. Рассмотрим семейство методов следующего вида:
т
Уп+1=Уп + ^'Ргкг, В>пк\ = к/(уп), Опкг = кг-1, 2 ^ г ^ т, (5)
г=1
где матрица Вп определена в (4). Для изучения (5) введём в рассмотрение матрицу В с элементами Ь^:
Ьц — Ьц — 1, 2^1, Ьц — ^ Ьг^— 1, I,] ^ 2. (6)
Лемма 1. Элементы матрицы В представимы в виде
з г
Ьц = ^ 1,к1 Ьц = ^ ' Ък^ — \, ъ,] ^ 2. (7)
к=1 к=1
Доказательство. Для доказательства достаточно расписать второе рекуррентное соотношение (6), используя первое условие. □
Ниже через В3гк будем обозначать матрицу с элементами (6), составленную из первых 8 строк и к столбцов матрицы В.
Лемма 2. Матрица Вт>т невырожденная.
Доказательство. Для доказательства введем в рассмотрение матрицы Кк, 2 ^ к ^ т, с элементами гг£, у которых все элементы равны нулю, за исключением следующих:
г%к = 1, 1 ^ г ^ т, гг£г~1 = —1, к ^ г ^ т, (8)
и матрицу К вида К = ИтИт-1... К.2- Очевидно, матрица К невырожден-
ная. Покажем, что ЯВт>т есть верхняя треугольная матрица с единицами на главной диагонали, что доказывает лемму 2.
Сначала умножим _Й2 на Вт>т. Учитывая (8), для этого нужно из второй строки матрицы Вт^т вычесть первую, из третьей вторую и т. д., а результат
следует поместить на место соответственно второй, третьей и т.д. строк. Тогда с использованием первого соотношения (6) получим, что в первом столбце матрицы -Й2-Вт,т на первом месте стоит единица, а на остальных нули. Из второго равенства (6) имеем
Ьг^-1 = Ьц - Ьг-1^', 2^г^т, 3 ^ ^ т,
ТО есть, если В матрице ДгАта.т. вычеркнуть первую строку и столбец, то она совпадает с Дп,т> У которой вычеркнуты первая строка и последний столбец. Умножая последовательно матрицу ЕгАп,™ на Из, ..., Кт и проводя аналогичные рассуждения, получим доказательство леммы. □
Заметим, что так как матрицы Вт^т и Кь, 2 ^ к ^ т, целочисленные, то К и КВт>т также целочисленные.
Лемма 3. Стадии метода (5) представимы в виде
ОС
к, = Ху-^/п^/п. (9)
г=1
Доказательство. Доказательство проведем с использованием метода математической индукции. Для этого запишем И~1 в виде ряда Тейлора
Оп 1 = Е + ак/'п + а2Н2/'2 + а3И3 + — (10)
Представление (9) при ] = 1 очевидно. Пусть (9) выполняется при некотором ]. Для определения ^+1 умножим (10) на (9) и соберем слагаемые при одинаковых степенях Н. В результате имеем
Воспользовавшись вторым равенством (7), завершим доказательство. □
С применением (9) приближенное решение по схеме (5) представимо в виде
сю / т \
Уп+1 = Уп + ^ а%~х ( Е ЪаРз ) Кй1-Х)!п. (11)
1=1 \? = 1 '
Теорема. При любом значении т, нельзя построить т-стадийную схему
(5) выше второго порядка точности.
Доказательство. Запишем ряд Тейлора для точного решения у (^+1) задачи (1) в окрестности точки Ьп по степеням Н до членов с Л3 включительно, то есть
у^п+о = у(ьп) + л,/ + ^2/7 + ^3(//2/ + /;72) + 0(/г4), (12)
где элементарные дифференциалы вычислены на точном решении у(Ьп). Положим уп = у(Ьп). Тогда из сравнения (11) и (12) следует, что в разложении
точного решения в ряд Тейлора есть слагаемое вида //;/2> а в представлении (11) оно отсутствует. □
Отметим, что в случае линейной задачи (1) вида
где А и Ь — соответственно матрица и вектор с постоянными коэффициентами, имеет место соотношение
с невырожденной матрицей Вт<т, которое можно применять для построения методов заданного порядка точности. Здесь т —число стадий в методе (5),
4. Метод второго порядка точности. Пусть т = 2, то есть рассмотрим схему вида
Подставляя ряды Тейлора для к\ и к2 в первую формулу (13), получим
Полагая уп = у(£га) и сравнивая полученное соотношение с (12) до членов с Л2 включительно, получим условия второго порядка точности (13), то есть
Исследуем устойчивость (13). Для этого применим его для решения скалярного тестового уравнения у' = Ху. Учитывая условия порядка, получим уп+1 = = Я{х)уп, где функция устойчивости С}(х) имеет следующий вид:
Тогда схема (13) будет Ь-устойчивой, если р\ = а. Решая систему р\ = а и (14), имеем набор коэффициентов р\ = а и р2 = 1 — а, где а удовлетворяет условию Ь-устойчивости:
Данное уравнение имеет два корня: а\ = 1 — л/2/2 и й2 = 1 + л/2/2. Обычно в расчётах применяется корень а = а\, потому что в этом случае меньше коэффициент в локальной ошибке.
Контроль точности вычислений численной схемы (13) построим по аналогии [5]. Для этого введём обозначение
у' = Ау + Ь, у(і0) = уо,
Уп+1 = Уп + Рікі + р2к2, Опкі = к/(уп), Опк2 = кі. (13)
Уп+і =уп + (Рі +Р2)Ь/п + а(рі + 2 р2)Л.2/4/п+
+ а2(рі + Ър2)Ъ?/п/п + 0(Л,4).
£>і+£>2 = 1, а(рг+2р2) = 0,5.
(14)
Сі(х)
1 + (1 — 2 а)х + а(а — р\)х2 (1 — ах)2
а2 — 2а + 0,5 = 0.
є(.іп) = 01п ]п(к2-кі)
(15)
где к\ и к2 вычисляются по формулам (13). Тогда согласно [5] для контроля точности вычислений на каждом шаге нужно проверять неравенство
1к0п)1Ке, 1<;'п<2, (16)
где е — требуемая точность расчётов, || • || —некоторая норма в Мм, а целочисленная переменная ]п выбирается наименьшей, при которой выполняется неравенство (16).
Отметим одну важную особенность построенной оценки ошибки (15). Схема (13) Ь-устойчивая, то есть для ее функции устойчивости С,)(х) имеет место соотношение С}(х) —> 0 при х У оо. Так как для точного решения у(Ьп+1) = = ехр(ж)у(£га) задачи у' = Ху, у{Ъо) = уо выполняется аналогичное свойство, то естественным будет требование стремления к нулю оценки ошибки при х —>■ —оо. Однако для е(1) это свойство не выполняется — данная оценка ведет себя А-устойчивым образом. С целью исправления асимптотического поведения ошибки вместо е(1) введена оценка £^п), 1 ^ ]п ^ 2. В этом случае поведение оценки ошибки при ]п = 2 будет согласовано с поведением точного решения тестовой задачи у' = Ху, у(£о) = Уо при х —> — оо.
Подчеркнем, что в смысле главного члена оценки е(1) и е(2) совпадают. Использование е{]п) фактически не приводит к увеличению вычислительных затрат. Это связано с тем, что е{]п) при ]п = 2 проверяется только в том случае, если оно нарушено при ]п = 1. Такая ситуация встречается достаточно редко, в основном при быстром росте величины шага интегрирования. Однако это позволяет выбирать шаг более правильно и тем самым уменьшается количество неоправданных повторных вычислений решения (возвратов).
В расчётах норма ||£|| в левой части неравенства (16) вычисляется по формуле
||£|| = шах —.
\угп\ + IX
В случае выполнения неравенства |у^| < ц по г-той компоненте решения контролируется абсолютная ошибка /хе, в противном случае контролируется относительная ошибка е. Иногда с целью повышения надежности расчётов задают набор параметров /^, 1 ^ г ^ Ж, что позволяет более квалифицированно контролировать точность расчётов.
5. Заключение. Из сравнения численной формулы типа Розенброка (3) и (2,1)-метода (13) следует, что по вычислительным затратам схема (13) отличается от (3) на один дополнительный обратный ход метода Гаусса. Однако (14) имеет второй порядок точности и, как показывают расчёты, примерно в три раза эффективнее (3) по времени счета. Кроме того, для (14) построено неравенство для контроля точности вычислений, что позволяет проводить расчёты с переменным шагом интегрирования. Это приводит к дополнительному повышению эффективности.
Численную схему (13) можно рассматривать как способ реализации неявного одностадийного метода типа Рунге—Кутта, причем в (13) нет необходимости применять итерационный процесс Ньютона, что позволяет оценить вычислительные затраты на шаг интегрирования до начала расчётов и значительно упрощает реализацию алгоритма интегрирования.
Работа выполнена при поддержке РФФИ (проект № 11-01-00106-а).
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Rosenbrock Н. Н. Some general implicit processes for the numerical solution of differential equations// Computer, f963. Vol. 5, no. 4. Pp. 329-330.
2. Hairer E., Wanner G. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems / Springer Series in Computational Mathematics. Vol. 14. Berlin: Springer-Verlag, 1996. 614 pp.; русск. пер.: Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жёсткие и дифференциально-алгебраические задачи. М.: Мир, 1999. 685 с.
3. Новиков Е.А., Двинский А. Л. Замораживание матрицы Якоби в (3, 2)-методе решения жёстких систем/ В сб.: Совместный выпуск журналов «Вычислительные технологии» и «Региональный вестник Востока»: Труды международной конференции «Вычислительные и информационные технологии в науке, технике и образовании». Часть II. Новосибирск, Алматы, Усть-Каменогорск, 2003. С. 272-278. [Novikov Е.А., Dvinskiy A.L. Freezing of the Jacobi matrix in (3,2)-method of solving stiff systems/ In: Joint issue of “Computational Technologies” and “Regional Bulletin of the East Proceedings of International Conference “Computational and Informational Technologies for Science, Engineering and Education”. Part II. Novosibirsk, Almaty, Ust’-Kamenogorsk, 2003. Pp. 272-278].
4. Новиков E.A., Шитов Ю.А., Шокин Ю. И. Одношаговые безытерационные методы решения жёстких систем// ДАН СССР, 1988. Т. 301, №6. С. 1310-1314; англ. пер.: Novikov Е.А., Shitov Yu. A., Shokin Yu. I. One-step noniterative methods for solving stiff systems // Soviet Math. Dokl., 1989. Vol. 38, no. 1. Pp. 212-216.
5. Новиков E. А. Явные методы для жёстких систем. Новосибирск: Наука, 1997. 197 с. [Novikov Е.А. Explicit methods for stiff systems. Novosibirsk: Nauka, 1997. 195 pp.]
Поступила в редакцию 28/1/2011; в окончательном варианте — 17/VIII/2011.
MSC: 65L20; 65L05, 34A34
MAXIMAL ORDER OF ACCURACY OF (m, l)-METHODS FOR SOLVING STIFF PROBLEMS
E. A. Novikov1,2
1 Institute of Computational Modelling, Siberian Branch of the Russian Academy of Sciences, Akademgorodok, Krasnoyarsk, 660036.
2 Siberian Federal University,
79, Svobodniy, Krasnoyarsk, 660041.
E-mail: novikovSicm.krasn.ru
We investigate (to, 1)-methods for solving stiff problems in which the right part of system of the differential equations is calculated one times on each step. It is shown that the maximal order of accuracy of the L-stability (to, 1)-method is equal to two, and the method of the maximal order is constructed.
Key words: stiff problems, Rosenbrock schemes, (jn,k)-methods, A-stability, L-sta-bility.
Original article submitted 28/1/2011; revision submitted 17/VIII/2011.
Evgeniy A. Novikov (Dr. Sci. (Phys. & Math.)), Chief Research Scientist, Dept, of Computational Mathematics1; Head of Dept., Dept, of Mathematical Software for Systems and Discrete Devices2.
