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

ISSN 2074-1863 Уфимский математический журнал. Том 1. № 4 (2009). С. 125-132.

УДК 519.61

ВЕРХНИЕ ОЦЕНКИ НОРМ РЕШЕНИЙ ЭРМИТОВЫХ СИСТЕМ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ И АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ

Ю.М. НЕЧЕПУРЕНКО, Г.В. ОВЧИННИКОВ

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

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

1. Введение

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

dx

x(0) = x0, E— = Ax + f, (1.1)

где f (t) — достаточно гладкая n-компонентная векторная функция, A и E, соответственно,

— отрицательно и неотрицательно определенные эрмитовые матрицы порядка n и пучок AE — A регулярный. Конечные собственные значения такого пучка, т.е. корни уравнения det(AE — A) = 0, являются вещественными отрицательными [1]. Такие системы возникают в различных приложениях, например, при анализе влияния неидеальности межсоединений в интегральных схемах на прохождение сигнала [2].

Обозначим через P правый проектор Рисса на понижающее подпространство, отвечающее множеству конечных собственных значений пучка AE — A [3]:

P = —i(AE — A)-1EdA, (1.2)

2гп Jг

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

Yu.M. Nechepurenko, G.V. Ovchinnikov, Upper bounds for the solution norms of the Hermitian ODAE systems.

© Нечепуренко Ю.М., Овчинников Г.В. 2009.

Работа поддержана РФФИ (грант 07-01-00658-a) и программой РАН "Современные проблемы теоретической математики", проект "Оптимизация вычислительных алгоритмов решения задач математической физики".

Поступила 28 октября 2009 г.

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

Р* = — [ Е(АЕ - А)-1&А,

2т Jг

и для этих проекторов справедливы следующие равенства:

а) Р*Е = Е = ЕР, Ь) Р*А = АР. (1.3)

Решение рассматриваемой задачи Коши существует и единственно в том и только том случае, если

(1 -Р*)(Ах° + /(0)) = 0. (1.4)

Введем в рассмотрение матрицу

Е + = 2^/(АЕ - А)-1АА, (1.5)

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

а) Е+Е = Р, Ь) ЕЕ + = Р*, с) ЕЕ+Е = Е, &) Е+ЕЕ+ = Е+, (1.6)

первые два из которых очевидны, третье следует из первого и (1.3а), а четвертое будет

доказано в разделе 2. Таким образом, Е + совпадает с традиционной псевдообратной матрицей [4], если Р = Р*, т.е. когда Р - ортопроектор.

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

Теорема 1. Пусть х - решение задачи Коши (1.1), (1-4)- Тогда при всех £ > 0 справедливы неравенства

||Рх(£)||2 ^ се-к*||х°||2 + & Ге-к(*-т)||/(т)||2&т (1.7)

11(1 -Р)х(*)||2 ^ ||А-1|Ы|/(*)||2, (1.8)

где

с = V||Е|Ы|Е +1|2, & = ||Е+||2, к = -Атах(А,Е) < 1/(||Е|Ы|А-1||2),

а Атах(А,Е) означает максимальное конечное собственное значение пучка АЕ - А.

Наряду с этими оценками, мы выведем нижеследующую оценку нормы вектора

у(£) = х(£) + А-1/СО, (1.9)

полезную для анализа установления решения.

Теорема 2. Пусть х — решение задачи Коши (1.1), (1-4) и у — вектор (1.9). Тогда при всех £ > 0 справедливо неравенство

