О блочных разностных схемах для дифференциально-алгебраических уравнений Текст научной статьи по специальности «Математика»

Серия «Математика»

2010. Т. 3, № 2. С. 2-12

Онлайн-доступ к журналу: http://isu.ru/izvestia

УДК 518.517

О блочных разностных схемах для дифференциально-алгебраических уравнений

М. В. Булатов

Институт динамики систем и теории управления СО РАН

Л. С. Соловарова

Институт динамики систем и теории управления СО РАН

Аннотация. В работе предложены блочные разностные схемы для численного решения начальной задачи линейных дифференциально-алгебраических уравнений. Проведен детальный анализ таких схем первого и второго порядка точности на модельных примерах и показано их преимущество над известными неявными методами первого и второго порядков. Приведено описание оБщго вида Блочных разностных схем.

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

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

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

1. Введение

ренцирования назад (ФДН методы) [1, 2], неявные методы Рунге-Кутта (РК) [1, 2] и блочные методы [3, 4].

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

2. Постановка задачи и примеры

Рассмотрим задачу

A(t)x'(t) + B(t)x(t) = f (t), t € [0,1], (1)

x(0) = xo, (2)

где A(t) и B(t) — (n x п)-матрицы, f (t) — заданная, x(t) — искомая n-мерные вектор-функции, xo —заданный n-мерный вектор.

Системы вида (1) принято называть ДАУ, если

detA(t) =0 Vt € [0,1]. (3)

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

Характеристикой сложности рассматриваемого класса задач является понятие индекса, которое имеет несколько определений. Приведем одно из них (см., например, [3]).

Определение 1. Система (1) имеет решение типа Коши индекса г, если существуют гладкие (n x п)-матрицы $(t), K0(t,T), Ki(t), i = 1, 2,... ,r, причем rank^(t) = k = const Vt € [0,1] такие, что линейная комбинация

ft r

^(t,c) = ^(t)c + Ko(t,T)f(T)dT + J2 Ki(t)f(i-1)(t) (4)

Jo J=i

является решением системы (1) и на любом отрезке [а, в] ^ [0,1] нет решений, отличных от (4). Здесь c € Rn, r < maxte[01] rankA(t).

Приведем примеры, которые иллюстрируют данное определение. Задача

ооН V I Но 1) [V) -( т)' £€ 101]-

имеет единственное решение, которое не зависит от начальных данных: и(£) = 0^(£) = / (£), однако на отрезке [а, в] : а > 0 этот пример имеет решение и(£) = с/£^(£) = 0, где с—любое число.

Система, похожая на предидущую

£ 0 ) /" и \ , ( 01 ) ( и ) ( 0 )

V

00Н V 1 + {10) Ы т)' £ € 101]-

имеет единственное решение, которое не зависит от начальных данных и отрезка интегрирования: v(t) = /(£), и(£) = —£/ (£).

Если коэффициенты матриц А(£) и В(£) являются вещественно-аналитическими функциями, и система (1) имеет общее решение индекса г, то известно (см., например, [2]), что существуют невырожденные для любого £ € [0,1] матрицы Р(£) и Q(t) с вещественно-аналитическими коэффициентами такие, что

Р (£М(£М£) = ( Е №{,.)),

Р(£)А(£)<3(£) + Р(t)B(t)Q(t) = ( J0£) Е0_ к ) , (5)

где Ек и Еп_к — единичные матрицы размером к и п — к соответственно (к то же самое число, что и в определении 1), .](£) — (к х к)-матрица, N (£) = diag{N1, , М3}, N3, ] = 1, 2,... ,в, — верхне-

треугольные матрицы с нулевыми квадратными блоками на диагонали, максимальное число которых равно г — индексу (1).

Существуют и другие определения индекса (1), однако в случае вещественно-аналитических коэффициентов матриц А(£) и В(£) в монографии [5] показана их эквивалентность.

Если ДАУ имеет индекс 2 и выше, то принято называть их ДАУ высокого индекса. Для численного решения таких задач ряд стандартных неявных методов, эффективных для численного решения жестких ОДУ, явялется принципиально неприменим или порождает неустойчивый процесс. Для наглядности иллюстрацию проведем на неявной схеме Эйлера.

Пример 1 (см., напр., [2])

1 а)(и:)+(11 г)(:)=(I)• **ю,ч, ^

где а — числовой параметр.

Опуская элементарные выкладки, легко показать, что данная система имеет единственное решение и = д — а£и, V = / — д и ее индекс равен двум. Применим неявный метод Эйлера для данного примера

