CC BY
18
Серия «Математика»
2011. Т. 4, № 4. С. 2-11
Онлайн-доступ к журналу: http://isu.ru/izvestia
УДК 517.518
О методе матричной прогонки для одного класса
дифференциально-алгебраических
уравнений второго порядка *
М. В. Булатов
Институт динамики и теории управления СО РАН
Н. П. Рахвалов
Восточно-Сибирская государственная академия образования
Та Зуй Фыонг
Институт математики вьетнамской академии наук и технологий
Аннотация. В работе рассмотрены численные методы решения краевой задачи для дифференциально-алгебраических уравнений второго порядка. Выделены условия, при выполнении которых предложенные алгоритмы являются устойчивыми и сходятся к точному решению. Приведены результаты численных расчётов.
Ключевые слова: линейные дифференциально-алгебраические уравнения; матричный полином; краевая задача; метод матричной прогонки.
1. Постановка задачи и вспомогательные сведения
Математические модели, встречающиеся при описании различных механических процессов включают в себя взаимосвязанные обыкновенные дифференциальные уравнения первого и второго порядков и алгебраические связи [6, с. 10-24], [8, е. 140-144, 150-157]. Если все эти уравнения линейные, то их можно объединить и записать в виде системы обыкновенных дифференциальных уравнений
где А(£), £(£), С(£) — (п х п)-матрицы, ж(£) — искомая, f (£) — заданная п-мерная вектор-функции, причём ёе! А(Ь) = 0.
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, проекты № 11-01-93005 Вьет_а, № 11-01-00639_а
A(t)x"(t) + B(t)x'(t) + C(t)x(t) = f (t), t e [0,1], (1.1)
Системы вида (1.1), с вырожденной в области определения матрицей A(t), принято называть дифференциально-алгебраическими уравнениями (ДАУ) второго порядка.
Для таких систем, как правило, задают начальные или краевые условия. К настоящему времени для численного решения таких задач применяют следующий подход: вводя обозначение z(t) = (x/T(t),xT(t))T систему (1.1) сводят к системе первого порядка
(A(t) B(t)) z (t)+( -e C(t)) z(t) = (0(t))
и решают ее методами, разработанными для численного решения дифференциально-алгебраических уравнений первого порядка. В данной работе рассмотрена система (1.1) с краевыми условиями
x(0) = x0, x(1) = xN, (1.2)
в предположении, что (1.2) согласованы с правой частью (1.1).
Для численного решения задачи (1.1), (1.2) предлагается трехточечная аппроксимация, которая несколько отличается от стандартной [4, с. 138-139]. В итоге мы будем иметь систему линейных алгебраических уравнений (СЛАУ) с блочно-трехдиагональной матрицей, которую будем решать методом матричной прогонки [4, с. 522-523]. Будем предполагать, что на отрезке интегрирования решение задачи (1.1), (1.2) существует и единственно.
В дальнейшем изложении нам потребуется ряд определений и вспомогательных утверждений [1], [2], [7].
Определение 1. Матричный пучок AA(t) + B(t), где A — скалярный параметр, удовлетворяет критерию «ранг-степень» (имеет индекс
1), если
rankA(t) = deg(det(AA(t) + B(t))) = k = const Vt e [0, 1], (1.3)
где deg() — символ степени многочлена.
Определение 2. Матричный полином AA(t)+^B(t) + C(t), где А и ц - скалярные параметры, имеет простую структуру, если выполнены условия:
rankA(t) = k = const Vt e [0, 1]; (1.4)
rank(A(t)|B(t)) = k + l = const Vt e [0, 1]; (1.5)
det(AA(t) + yB(t) + C (t)) = a0(t)Ak Ц + ..., (1.6)
причем a0(t) =0 Vt e [0, 1].
Если матричный пучок АА(£) + В(£) удовлетворяет критерию «ранг-степень», то существуют [7] невырожденные для любого £ Є [0, 1] матрицы Р(£) и ^(£) такие, что умножая слева исходную систему на матрицу Р (£) и заменяя переменную х(£) = ^(£)у(£), (1-1) приводится к виду
Р (£)А(£)(ф(£)у(£))" + Р (£)£(£)(д(%(£))' + Р (£)С (£)£(% =
(Р (£)А(£)ф(£))у"(£) + (2Р (£)А(£)ф' (£) + Р (£)£(£)д(£))у"(£) +
(Р (£)Л(£)д"(£) + Р (£)В(£)^' (£) + Р (£)С (£)^(£))у(£) = Р (£)/(£), (1-7)
где
Р(£)А(£)^(£) = ( ^ 0 ) ,
2Р(£)А(£)0'(£) + Р(£)В(£)^(£) = ( ^(£) Еп_к ) ,
Р(£)Л(£)д"(£) + Р(£)В(£)$'(£) + Р(£)С(£)$(£) = ( ^(£) ЗД) ) ,
здесь и всюду в дальнейшем Ет единичные матрицы размерности т, 7(£) и Сі(£), і = 1, 2, 3, 4 матричные блоки соответствующего размера.
Таким образом, исходную систему (1-1) путём неособенных преобразований можно привести к виду
(? I) >"+("и •'*( §ЯЙ !Ю •="«»■ »■*■
Если матричный полином АА(£) + ^В(£) +С(£) имеет простую структуру, то известно [1], [2], что существуют невырожденные для любого £ Є [0, 1] матрицы Р(£) и ф(£) такие что, умножением слева на Р(£) и заменой х(£) = <Э(£)у(£), система (1-1) приводится к виду (1.7) с матрицами
/ Е 0 0 \
Р(£)А(£)^(£) И 0 0 0 , (1.9)
\ 0 0 0/
/ 7і(£) 0 ^(£) \
2Р(£)А(£)д'(£) + Р(£)В(£)£(£) И 0 Е 0 , (1.10)
\ 0 0 0 /
Р (£)А(£)^"(£) + Р (£)£(£)д' (£ ) + Р (£ )С (£)^(£ ) =
/ Сі(£) С2(£) 0 \
Сз(£) С4(£) 0 , (1.11)
\ 0 0 Еп_к_1 /
здесь 71(£),72(£) и С^(£),г = 1, 2, 3, 4 матричные блоки соответствующего размера.
Таким образом, исходную систему путём неособенных преобразований можно привести к виду
Ек 0 0 0 0 0 | у" + 0 0 0
7і(£) 0 72(£)
0 Еі 0 І у' + 0 0 0
Сі(£) С2(£) 0
Сз(£) С4(£) 0 0 0 Еп_к_і
У = Р (£)/(£).
(1.12)
2. Численное решение
Опишем численные алгоритмы решения задачи (1.1), (1.2). Зададим на отрезке [0,1] равномерную сетку ^ = гН, і = 1, 2,...,Ж, Н = 1/Ж и применим для численного решения задачи (1.1), (1.2) разностную схему вида
АіДіЖі+і + ЛДД2Жг+1 + Н2СіДзЖі+1 = Н2/і, г = 1,2,...,Ж - 1,
жо = ж(0), зд = ж(1),
где Аі, Ві, Сі и /і, матрицы А(£), В(£), С(£) и свободный член /(£),
вычисленные в точке £і Є [£і_і , £і+і], а разностные операторы Ді^і+і, Д2^і+і, Дз^і+і, где $(£) - некоторая функция, определены по правилам:
Ді$і+і = 9і+і — 29і + 5^_ъ
Д2$і+і = Ро $і+і + Рі#і + Р25'і_і,
Дз^і+і = ^о^і+і + Оі£і + 02£і_Ь
Конкретный выбор точки £і и разностных коэффициентов р^-, о, j = 0,1, 2 будет указан ниже. Такая аппроксимация приводит к системе линейных алгебраических уравнений (СЛАУ) блочно-трехдиагонального вида
Еіжі_і + ^іхі + ^^іхі+і — Я (2.1)
с матрицами
Еі = Аі + Р2НВі + 02Н2Сі,
^і = — 2Аі + ріНВі + о і Н2 С?і , (2.2)
Мі = Аі + ро нВі + °о ^
6 М. В. БУЛАТОВ, Н. П. РАХВАЛОВ, ТА ЗУЙ ФЫОНГ
правой частью ^ = Н2/і, г = 1, 2,...,Ж — 1 и
жо = ж(0), жм = ж(1).
Система (2.1) может быть решена методом матричной прогонки [4, стр. 522]. Суть метода в том, что исходная система (2.1) приводится к блочно-двудиагональному виду:
Е 0 0 . 0
0 Е — а2 . . 0
0 0 Е . 0
0 0 0 . Е
\ 0 0 0
\ / жо
жі Ж2
Е —ам 0Е
\
жм _і
) \ЖМ )
( ві
в2
вз
\
вм
V вм+і )
где матрицы а^+1 и вектора вг+1 определены по рекуррентным соотношениям
ОД+1 = —+ Дга^) 1 Мг, 2 = 1, 2,...,Ж — 1, «1 =0,
А+1 = (Дгаг + ^г) 1 № — )) ^1 = х(0) 2 = 1) 2,...,^ — 1 (2-3)
Применяя обратный ход метода Гаусса, получим
X = о?+1 ж,-+1 + в'+ъ вм+1 = х(1), ] = N — 1,Ж — 2,...,1. (2.4)
Приведём известные определения [5].
Определение 3. Алгоритм метода матричной прогонки корректен, если матрицы (Дгаг + £г) обратимы для 1 < г < N — 1.
Определение 4. Алгоритм метода матричной прогонки называется устойчивым, если ||аг|| < 1 для г = 1, 2,...^ — 1.
Здесь и всюду далее норму матрицы будем определять по формуле
МАИ = тах ( тах |агк|).
1<г<т 1<к<П
Как правило, для численного решения обыкновенных дифференциальных уравнений применяют аппроксимацию в точке £г. В этом случае д2# = 2(#г+1 — £г-1), дз# = £г, Ж£) = Л2/(^), а матрицы ^ из формулы (2.1) имеют вид = —2А(£г) + Л,2С(£г). Данные матрицы будут невырожденными только в случае регулярности пучка АА(£) + С(£) [3], т.е.
^в£(АА(£) + С(£)) = 0 VIА| < Ао V^ € [0, 1].
В данной работе предлагается другая аппроксимация, которая сохраняет второй порядок аппроксимации, но охватывает более широкий
класс задач, а именно данная аппроксимация применима в случае сингулярности матричного пучка АА(£) + С(£). Мы будем выбирать точки = £г-1 или = ^г+1.
Рассмотрим аппроксимацию в точке £г-1. В этом случае разностные коэффициенты р0 = — 1, р1 =2, р2 = — 3, которые соответствуют формуле дифференцирования вперёд, а коэффициенты ст0, ст1, ст2 будем задавать следующим образом ст0 = —1, ст1 =2, ст2 = 0. Такой выбор можно интерпретировать как экстраполяцию функции $(£) в точке £г_1 по заданным значениям функции в точках ^ и £г+ь В этом
случае матрицы Дг, Ьг, Мг (смотри формулы (2.1)- (2.2)) имеют вид:
3
Дг = А(£г-1) — 2 ^Б(£г-1), (2.5)
Ьг = —2А(£г-1) + 2ЛВ(£г_ 1) + 2Л,2С (^г— 1), (2.6)
Мг = А(^г-1) — ^ ^В(^г-1) — Л-2С (/г-^ (2.7)
^ = й2/(£г-1), г = 1, 2,...,N — 1. (2.8)
Если аппроксимировать задачу (1.1), (1.2) в точке £г+1, то по аналогии получим разностные коэффициенты р0 = 2, р1 = —2, р2 = 1, оо = 0, ст1 = 2, ст2 = —1. В этом случае
Дг = А(^г+1) + ^ йВ(^г+1) — Л-2С (^г+1), (2.9)
Ьг = —2А(^г+1) — 2йВ(^г+1) + 2^2с (/г+1^ (2.10)
3
Мг = А(^г+1) + ^ йВ(^г+1), (2.11)
^ = й2/(^г+1), г = 1, 2,...,N — 1. (2.12)
Таким образом, при выполнении условия (1.3) матрицы Ьг будут невырожденны и формально мы можем применить метод матричной прогонки.
Утверждение 1. Пусть:
1) Элементы матриц А(£), В(£), С(£) и вектор-функция /(£) уравнения (1.1) принадлежат классу С[0 1],
2) Матричный пучок АА(£) + В(£) удовлетворяет критерию «ранг-степень», или матричный полином АА(£) + рВ(£) + С(£) имеет простую структуру.
Тогда методы матричной прогонки (2.3), (2.4) определённые по формулам (2.5)-(2.8), или (2.9)-(2.12) корректны, устойчивы и справедлива оценка ||ж./ — ж(^)| = 0(й2), ] = 1, 2,...,N — 1.
Доказательство основано на блочном представлении (1.8) и (1.12). Оно весьма громоздко и поэтому не приводится.
3. Численные эксперименты
Пример 1. Рассмотрим задачу
1 і 0 0
ж''(і) +
0 0 1 2
ж' (Ь) +
0 0 1 і
ж(і) =
2 + 2Ь
Ьз + Ь2 + 6і
с краевыми условиями ж(0) = (0, 0)т, ж(1) = (1,1)т. Матричный пучок АА(£)+В(£) удовлетворяет критерию «ранг-степень», этот пример имеет единственное решение ж(£) = (^^2)т.
Результаты расчётов примера приведены в табл. 1.
Погрешность методов
Таблица 1.
ь 0.1 0.05 0.025 0.0125 0.00625
вГі 0.01630 0.00575 0.00176 0.00049 0.00013
ЄГ2 0.00136 0.00371 0.00097 0.00025 0.00004
Здесь и всюду далее ег, = тах ||ж(Ьі) — ж^|, N = 4, j = 1, 2 погреш-
і<і<м
ность методов (2.1), у которых матрицы Еі, Ьі, Мі вычислены в точках Ьі_ і и Ьі+і, соответственно.
Пример 2. Рассмотрим задачу
А(Ь)у''(Ь) + в(ь)у'(ь) + с(Ь)У(Ь) = Т(ь)
с матрицами
А(Ь) = Р (ь)А(ь)^(ь),
А(ь) = 2Р (ь)А(ь)ф' (Ь) + Р (ь)в(ь)^(ь),
с (ь) = р (*)А(*)д"(*) + р (ь)в (і)д' (ь) + р (ь)с (*)$(*),
где Р (Ь) =
А(Ь) =
0 2 е_*
—і 2(ь + 1 — 6е*)е
і
2
— і (Ь + 1)е_21 е_
В (ь) =
_і
, Ф(і) =
и С (і) =
тЬ 0 0
I е_
_21 е 21 2 2е
(3.1)
_1
і„ _|
краевые
условия определены следующим образом
"(0) = ^-1(0)ж(0), у(1) = д-1(0)ж(1), (3.2)
где ж(0) = (1,0,0)т, ж(1) = (2,1,1)т.
Если систему (3.1) умножить на матрицу Р-1(*) и произвести замену переменной у(*) = ^-1(*)ж(*), то получим систему
/ 10 0 \ / 00 0 \ / 000 \ / 0 \
0 0 0 ж''(*) + 010 ж'(*) ^ 0 0 0 ж(*) = 2*
\0 0 ^ \0 0 ^ \0 0 ^ \ *3 /
с краевыми условиями ж(0) = (1, 0, 0)т, ж(1) = (2,1,1)т, которая имеет единственное решение ж(*) = (* + 1 ,^2, ^3)т. Конкретный вид матриц А(*), В(*), С(*) и векторов у0, не приводится в виду их громоздкости. Результаты расчётов задачи (3.1), (3.2) приведены в табл. 2.
Таблица 2.
Погрешность методов
ь 0.1 0.05 0.025 0.0125 0.00625
&Г\ 0.07051 0.02005 0.00541 0.00141 0.00036
еГ2 0.08614 0.02319 0.00600 0.00153 0.00038
Вышеприведённые результаты расчётов хорошо согласуются с теоретическими выкладками.
Приведём результаты расчётов, когда матричный полином АА(*) + рВ(*) + С(*) не обладает простой структурой.
Пример 3. Рассмотрим задачу
1 * А гг/.\ , / 0 а + 1 \ /, •. , / 0 0 \ , ч ( (2 + * + а)е* \
0^у«+^1 * у"(()Ч»0у(г)Ч (2+()е' ;
с краевыми условиями у(0) = (1,1)т, у(1) = (е, е)т, которая имеет единственное решение у(*) = (е*,е*)т.
Здесь гапк(А) = 1, топк(А|В) = 2, а ^е*(АА(*) + рВ(*) + С(*)) =
ае*( Р р(1 + а)) ^ = А—р2(1+а), т.е. а0(*) = 0. Таким образом условие
(1.6) нарушено.
В табл. 3 приведены результаты численных экспериментов с параметром а = 10.
Таблица 3.
Погрешность методов
ь 0.1 0.05 0.025 0.0125 0.00625
ег 1 0.0206 0.0110 0.0060 0.0033 0.0017
еГ2 0.0201 0.0124 0.0069 0.0036 0.0018
Из табл. 3 видно, что предлагаемые алгоритмы для данного примера имеют первый порядок. Отметим, что при некоторых значениях а наблюдалась неустойчивость алгоритмов (2.3), (2.4).
Пример 4.
/ 0 10 \ / 0 10 \ / 10 0 \ / 0 \
0 0 1 ж"(£) +001 ж'(*) + 010 ж(£) = 0 .
\0 0 ^ \0 0 ^ \0 0 ^ \ ф(£) )
Система имеет единственное решение, которое не зависит от краевых условий ж = (ф/у(£)+2фт(£)+ф"(£), -ф"(£)-ф'(*), ф(*)). Методы (2.3), (2.4) определённые по формулам (2.5)-(2.8) или (2.9)-(2.12) порождают неустойчивый процесс для данного примера.
Таблица 4.
Погрешность методов
ь 0.1 0.05 0.025 0.0125 0.00625
ег1 127.94 553.52 2303.9 3079.15 9004.1
ег 2 24.20 94.36 384.44 1564.5 6324.5
Итак, при нарушении второго условия утверждения нельзя гарантировать сходимость предлагаемых методав к точному решению.
Список литературы
1. Булатов М. В. Об одном семействе матричных троек / М. В. Булатов // Ляпуновские чтения и презентация информационных технологий : материалы конф. - Иркутск, 2002. - С.10.
2. Булатов М. В. Применение матричных полиномов к исследованию линейных дифференциально-алгебраических уравнений высокого порядка / М. В. Булатов, Минг-Гонг Ли // Дифференциальные уравнения. - 2008. - Т.44, № 10. -С. 1299-1306.
3. Гантмахер Ф. Р. Теория матриц / Ф. Р. Гантмахер. - М. : Наука, 1967. - 576 с.
4. Самарский А. А. Теория разностных схем / А. А. Самарский.- М. : Наука, 1989. - 616 с.
5. Самарский А. А. Методы решения сеточных уравнений / А. А. Самарский, Е. С. Николаев. - М. : Наука, 1978. - 590 с.
6. Хайрер Э. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи / Э. Хайрер, Г. Ваннер. - М. : Мир, 1999. - 685 с.
7. Чистяков В.Ф. Алгебро-дифференциальные операторы с конечномерным ядром / В. Ф. Чистяков. - Новосибирск : Наука. Сиб. издат. фирма РАН, 1996. - 278 с.
8. Brenan K. E. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations / К. Е. Brenan, S. L. Campbell, L. R. Petzold. - SIAM. Philadelphia, 1996.
M. V. Bulatov, N. P. Rakhvalov, Ta Duy Phuong Numerical methods of solution of boundary-value problem for differential-algebraic equations of the second order
Abstract. In this paper the numerical methods of solution of boundary-value problem for differential-algebraic equations of the second order are considered. We found conditions fulfillment of which ensures stability and convergence to exact solution of proposed algorithms. The results of numerical calculations are given.
Keywords: linear differential-algebraic equations, matrix polynomial, boundary-va-lue problem, matrix sweep method
Булатов Михаил Валерьянович, доктор физико-математических наук, главный научный сотрудник ИДСТУ СО РАН, 664033, г. Иркутск, ул. Лермонтова, д. 134 (mvbul@icc.ru)
Рахвалов Николай Петрович, аспирант, Восточно-Сибирская государственная академия образования, 664011, г. Иркутск, ул. Нижняя Набережная, д. 6 (nikolar@mail.ru)
Та Зуй Фыонг, главный научный сотрудник, Институт математики вьетнамской академии наук и технологий, Ханой, Вьетнам (tdphuong@math.ac.vn)
Bulatov Mikhail, Chif reseacher, ISDCT SB RAS Irkutsk, Lermontov st. 134 (mvbul@icc.ru )
Rakhvalov Nikolay Ph.D. student Irkutsk State Pedagogical Academy Irkutsk, Lower Quay st., 6 (nikolar@mail.ru )
Ta Duy Phuong, Assoc.Prof. Dr.,Institute of Mathematics, Vietnam Academy of Science and Technology 18 Hoang Quoc Viet road, CauGiay district, 10307, Hanoi, Vietnam (tdphuong@math.ac.vn)