||у(()||2 < се-“ { ||у(0)|2 + 1‘ ||А-1/ '(т )||2&т1 , (1.10)

где с и к означают константы, определенные в теореме 1.

Доказательству теорем 1 и 2 посвящен раздел 2 данной работы. В разделе 3 обсуждается численная реализация полученных оценок, а в разделе 4 проводится их сравнение с известной оценкой нормы решения задачи Коши, основанной на уравнении Ляпунова.

°

и

Отметим, что согласно формуле Коши [5] решение задачи Коши (1.1) с единичной матрицей Е представимо в следующем виде:

х(£) = е*Ах° + [ е(*-т)А/(т)&т,

°

откуда можно вывести формулу

у(£) = е“у(0) + /‘ е(‘-т )АА-1 /' (т )*-.

°

Используя эти две формулы, несложно показать, что оценки (1.7) и (1.10) достижимы. Например, неравенство (1.10) превращается в равенство, если Е и, следовательно, Е + — единичные матрицы, х° — собственный вектор матрицы А, отвечающий ее минимальному по модулю собственному значению, а /(£) = х°<^(£), где ^ — скалярная функция такая, что ^(¿) < 0 и ^'(¿) < 0. Оценка (1.8), очевидно, достижима в случае нулевой матрицы Е.

2. Доказательства

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

Лемма 1. Матрица (1.5) удовлетворяет равенствам (1.6). Никакая другая матрица одновременно всем этим равенствам удовлетворять не может.

Доказательство. Поскольку матрица А отрицательно определенная, а Е — неотрицательно определенная, найдется такая невырожденная матрица Р, что

а) Р * АР = Д Ь) Р *ЕР

1Г 0

0 0

(2.1)

где О — отрицательно определенная диагональная матрица, а 1г — единичная матрица порядка г = гапкЕ [1]. Используя эти равенства и учитывая, что первые г диагональных элементов матрицы О образуют множество конечных собственных значений пучка АЕ - А, можно показать, что для проектора (1.2) справедливо представление

1Г 0

00

=Р

Р-1, (2.2)

а для псевдообратной матрицы (1.5) - представление

/г 0

Е + = Р

1Т 00

Р*. (2.3)

Опираясь на (2.1), (2.2) и (2.3), несложно непосредственно проверить справедливость равенств (1.6). Таким образом, мы доказали, что матрица (1.5) удовлетворяет (1.6). Осталось показать, что никакая другая матрица одновременно всем этим равенствам удовлетворять не может.

Очевидно, достаточно убедиться, что система уравнений

а) (Е + + Д)Е = Р, Ь) Е(Е + + А) = Р*, (2.4)

с) (Е + + Д)Е(Е + + А) = Е + + А имеет только тривиальное решение А = 0. Из (2.4а,Ь) и (1.6а,Ь) следует, что ДЕ = ЕА = 0. Раскрывая скобки в (2.4с) и используя эти равенства и равенство (1.6^, получаем требуемое равенство А = 0.

Применяя (2.1), (2.2) и (2.3), нетрудно проверить (1.3) и следующие равенства:

а) РЕ + = Е + = Е+Р*, Ь) РА-1 = А-1Р*, (2.5)

которые нам потребуются в дальнейшем. Отметим, что (2.5а) можно также вывести из (1.6а,Ь^), а (2.5Ь) непосредственно следует из (1.3Ь).

Лемма 2. При всех t > 0 справедливы неравенства

a) ||eie+AP||2 ^ ce-Kt, b) ||etE+AE+1|2 ^ de-Kt, (2.6)

12

где к, с и d — константы, определенные в теореме 1. Доказательство. Из (2.1), (2.2) и (2.3) следует, что

etDr 0 00

etE+Av = р

P-1, etE+AE+ = P

etDr 0 00

р *

где Dr означает главную подматрицу порядка r матрицы D, диагональные элементы которой образуют множество конечных собственных значений пучка AE — A. Поэтому неравенства (2.6) справедливы с константами с = ||Pr||2||(P-1)r||2 и d = ||Pr||2, где Pr и (P-1)r означают, соответственно, n х r-матрицы первых r столбцов матриц P и P-1, и

к = — Amax(A> E) = Amin(A — E) = Amin(1ra, — EA 1) =

= 1/Amax( —EA-1 ,/„) > 1/||E || 2 || A-11| 2.

Но, в силу (2.1b), ||E||2 = ||(P-1)r 112, а в силу (2.3), ||E+||2 = ||Pr||2.

Перейдем к доказательству теорем 1 и 2. Введем следующие обозначения: x1 = Px, x2 = (I — P)x. Умножая (1.1) слева на E + и используя равенства (1.6a), (1.3b) и (2.5a), для x1 можно вывести следующее уравнение:

dXt1 = E+Ai, + E +/.

dt

Поэтому в соответствии с формулой Коши