А^(Х+ — Хг) + НВг+1Хг+1 = ^¿+1, (7)

гдех = (и, v)т, Аг+1 = А(£г+1), В+ = В(£¿+0,^+1 = Л(/(£¿+1),д(£г+1))т, V0 = /(0) — д/(£)|*=о, ио = д(0), £¿+1 = (г + 1)л, г = 0,1— 1, N = 1/Л,.

Опуская элементарные выкладки,получим, что

Хг+1 = ( 0 а^"+^ Хг + (Аг+1 + ЛВг+1)_^+1,

т.е. при а < —1/2 метод неустойчив, а при а = —1 принципиально не применим из-за сингулярности матрицы (Аг+1 + ЛВг+1).

Пример 2 [6, 7]. Рассмотрим ДАУ

01КV)+(а —)(и)=(дм)■ £^м- <8>

где а и в — скалярные параметры. Обозначим д(£) = а + @£ — 1. Если д(£) = 0 Vt € [а, Ь], то ДАУ (8) имеет решение типа Коши индекса 2 (см. формулу (4))

и(іП _ 1 у(і) ) д(і)

(9)

Обозначим det(A(£) + ЛВ(£)) = Л2«(£), где «(£) = а + в£. Применяя для данного примера неявный метод Эйлера и опуская выкладки, получим

Х+ = ж+Г)( 0)Х^ (10)

где ^>¿+1 некоторые вектор-вункции, конкретный вид которых для нас, в данном случае, не важен.

Если а и в выбраны таким образом, что |s(t)| < 1 Vt € [а, Ь], то получим неустойчивый процесс. С другой стороны, для любых фиксированных параметров а и в можно выбрать отрезок [а, Ь] так, что |s(t)| < 1 V£ € [а, Ь] и д(£) = 0 Vt € [а, Ь]. Таким образом, для данного примера можно по заданным параметрам а и в можно указать отрезок [а, Ь] или по заданному отрезку [а, Ь] указать параметры а и в таким образом, что неявный метод Эйлера будет неустойчивым.

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

метода Эйлера (аналог того, что для численного решения жестких ОДУ применять явный метод Эйлера).

Пример 3 [8].

а—1 ?)(и: И—ва-—1) ——о о- £ € ^

где а = 1 и в — скалярные параметры. Данный пример имеет общее решение

и =1—а V, V = e(a_в)tv(0).

а—1

При а < в и £ ^ то имеем и ^ 0, V ^ 0.

Применяя для данного примера неявный метод Эйлера, получим

1 — а£г+1 1 + аЛ

и+ = ^ = ТТвн"'- (12)

Таким образом, данный метод будет устойчив тогда и только тогда, когда выполнено весьма жесткое ограничение на шаг интегрирования |1 + hа| < |1 + Лв|.

Пример 4 [9].

1 —а£ ) ( и \ = / Л а(1 — А£) \ / и 0 0 ) \ V = I —1 1 + а£ М V

(—Л1 {и)' и<0> = V«»“1. (13>

где а и Л — скалярные параметры. Точное решение данного примера и = (1 + а£)ехь, V = ехь. При А < 0 и £ ^ то имеем и ^ 0 и V ^ 0.

По аналогии с предыдущим примером, неявный метод Эйлера для данного случая дает

(и+1)=(01+а*,+1)(V:), <14)

где Я(г,и) = , и = аЛ, г = АЛ.

Я(г, и) называют функцией устойчивости [9] для ДАУ по аналогии с функцией устойчивости для модельного уравнения Далквиста Х = Ах. Видно, что для любого г < 0 можно подобрать такое и, что |Я(г, и)| > 1, т.е. метод будет неустойчивым. В работе [9] показано, что ряд неявных методов Рунге-Кутта (РК) для данного примера также неэффективен: для устойчивости требуется достаточно малый шаг интегрирования.

Вышеприведенные примеры показывают необходимость разработки эффективных разностных схем для ДАУ.

3. Схемы первого и второго порядков

Одним из подходов к построению эффективных алгоритмов для ДАУ является разработка блочных методов [3, 4], которые включают в себя как частные случаи линейные многошаговые методы (ММ) и некоторые методы РК. Принципиальным отличием таких схем от ММ и методов РК является различная аппроксимация матриц А(£) и Б(£) в узлах сетки. Приведем две таких схемы

Аг+1/2(хг+1 — хг) + 2(В+1Х+1 + Вгхг) — 2(/г+1 + /¿)> (16)

