CC BY
49
Серия «Математика»
2015. Т. 11. С. 3—11
Онлайн-доступ к журналу: http://isu.ru/izvestia
УДК 518.517
О потере Ь-устойчивости неявного метода Эйлера
О О О -А-
для одной линейной задачи *
М. В. Булатов
Институт динамики систем и теории управления им. В. М. Матросова СО РАН
Л. С. Соловарова
Институт динамики систем и теории управления им. В. М. Матросова СО РАН
Аннотация. Ряд важных прикладных задач из химической кинетики, биофизики, теории электрических схем описываются системами жестких обыкновенных дифференциальных уравнений с заданными начальными условиями. Одним из подходов для их численного решения являются одношаговые методы Рунге - Кутта. Для задач небольшой размерности применяют неявные методы Рунге - Кутта. Среди таких алгоритмов выделяют так называемые А- и Ь-устойчивые. Как правило, Ь-устойчивые гораздо лучше справляются с данными задачами. А именно, при реализации Ь-устойчивых методов шаг интегрирования можно выбрать значительно большим, чем при реализации А-устойчивых методов. Самым простым и хорошо себя зарекомендовавшим из данных алгоритмов является неявный метод Эйлера.
В данной статье приведен пример линейной автономной системы обыкновенных дифференциальных уравнений, зависящей от параметров, выбирая которые можно получить сколь угодно жесткую задачу. Показано, что при определенном выборе этих параметров неявная схема Эйлера оказывается неэффективной. Данный алгоритм будет устойчив только при существенном ограничении на шаг интегрирования. Построение данного примера основано на некоторых фактах из теории численного решения дифференциально-алгебраических уравнений высокого индекса. Приведены детальные выкладки.
Ключевые слова: жесткие ОДУ, дифференциально-алгебраические уравнения, разностные схемы, Ь-устойчивые методы.
Численное решение жестких обыкновенных дифференциальных уравнений (ОДУ) является классическим разделом вычислительной математики. К настоящему времени вышло и продолжает выходить огром-
* Работа выполнена при финансовой поддержке РФФИ, грант 15-01-03228.
1. Введение
ное количество работ по различным аспектам данной тематики. Весьма полная библиография приведена в монографиях [1] - [6], [9].
Для предсказания свойств разрабатываемых методов на протяжении десятилетий используется тестовое уравнение [10]
х = Лх(г),х(0) = хо,г е [0,1], (1.1)
где Л — скаляр и ЯеЛ << 0. Введем обозначения Н = 1/Ж, и = гН, г = 0,1, ...,Ы, XI & x(íi), г = ЛН. Тогда любой одношаговый метод для (1.1) можно представить в виде
Хг+1 = Я(г)Хг.
Я(г) принято называть функцией устойчивости (полином или дробно-рациональная функция).
Принята следующая классификация одношаговых методов. Если для любого г с Яег < 0 |Я(г)| < 1, то метод называют А-устойчивым. Если при этом Я(г) ^ 0 при г ^ те, то метод называется Ь-устойчивым.
Для численного решения начальной задачи жесткой системы ОДУ небольшой размерности, как правило, используют Ь-устойчивые разностные схемы. Простейшей из этих схем является неявный метод Эйлера, функция устойчивости которого
вд = 1
1- г
Возникает естественный вопрос: всегда ли эта разностная схема успешно справляется с жесткими задачами?
Ниже приведено построение жесткой системы ОДУ, для численного решения которой Ь-устойчивые методы, в частности, неявная схема Эйлера и двухстадийный метод Лобатто 111С требуют ограничения на шаг интегрирования. Построение данного примера основано на теории численного решения дифференциально-алгебраических уравнений (ДАУ). Необходимые сведения приведены в следующем параграфе.
2. Дифференциально-алгебраические уравнения
Рассмотрим задачу
А(г)х'(г) + в(г)х(г) = /(г), г е [0,1], (2.1)
х(0) = хо, (2.2)
где А(г), В(г) — (п х п)-матрицы, f (г) и х(г) — заданная и искомая п-мерные вектор функции.
Если йегА(г) = 0, то систему (2.1) принято называть ДАУ.
Предполагается, что для ДАУ (2.1) начальное условие (2.2) согласовано с правой частью, т. е. рассматриваемая задача имеет решение. Под решением здесь и всюду в дальнейшем понимается непрерывно-дифференцируемая вектор-функция, обращающая (2.1) в тождество и удовлетворяющая условию (2.2).
Характеристикой сложности ДАУ является понятие индекса, которое имеет несколько определений. В случае вещественно-аналитических коэффициентов матриц A(t), B(t) эти определения эквивалентны [7]. Приведем одно из них.
Определение 1. [7] ДАУ (2.1) имеет общее решение индекса r, если существуют (n х n)-матрицы &(t), K0(t,r), Ki(t),..., Kr(t), причем rank^(t) _ r _ const Vt £ [0,1], такие, 'что
x(t,c) = mc + J Ko(t,T)f (r)dr + Kj (t)f (-l)(t), (2.3) 0 j=l
является решением (2.1) и на любом интервале [а, в] ^ [0,1] нет решений, отличных от (2.3).
Например, ДАУ
w(t) 0\( U(t) \ , ( 0 1 \( u(t) \ _ ( g(t) \ 0 0) \v' (t)) + \l0j {v(t))_{ q(t)J
имеет решение вида (2.3) u _ q, v _ g — wq при любом w(t),q(t) £ Cl, g(t) £ C2 на любом интервале. Здесь $(t) _ K0(t,r) = 0 — нулевые матрицы, а ( \ ( \
ki_(01),к_i0 0
ч 1 0)' 2 V 0 -ш
ДАУ, похожее на предыдущее
о о)()+(V)=(о - о ч
имеет только тривиальное решение, однако на интервале [а, в], а > 0, в < 1, данное ДАУ имеет общее решение
V ) ^ 0 0) у С2
где С1,С2 — произвольные числа. Здесь нарушено условие постоянства ранга матрицы Ф(£) на отрезке [0,1]:
таикФ(Ь) = | | > 0°'1]
Численное решение ДАУ индекса 2 и выше достаточно сложная задача. Многие разностные схемы для таких задач либо принципиально неприменимы, либо неустойчивы. В качестве иллюстрации приведем пример [7]:
00 )( Я + ( 0 ° )( и ) = ( I > е 0 я (2.4)
где а - скалярный параметр, а = 1. Данное ДАУ имеет единственное решение индекса 2, которое не зависит от начальных данных: и = д—Ьь,
V = (д' — о)/(1 — а).
Проанализируем поведение неявного метода Эйлера для данного примера. Введем равномерную сетку на [0,1] tj = ]Н, Н = 1/И, ] = 0,1,...,И и обозначим Аj = А(tj),Bj = В(tj),fj = f(tj), Xj & x(tj),
xj = (и,vj)Т & (u(tj),ь(Ь))Т = x(tj).
Опуская несложные выкладки, получим, что
1! , дг+1 — д% \
а Н
Итак, при |а| < 1 данный метод неустойчив, а при а = 0 принципиально неприменим в силу того, что матрица А^+1 + НВ^+1 тождественно вырожденна.
Рассмотрим однородное ДАУ, похожее на (2.4), с параметром в, 0 < в << 1
Ч?t:,){:) ==0^е0и (2,)
Данная задача имеет индекс 1, и ее общее решение
и = + в)ь, V = С ехр((а — 1^/в),
где С — произвольное число.
Таким образом, при а < 1 и 0 < в << 1 получим жесткое ДАУ индекса 1. Неявный метод Эйлера для данного примера дает следующие рекуррентные соотношения:
иг+1 = —(и+1 + в ^+1, Н—в
=
Последнее соотношение будет устойчиво, если
1(Н — в)/(На — в)1 < 1, т.е.в этом случае мы получим ограничение на шаг
Л 6(0;-^-).
1 + а
Замечание 1. Пример жесткого ДАУ индекса 1, для которого неявная схема Эйлера требовала ограничения на шаг интегрирования, приведен в [11]. Данное явление детально проанализировано в статье [8].
В дальнейшем нам потребуется следующая
Теорема. [7] Если задача (2.1)-(2.2) имеет индекс 1, то \хе — х(£)\\ = О(е), где хе(Ь) является решением ОДУ
(Л(г) + ев(г))хе(г) + в(г)хе(г) = /(г),г - [0,1],х€(0) = хо,0 < е << 1.
В следующем разделе приведем анализ неявного метода Эйлера для одной жесткой линейной системы ОДУ. Данная система была построена на основе вышеизложенных результатов.
3. Поведение неявного метода Эйлера для одной жесткой
задачи
Применим для примера (2.5) результат теоремы. Получим
)(V)+(0«и)(V)=(0>-м ^
Данную системы можно переписать в виде
х + Бх = 0, (3.2)
т. е. в виде классической системы ОДУ. Здесь матрица
-1
. 1 * + в
Элементарные выкладки дают
и = —(ев — е2а)v — (I + в — еа)и, а компонента V удовлетворяет уравнению
(ев — е2а)ю' + (в — 2еа)у + (1 — а> = 0, (3.3)
где (и(г)^(г))Т.
Собственные числа матрицы Б(Ь,е, в, а) постоянны, удовлетворяют уравнению
(ев — е2а)Х2 + (2еа — в)Х — а = 0
Б = ^е-в.а) = (1 е^У (0 а
и равны
Ai = -,A2 = --^—. (3.4)
е р — еа
Выберем параметры а, (3 и е таким образом, чтобы корни характере-стического уравнения (3.3) были действительными и отрицательными, например а = — 1/4,в = 5е. В этом случае общее решение (3.3) имеет вид
u = —5.25е2у — (t + 5.25e)v,
15 7
v = C\ ехР(-7^) + c<i и собственные числа матрицы D(t,e, в, а) положительны
Ai = -, А2 = —.
е 21е
Неявный метод Эйлера для ОДУ (3.2) дает рекуррентные соотношения
Uj+i = —P (vi+i — vi) — (ti+i + /3 — еа^г+1, (—P + 2еа — в + ho^vi+i + (2P — h — 2еа + (3)vi — Pv—i = 0, (3.5)
г, ев—е2а
где P = .
Применим для разностного соотношения (3.5) критерий Рауса - Гур-
2 \ 2
вица, согласно которому корни характеристического уравнения a2л2 + bX + c = 0 лежат в единичном круге, если a < 0 и
1) a+b+c<0
2) a-c<0
3) a-b+c<0.
Первое условие данного критерия дает неравенство (а — l)h < 0, которое всегда справедливо при а < 1. Из второго условия следует неравенство —в + 2еа + hŒ < 0, которое также всегда справедливо при а < 0 и /в, h > 0. Третье неравенство критерия дает
(1 + а)h2 + (4ае — 2^)h — 4(е@ — е2а) < 0, (3.6)
Последнее неравенство при а G (—1,0) будет справедливо при
2/3 - 4ае + у/402 + 16е/3 - 16 Л€( 0,--). (3.7)
Таким образом, система ОДУ (3.2) при а = —1/4, /3 = 5е, 0 < е << 1 жесткая, и собственные числа матрицы D(t^, /3, а) действительные и
много больше нуля. Однако, из формулы (3.7) при таких значениях параметров вытекает ограничение на шаг
h 6 (0,22t
Если а < — 1, то при 0 < e,ß << 1 неявная схема Эйлера прекрасно справляется с данным примером. Это вытекает из того факта, что данный алгоритм сходится к точному решению с первым порядком точности для примера (2.4) при а < 1 и согласованности начальных условий.
В заключении, отметим, что многие двухстадийные методы Рунге -Кутта для ДАУ (2.4) при а = —1 дают вырожденную систему линейных алгебраических уравнений. Если для исходного ОДУ применить двух-стадийный метод Лобатто III C, то при а ^ 1 и ß = 6е (общее решение не содержит в этом случае осциллирующих компонент) метод теряет свойство L-устойчивости. Детальный анализ этого явления планируется провести в дальнейшем.
Список литературы
1. Арушанян О. Б. Численное решение обыкновенных дифференциальных уравнений на Фортране / О. Б. Арушанян, С. Ф. Залеткин. - М. : Изд-во Моск. ун-та. - 336 с.
2. Деккер К. Устойчивость методов Рунге - Кутты для жестких нелинейных дифференциальных уравнений / К. Деккер, Я. Вервер. - М. : Мир, 1988. -332 с.
3. Новиков Е. А. Компьютерное моделирование жестких гибридных систем / Е. А. Новиков, Ю. В. Шорников. - Новосибирск : НГТУ, 2012. - 450 с.
4. Ракитский Ю. В. Численные методы решения жестких систем / Ю. В. Ракит-ский, С. М. Устинов, И. Г. Черноруцкий. - М. : Наука. Гл. ред. физ.-мат. лит., 1979. - 208 с.
5. Хайрер Э. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи : пер. с англ. / Э. Хайрер, Г. Ваннер. - М. : Мир, 1999. - 685 с., ил.
6. Холл Дж. Современные численные методы решения обыкновенных дифференциальных уравнений / Дж. Холл, Дж. Уатт. - М. : Мир, 1979. -312 с.
7. Чистяков В. Ф. Алгебро-дифференциальные операторы с конечномерным ядром : монография / В. Ф. Чистяков. - Новосибирск : Наука, 1996. - 278 с.
8. Чистяков В. Ф. О сохранении типа устойчивости разностных схем при решении жестких дифференциально-алгебраических уравнений / В. Ф. Чистяков // Сиб. журн. вычисл. математики. - 2011. - Т. 14, № 4. - С. 443-456.
9. Butcher J. C. Numerical Methods for Ordinary Differential Equations / J. C. Butcher. - Wiley, 2008.
10. Dahlquist G. Convergence and stability in the numerical integration of ordinary differential equations / G. Dahlquist // Math. Scand. - 1956. - Vol. 4. - P. 33-53.
11. Marz R. Differential-algebraic systems anew / R. Marz // Appl. Numer. Math. -2002. - Vol. 42. - P. 315-335.
Булатов Михаил Валерьянович, доктор физико-математических наук, Институт динамики систем и теории управления им. В. М. Мат-росова СО РАН, 664033, Иркутск, ул. Лермонтова, 134 тел.: (3952)453018 (e-mail: mvbul@icc.ru)
Соловарова Любовь Степановна, аспирант, Институт динамики систем и теории управления им. В. М. Матросова СО РАН, 664033, Иркутск, ул. Лермонтова, 134 тел.: (3952)453029 (e-mail: soleilu@mail.ru)
M. V. Bulatov, L. S. Solovarova
On the Loss of L-stability of the Implicit Euler Method for a Linear Problem
Abstract.A number of important applied problems of chemical kinetics, biophysics, theory of electrical circuits are described by systems of stiff ordinary differential equations with given initial conditions. The one-step Runge-Kutta method is one of approaches for their numerical solution. The implicit Runge - Kutta methods are used for problems of small dimension. The so-called A- and L-stable methods are singled out among these algorithms. Usually, L-stable methods much better cope with these problems. Namely, when we implement L-stable methods we can choose the integration step much greater than in the implementation of A-stable methods. The implicit Euler method is the simplest of these algorithms and it well proved itself.
In the article we consider an example of a linear autonomous system of ordinary differential equations depending on parameters. By choosing these parameters, as is wished stiff problem can be obtained. It is shown that for a particular choice of the parameters Euler implicit method will be ineffective. It is stable only under significant restrictions of the integration step. The construction of this example is based on some facts of the theory of the numerical solution of differential-algebraic equations of high index. The detailed computations are shown.
Keywords: stiff ODE, differential-algebraic equations, difference schemes, L-stable methods.
References
1. Arushanyan O.B., Zaletkin S.F. Chislennoe reshenie obyknovennyx differentsial'nyh uravneniy [Numerical Solution of Ordinary Differential Equations Using FORTRAN]. Moscow, Mos. Gos. Univ., 1990. 336 p.
2. Dekker K., Verwer J.G. Stability of Runge - Kutta Methods for Stiff Nonlinear Differential Equations. North-Holland, Amsterdam, 1984. 332 p.
3. Novikov E.A., Shornikov Yu.V. Komp'yuternoe modelirovanie zhestkikh gibridnyh sistem[Computer modeling of stiff hybrid systems. Novosibirsk, Publishing house NGTU, 2012. 450 p.
4. Rakitskii Yu.V. Ustinov S.M., Chernorutskii I.G. Chislennye metody resheniya zhestkikh sistem [Numerical Methods for Solving Stiff Systems]. Moscow, Nauka, 1979. 208 p.
5. Hairer E., Wanner G. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Springer-Verlag, Berlin, 1996. 385 p.
6. Hall G., Watt J. Modern Numerical Methods for Ordinary Differential Equations. Oxford Univ., Oxford, 1976. 312 p.
7. Chistyakov V.F. Algebro-differentsial'nye operatory s konechnomernym yadrom[Algebraic Differential Operators with a Finite-Dimensional Kernel]. Novosibirsk, Nauka, 1996. 278 p.
8. Chistyakov V.F. Preservation of stability type of difference schemes when solving stiff differential algebraic equations. Numerical Analysis and Applications, 2011, vol. 4, issue 4, pp. 363-375.
9. Butcher J.C. Numerical Methods for Ordinary Differential Equations. Wiley, 2008.
10. Dahlquist G. Convergence and Stability in the Numerical Integration of Ordinary Differential Equations. Math.Scand., 1956, vol. 4, pp.33-53.
11. Marz R. Differential-algebraic Systems Anew. Appl. Numer. Math., 2002, vol. 42, pp. 315-335.
Bulatov Mikhail Valeryanovich, Doctor of Sciences (Physics and Mathematics), Matrosov Institute for System Dynamics and Control Theory SB RAS, 134, Lermontov st., Irkutsk, 664033 tel.: (3952)453018 (e-mail: mvbul@icc.ru)
Solovarova Liubov Stepanovna, Postgraduate, Matrosov Institute for System Dynamics and Control Theory SB RAS, 134, Lermontov st., Irkutsk, 664033 tel.: (3952)453029 (e-mail: soleilu@mail.ru)