X1(t) = etE+AP x0 + /t e(t-T)E+AE+/(т )dr. (2.7)

Jo

Неравенство (1.7) непосредственно следует из (2.7) и леммы 2.

Умножая теперь (1.1) слева на (I — P)A-1 и учитывая, что, в силу (2.5b) и (1.3a) (I — P)A-1E = 0, получим следующее равенство:

0 = Х2 + (/ — P )A-1/. (2.8)

Поскольку A — эрмитовая отрицательно определенная матрица и каждая из двух матриц в правой части тождества

A-1 = P A-1P * + (I — P )A-1(1 — P *)

— эрмитовая неположительно определенная, справедливо неравенство

||A-1||2 > ||(I — P)A-1(1 — P*)12 = ||(I — P)A-1||2.

Следовательно,

l|X2(t)||2 = ||(1 — P )A-1/(t) || 2 ^ ||A-1|2|/(t) || 2, что доказывает неравенство (1.8). Теорема 1 доказана.

Отметим, что из равенств (2.7) и (2.8) непосредственно следует, что условие (1.4) является необходимым и достаточным для существования и единственности решения задачи Коши (1.1).

Для доказательства теоремы 2 перепишем (1.1) следующим образом:

y(0) = x0 + A-1/(0), Edt = Ay + EA-1/'. (2.9)

Умножая (2.9) слева на (I — P)A-1, можно установить, что y = Py. Поэтому, умножая

(2.9) слева на E + и используя формулу Коши, получим

y(t) = etE+APy(0) + Г e(t-T)E+APA-1/'(т)dT. (2.10)

Неравенство (1.10) непосредственно следует из (2.10) и (2.6а). Теорема 2 доказана.

3. Численная реализация

Учитывая, что вторая норма квадратной эрмитовой неотрицательно определенной матрицы Е равна ее максимальному собственному значению, а вторая норма обратной матрицы для эрмитовой отрицательно определенной матрицы А является величиной, обратной минимальному собственному значению матрицы - А, вычисление этих норм сводится к вычислению экстремальных собственных значений эрмитовых матриц и является, таким образом, стандартной задачей матричного численного анализа, имеющей эффективные алгоритмические решения как в случае матриц небольшого размера общего вида, так и в случае больших матриц, допускающих быстрое умножение на вектор, например, разреженных. Величину к удобно вычислять, как обратную максимальному собственному значению пучка А(-А) - Е с эрмитовой положительно определенной матрицей - А и эрмитовой неотрицательно определенной матрицей Е. Это задача также эффективно решается известными методами, описанными, например, в [1] и [6] и реализованными в различных численных пакетах. Поэтому поиск эффективной численной реализации оценок из теорем

1 и 2 сводится к поиску методов вычисления ||Е+||2 или ее оценки сверху.

Как уже было отмечено в доказательстве леммы 2, Е + = Рг Рг*, где Рг означает п х г-матрицу первых г столбцов матрицы Р, дающую разложение (2.1). Матрицу Р можно вычислить следующим образом. Выполнить разложение Холецкого: -А = ЬЬ*. Сформировать матрицу Е = Ь-1ЕЬ-* и найти ее спектральное разложение:

Л. 0 00

Е = Я

Я*,

где Лг — диагональная матрица порядка г ненулевых собственных значений матрицы Е, а Я унитарная матрица ее собственных векторов. Тогда

Л-1/2 0

0 г

Р = £-*Я

т.е. требуемую матрицу Рг можно будет вычислить по формуле

рг = ¿-*а л-1/2 ,

где означает п х г-матрицу первых г столбцов матрицы Я.

Возможна также более тонкая реализация описанного выше подхода, значительно более устойчивая к накоплению погрешностей округления, которые могут существенно исказить минимальные ненулевые собственные значения матрицы Е и, как следствие, вычисленную норму матрицы Рг. Вместо матрицы Е вычислим п х г-матрицу У. такую, что Е = У. У.*. Для этого найдем п х г-матрицу Уг, дающую скелетное разложение Е = У.У.* (например, на основе спектрального разложения матрицы Е или разложения Холецкого с перестановками). Тогда У. = ¿-1У.. После этого выполним ^Д-разложение Уг = Ц.Дг, где Ц. — унитарная прямоугольная п х г-матрица, а Дг — квадратная верхняя треугольная порядка г, и найдем