А:(х:+1 — Х:) + ЛБ:+1Х+ = Л/:+1, (15)

ЛЛ

2(Б:+1Х:+1 + Б:Х:) — 2

где А:+1/2 = А(£: + Л/2), х: — аппроксимация х(£:).

При реализации схем (15) и (16) на каждом шаге нам необходимо решать системы алгебраических уравнений (СЛАУ) с матрицами А: + ЛБ:+1 и А:+1/2 + (Л/2)Б:+1 соответственно. Несмотря на то, что матрица А(£) + ЛБ(£) может быть тождественно вырожденной, видно, что матрицы А: + ЛБ:+1 и А:+1/2 + (Л/2)Б:+1 могут быть невырожденными. Опуская выкладки, приведем анализ схемы (15) для примеров 1-4.

Для примера 1 схема (15) дает

, е (5':+1 — 9:)

и:+1 = 9:+1 — аt:+lv:+l, v:+l = /:+1---л-----.

Для примера 2 получим

и:+1 = а + в*1+1 — 1 (/:+1 + в9:+1 — (9:+1Л 9г)) , V:+1 = 9: — £и:+1.

Таким образом, схема (15) для данных примеров сходится к точному решению с первым порядком при любых (за исключением для примера

2 а + в* — 1 = 0) значениях параметров.

Для жестких ДАУ индекса 1, примеры 3 и 4, соответственно, получим

1 а£:+1 1

и:+1 = ^—^ ^ ^ =1 — Л(а — в) ^

и:+1 = (1 + а£:+1 )v:+l, v:+l = ^ 'h\V:.

Сравнивая последние формулы с формулами (12) и (14), замечаем принципиальное отличие и очевидную предпочтительность метода (15) над неявным методом Эйлера.

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

Лемма 1. Если система (1) имеет индекс 2 и матрица Р в формуле (5) является постоянной, то

І

II П (А + НБ^+1)-1А^ ||< К < то, і — 1,2,...,Ж - 1.

¿=о

Лемма 2. Если система (1) имеет индекс 1, то І

II П (А^- + ЛБі+1)-1А^ II < К < то, і — 0,1,...,Ж - 1.

¿=о

Лемма 3. Если система (1) имеет индекс 1, то І

II П (А?+^/2 + ^/2БІ+1) 1 (Аі+^/2 — К/2БІ) II < К < то,і — 0, 1,...,N — 1 ¿=0

Относительно сходимости схемы (15) справедлива

Теорема 1. Пусть задача (1)-(2) имеет индекс 2 и матрица Р(¿) — Р в формуле (5) постоянная. Тогда справедлива оценка

||хг — х(іг )|| — О (К), і — 1,2,..., N.

Доказательство. Доказательство этой теоремы основано на лемме 1 и том, что локальная погрешность аппроксимации имеет второй порядок.

□

Относительно сходимости схемы (16) справедливы следующие утверждения.

Теорема 2. Пусть задача (1)-(2) имеет индекс 1, элементы матриц А(і),Б(і) и вектор-функции f (і) являются непрерывно дифференцируемыми. Тогда справедлива оценка

||хг — х(и)|| — О(К), і — 1,2,..., N.

Теорема 3. Пусть задача (1)-(2) имеет индекс 1 и элементы А(і), Б(¿), f (¿) дважды непрерывно дифференцируемы. Тогда справедлива оценка

||хг — х(^г)I — О(К2), і — 1,2,..., N.

Доказательство данных теорем основано на леммах 2 и 3 соответственно и на оценке локальной погрешности формул Эйлера и трапеций.

Пока не получены конструктивные результаты о сходимости схемы (16) для ДАУ индекса 2. Тем не менее результаты расчетов позволяют надеяться на работоспособность этой схемы и для ДАУ индекса 2.

В заключение этого раздела приведем результаты расчетов модельных примеров по схеме (16), оформленные в виде таблиц.

Для примера 4 при следующих значениях входных данныхА = -20; а = 30, -и(0) = 1 имеем

ь 0.2 0.1 0.05 0.025

еГи 1.27 ■ 10-1 4.26 ■ 10-7 1.83 ■ 10-7 6.63 ■ 10-8

?5 <13 1 1 О 1 со 1.52 ■ 10-8 6.2 ■ 10-9 2.2 ■ 10-9

Здесь Н — шаг сетки, еги = |и(1) — |, ег^ = |-и(1) — зд|.

Приведем пример индекса 2 с сингулярным пучком АА(*) + В(*).

10) (и )Ч00)С) + (—е* — ■ * -10.11.

ж(0) = (1,1)т, который имеет точное решение ж(*) = (е*,е-*)т. Результаты расчетов по схему (16) представлены в таблице

ь 0.2 0.1 0.05 0.025

ег„ - 0 1 .5 2. 4. 5 1 О 1 1.1 ■ 10-4 2. 8 1 О 1 Сл

В этой таблице принято обозначение ег = тах{|и(1) — им|, |V(1) —

Ц.

4. Блочные методы

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

Аг+1 (Жг+1 — Ж*) + Н£*Жг+1 = Н/г+1, (17)

то опуская несложные выкладки легко убедиться, что данная схема, также как и стандартный неявный метод Эйлера, будет неустойчива для первого и второго примеров и неэффективной для примеров

3 и 4 (требуется малый шаг интегрирования) при тех же значениях

параметров. Поясним почему так происходит.

Перепишем исходную систему (1) в виде:

(А(*)ж(*))' + (В(*) — А(*))ж(*) = /(*), * € [0,1], (18)

и для задачи (18), (2) применим неявный метод Эйлера. Имеем

Аг+1Жг+1 — АгЖг + Н(Вг+1 — Аг+1)жг+1 =

(Аг+1 — Аг+1)жг+1 — Агжг + Н£г+1жг+1 — Н/г+1.

Учитывая, что Аг+1 — А^+1 ~ А из последней схемы получим метод (15).

Запись (18) позволяет строить методы типа блочных для исходной задачи (1), (2) высокого порядка точности, которые были исследованы и обоснованыдля задачи (1), (2) в статьях [3, 4, 6] Приведем описание частного случая таких схем для задачи (18), (2), а именно интерполяционный вариант.

На отрезке [0,1] зададим равномерную сетку — гЛ, г — 0,1,..., N. Н — 1/Ж. Предположим, что достаточно точно вычислены стартовые значения ж(Н), ж(2Н),..., ж(тН), которые обозначим как ж1,ж2, ...,жт

Используя точки Жг+з, Жг+8-1, ... жг-т ,г — т, Ш + 1, ..., N — 8 построим интерполяционный полином степени т+в . Подставим данный полином в (18) и потребуем чтобы в точках сетки ¿¿+8, ¿г+з-1,... £¿+1 данное уравнение обращалось в тождество. В этом мы будем иметь 8—стадийные, т—шаговые блочные схемы вида

Е^'=0 ^1Аг+«-.7 жг+'-,7 + НСг+1жг+1 — Н/г+1

Ху7=0 ^7Аг+«-.7жг+з—+ НСг+2жг+2 — Н/г+2

, (19)

„ Е^'=0 кЗАг+«-.7жг+'-,7 + НСг+«жг+« — Н/г+'

где Сч — £(¿4) — А (£д).

Как мы видим при реализации схем (19) необходимо вычислять А (¿) в точках сетки. Предлагается приближенно вычислять значения А (¿) в точках ¿г+', ¿г+'—1,... ¿¿+1 по тем же самым формулам, что и ж (¿). Таким образом мы получим разностные схемы вида

' Е ”+ Л)Аг+'_^жг+'-^ + (Н£г+1 — Ет+ Л)Аг+'-)ж*+1 — Н/г+1

Е”1+о' Л|Аг+'-^жг+'-^- + (Н£г+1 — Е”1+о' Л|Аг+'-^)жг+2 — Н/г+2

,

„ Ет+0' Й'Аг+'-^-жг+'-^- + (Н£г+1 — Е™+' й'Аг+з-)жг+' — Н/г+'

(20)

Если т — в, и матрица А(£) — А постоянная, то получим методы типа РК, конкретный вид которых приведен в [4].

Приведем несколько схем для случая, когда т > 8. Если в — 1, то получим методы, основанные на формуле дифференцирования назад

для задачи (18), (2) вида

т т

'У ^7 Аі+1— іХі+1- і + (НВг+1 ^ Аі+1— і )хі+1 — Н/г+1-

і=о і=0

Если § — 2, т — 3, то получим

Г (2Аі+2хі+2 + 3Аі+1Хі+1 6АіХі + Аі-1 хі-1) + Сі+1хі+1 — 6Н/г+1

\ (1іАі+2Жі+2 - 18Аі+1Жі+1 + 9АіХі - 2Аі-1Жі-1 ) + Сі+2Жі+2 — 6Н/і+2,

где Сі+1 — 6^Вг+1 — 2Аі+2 — 3Аі+1 + 6АІ — Аі-1 Сі+2 — 6^Вг+2 —

1іАі+2 + 18Аі+1 — 9 А + 2Аі-1 и предполагая, что х1 заранее задано, или вычислено в дополнение к начальному значению Хо — а. . Если в — 2, т — 4, то получим

Г А1Аі+2Хі+2 + Сг+1Хг+1 — 12Н/г+1 \ Д2Аі+2Жі+2 + Сі+2Жі+2 — 12Н/і+2,

где разностные операторы А1 Аі+2 и А1Аі+2, действуюшие на матрицу или вектор функцию, пределены по правилу

А1Аі+2Жі+2 — ЗАі+2Жі+2 + 10АтЖі+1 — 18АіХі + 6А*-1Ж*-1 — А*-2Жг-2,

А 2 Аі+2Жі+2 — 25Аі+2Жі+2 — 48Аі+1Жі+1+36А*ж* — 16А*-1Ж*—+ЗАі-2 ж^-2,

а Сі+1 — 12Бі+1 — А1 Аі+2 , Сі+2 — 12Бі+1 — А2^+2. В данной схеме предполагая, что ж1,ж2 заранее задано в дополнение к начальному значению х0 — а.

Детальное обоснование блочных схем (20) предполагается провести в дальнейшем.

Список литературы

1. Хайрер, Э. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи /Э. Хайрер, Г. Ваннер - М.: Мир, 1999. М.: Наука, 1964.

2. Бгепап, К.Е. Numerical solution of initial-value problems in differential-algebraic equations / К.Е. Бгепап, S.B. Campbell, L.R. Petzold - North Holland, New York, 1989.

3. Булатов, M.B. Об условиях сходимости разностных схем для систем ОДУ, не разрешенных относительно производных / М.В. Булатов, В.Ф. Чистяков // Методы численного анализа и оптимизации. - Новосибирск: Наука,- 1987.-С. 175-187.

4. Bulatov, M.V. Numerical solution of differential-algebraic equations by block methods / M.V. Bulatov // Computational Science-ICCS 2003, Springer Verlag, Part 2 P. 516-522.

5. Чистяков, В.Ф. Алгебро-дифференциальные операторы с конечномерным ядром / В.Ф. Чистяков - Новосибирск: Наука, 1996.

6. Булатов, М.В.Об одном численном методе решения дифференциальноалгебраических уравнений /М.В. Булатов, В.Ф. Чистяков // Журн. вы-числ.матем. и мат. физики. -2002. -Т. 42, No 4. - С. 459-470.

7. Данилов, В.А., Чистяков В.Ф. О препятствиях на пути построения эффективных численных методов решения алгебро-дифференциальных систем. Иркутск, 1990, Препринт No 5 ИрВЦ СО АН СССР.

8. Marz R. Differential-algebraic systems anew // Appllied Numer. Math., 42, 2002, pp. 315-335.

9. Kunkel, P., Mehrmann V. Stability properties of differential-algebraic equations and spin-stabilized diskretizations // Electr. Trans. Numer. Anal.., 26, 2007, pp.385-420.

M. V. Bulatov, L. S. Solovarova

On block difference schemes for differential algebraic equations

Abstract. In this paper, block difference schemes for numerical solution of the initial value problems in linear differential algebraic equations are proposed. The detailed analysis of such schemes of the first and second order have been studied using modal examples. It has been shown that the schemes considered possess some advatnages when compared with known implicit first and second order methods. A description for the general form of block difference schemes has been given.

Keywords: differential-algebraic equations, index, block difference schemes

Булатов Михаил Валерьянович, доктор физико-математических наук, главный научный сотрудник, Институт динамики систем и теории управления СО РАН, 664033, Иркутск, ул. Лермонтова, 134, тел.: (3952)453018 ([email protected])

Соловарова Любовь Степановна, аспирант, Институт динамики систем и теории управления СО РАН, 664033, Иркутск, ул. Лермонтова, 134, тел.: (3952)453018 ([email protected])

Bulatov Mikhail, Dr.Sc., Institute of System Dynamics and Control Theory, Irkutsk, Lermontov St., 134 Phone: (3952)453018 ([email protected])

Solovarova Lyubov, PhD Student, Institute of System Dynamics and Control Theory, Irkutsk, Lermontov St., 134 Phone: (3952)453018 ([email protected])