Р. = Ь-*и. Д-1.

Описанный алгоритм вычисления нормы матрицы Е +, основанный фактически на прямом вычислении этой матрицы, можно легко реализовать с помощью стандартных подпрограмм, выполняющих разложение Холецкого, спектральное разложение и т.п., например используя пакет [7]. Единственная проблема, которая при этом может возникнуть, — это численное определение ранга г матрицы Е, если г не известно заранее, а отношение минимального ненулевого собственного значения матрицы Е к ее максимальному собственному значению — величина порядка машинной точности или меньше. Но в последнем случае

сама задача Коши (1.1) становится слишком плохо обусловленной, чтобы ее можно было численно исследовать или решать в условиях приближенных вычислений.

Описанный алгоритм требует, вообще говоря, 0(п3) арифметических операций и 0(п2) ячеек оперативной памяти, что, очевидно, неприемлемо для матриц большого размера. Опишем один из возможных подходов к вычислению | Е+| 2 в случае больших матриц А и Е, допускающих быстрое умножение на вектор, например, сильно разреженных. Рассмотрим матрицу

Е+ = (Е - аА)-1Е(Е - аА)-1.

В силу (2.1)

E+ = P

(Ir - аД.)-2 0

0 0

р *

и, следовательно,

lim E+ = E +.

а—> 0

Причем, как нетрудно видеть, E+ < E+, если а > ß > 1/Amin(A,E).

Таким образом, норму матрицы E + , т.е. ее максимальное собственное значение, можно вычислять как предел максимального собственного значения матрицы E+ при а ^ 0. Для этого, выбрав некоторую монотонно сходящуюся к 0 последовательность а^ > 0, будем вычислять максимальное собственное значение матрицы E+ и отвечающий ему собственный вектор для этих значений а. Для вычисления собственной пары при фиксированном а можно использовать метод Ланцоша или одновременных итераций, или максимизировать отношение Релея методом сопряженных градиентов. При этом саму матрицу E+ формировать не требуется, а ее умножение на вектор сведется к решению системы с матрицей E — аА, умножению на вектор матрицы E и еще одному решению системы с матрицей E — аА. Эти системы можно решать методами, требующими лишь умножения матрицы на вектор, например, методом сопряженных градиентов, либо, если матрица E — аА ленточная с небольшой шириной ленты, на основе разложения Холецкого.

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

После того, как процесс сошелся с заданной точностью при некотором значении а = а* < 1/|Amin(A,E)|, можно найти максимальное собственное значение матрицы E+a. Тогда для искомой нормы будем иметь двустороннюю оценку

IIE+ «2 <||E +»2 <||E±„.||2,

позволяющую проконтролировать точность вычислений.

4. Сравнение с оценкой, основанной на уравнении Ляпунова

Качество оценок (1.7) и (1.10) из теорем 1 и 2, очевидно, определяется оценками нормы матричной экспоненты, полученными в лемме 2.

Пусть М — произвольная квадратная матрица порядка п с собственными значениями, лежащими строго в левой полуплоскости. Тогда, как это показано например в [5], при всех £ > 0 справедливо следующее неравенство:

Не"«2« VII*У*-112ехр(-эдцс-тй;). <41>

где С — произвольная эрмитовая положительно определенная матрица, а X — решение уравнения Ляпунова

ХМ + М *Х = -С.

Сравним оценку (2.6а) с оценкой (4.1), примененной к матрице М = Е+А в случае, когда матрица Е невырожденная. Тогда Р - единичная матрица, Е + = Е-1, а матрица Е удовлетворяет уравнению Ляпунова с М = Е-1 А и С = -2А. Таким образом, в рассмотренном случае эти две оценки совпадают.

Отметим, что если М квадратная матрица простой структуры порядка п и отрицательным спектром, то найдутся такие положительно определенная матрица Е и отрицательно определенная матрица А, что М = Е-1А. Это представление дают матрицы Е = V-*У-1 и А = V-* DV-1, где Д — диагональная матрица собственных значений матрицы М, V — матрица отвечающих этим собственным значениям линейно независимых собственных векторов, а (.)-* = ((.)-1)*. Таким образом, оценки (1.7) и (1.10) применимы к задаче Коши вида

х(0) = х0, — = Мх + /,

с матрицей М простой структуры, имеющей отрицательный спектр.

Рассмотренный случай невырожденной матрицы Е наводит на мысль, что и в общем случае можно оценить норму решения задачи Коши (1.1), основываясь на уравнении Ляпунова. Действительно, покажем например, что оценка (2.6а) справедлива с константами

С = VII*Ы* +«2, к = Ат1п(С,Х)/2 > 1/(2||Х||2НС-1Н2),

где С — произвольная эрмитовая положительно определенная матрица, X означает решение обобщенного уравнения Ляпунова

ХЕ+А + АЕ+Х = -Р *СР, X = Р *ХР, (4.2)

а Х+ — матрица, определенная следующими равенствами:

Х+Х = Р, ХХ+ = Р *, Х+ХХ+ = Х+. (4.3)

Опираясь на разложения (2.1а), (2.2) и (2.3), можно показать, что решение Х уравнения (4.2) существует, единственно и представимо в виде

Х. 0 00

где Х. — решение уравнения Ляпунова

Хг Хг = —С.

с матрицей Сг = РГ*СРГ. Матрица Х +, удовлетворяющая (4.3), существует, единственна и представима в виде

Р

-1

X + = Р

X-1 0 0 0

Р*

Пусть и(£) = е АРи0, где и0 произвольный п-компонентный вектор. Тогда и = Ри и, следовательно,

л

— (Хи,и) = (ХЕ+Аи,и) + (АЕ+Хи,и) =

= -(Си,и) ^ -Атщ(С,Х) (Хи,и).

Поэтому

(Хи,и) ^ е-2к‘(Х«°,и0), £ > 0.

Осталось воспользоваться оценками

НХ+Н-1(и,и) < (Хи,и) < ||Х||2(и,и),

справедливыми для всех и = Ри, где правое неравенство очевидно, а левое доказывает следующая цепочка

(и, и) = (Х+Хи, Х+Хи) = (Х+(Х+)1/2Хи, (Х+)1/2Хи) <

< ||X+||2((X+)1/2Xu, (X+)1/2Xu) = ||X+||2(Xu,X+Xu) = ЦХ+Ц2(Xu,u).

Выбирая различные матрицы C, мы будем получать оценки вида (2.6a) с различными значениями констант c, d и к. Однако, за исключением выбора C = — 2A, приводящего к константам из теоремы 1, возможность эффективной реализации полученных таким образом оценок для больших разреженных матриц представляется проблематичной.

СПИСОК ЛИТЕРАТУРЫ

1. G.H. Golub, C.F. Van Loan Matrix computations. The John Hopkins University Press. London. 1996. 664 p.

2. M. Celik, L. Pileggi, A. Odalasioglu IC interconnect analysis. Kluver Academic Publishes. Norwell. 2002. 310 p.

3. Годунов С.К. Современные аспекты линейной алгебры. Новосибирск: Научная книга. 1997. 390 с.

4. R.A. Horn, C.R. Johnson Matrix analysis. Cambridge Univercity Press. Norwell. 1990. 576 p.

5. Годунов С.К. Обыкновенные дифференциальные уравнения с постоянными коэффициентами. Т.1: Краевые задачи. Новосибирск: Изд-во Новосиб. ун-та. 1994. 264 с.

6. Templates for the solution of algebraic eigenvalue problems. Practical guide. Edited by Zh. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst. SIAM. Philadelphia. 2000. 411 p.

7. E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum,

S. Hammarling, A. McKenney,S. Ostrouchov, D. Sorensen LAPACK Users’ Guide, Third edition. SIAM. Philadelphia. 1999. 407 p.

Юрий Михайлович Нечепуренко,

Институт вычислительной математики РАН, ул. Губкина, 8,

119333, г. Москва, Россия E-mail: yumn@inm.ras .ru

Георгий Викторович Овчинников,

Московский физико-технический институт (Государственный Университет), Институтский пер., д.9,

141700, г. Долгопрудный, Россия E-mail: ovgeorge@yandex.ru