1611688890-f641c9ec8276824e4686da772eb56520 (826652), страница 64
Текст из файла (страница 64)
если матрица G положительно определена. Из этого неравенства следуеттакже существование такой константы µ > 0, что hGx, xi > µ hx, xi.Неравенство G ⊲ H будем понимать как hGx, xi > hHx, xi для всех x,что равносильно также G − H ⊲ 0.Достаточное условие сходимости итерационного процесса в формеСамарского (3.129) даётТеорема 3.12.1 (теорема Самарского) Если A — симметричная положительно определённая матрица, τ > 0 и B ⊲ 21 τ A, то стационарный итерационный процессBx(k+1) − x(k)+ Ax(k) = b,τk = 0, 1, 2, . .
. ,сходится к решению системы уравнений Ax = b из любого начальногоприближения.Доказательство. Пусть x∗ — решение системы уравнений Ax = b, такчтоx∗ − x∗+ Ax∗ = b.BτЕсли обозначить через z (k) = x(k) − x∗ — погрешность k-го приближения, то она удовлетворяет однородному соотношениюBz (k+1) − z (k)+ Az (k) = 0,τk = 0, 1, 2, . . . .(3.131)Исследуем поведение энергетической нормы погрешности. Покажем сначала, что в условиях теоремы числовая последовательность kz (n) kA =hAz (n) , z (n) i является невозрастающей.3923.
Численные методы линейной алгебрыИз соотношения (3.131) следуетz (k+1) = I − τ B −1 A z (k) ,и(3.132)Az (k+1) = A − τ AB −1 A z (k) .Таким образом,hAz (k+1) , z (k+1) i = hAz (k) , z (k) i − τ hAB −1 Az (k) , z (k) i− τ hAz (k) , B −1 Az (k) i + τ 2 hAB −1 Az (k) , AB −1 Az (k) i.Коль скоро матрица A симметрична,hAB −1 Az (k) , z (k) i = hAz (k) , B −1 Az (k) i,и потомуhAz (k+1) , z (k+1) i =hAz (k) , z (k) i − 2τB−12τ A B −1 Az (k) , B −1 Az (k) . (3.133)Учитывая неравенство B ⊲ 12 τ A, можем заключить, что вычитаемое вправой части полученного равенства всегда неотрицательно. По этойпричинеkz (k+1) kA ≤ kz (k) kA ,так что последовательность kz (k) kA монотонно не возрастает и ограничена снизу нулём.
В силу теоремы Вейерштрасса она имеет предел приk → ∞.Неравенство B ⊲ 21 τ A, т. е. положительная определённость матрицы(B − 12 τ A), означает существование такого η > 0, что для любых y ∈ RnB − 12 τ A y, y ≥ η hy, yi = η kyk22 .Окончательно получаем из (3.133) (k+1) 2 (k) 2z − z + 2ητ B −1 Az (k) 2 ≤ 0AA2для всех k = 0, 1, 2, . . . . Переходя в этом неравенстве к пределу поk → ∞, заключаем, что тогда kB −1 Az (k) k2 → 0. При неособенной матрице B −1 A это возможно лишь при z (k) → 0. Итак, вне зависимости от3.13. Вычисление определителей и обратных матриц393выбора начального приближения итерационный процесс в самом делесходится.Отметим, что из теоремы Самарского следует теорема ОстровскогоРайха (Теорема 3.9.5) о сходимости метода релаксации для СЛАУ ссимметричными положительно определёнными матрицами, а также,как её частный случай, Теорема 3.9.4 о сходимости метода ГауссаЗейделя.
В самом деле, пусть A = L̃ + D + Ũ в обозначениях §3.9д,т. е. L̃ и Ũ — строго нижняя и строго верхняя треугольные части матрицы A, а D — её диагональная часть. Если A симметрична, то L̃ = Ũ ⊤ ,и поэтомуhAx, xi = hL̃x, xi + hDx, xi + hŨ x, xi = hDx, xi + 2hL̃x, xi.ТогдаhBx, xi −12ωhAx, xi = h(D + ω L̃)x, xi − 12 ω hDx, xi + 2hL̃x, xi= 1 − 12 ω hDx, xi > 0при 0 < ω < 2.Дальнейшие результаты в этом направлении читатель может увидеть, к примеру, в [37, 81].3.13Вычисление определителей матрици обратных матрицПредположим, что для матрицы A выполняется LU-разложение.Как отмечалось, выполняемые в представленной нами версии метода Гаусса преобразования — линейное комбинирование строк — не изменяют величины определителя матрицы.
Следовательно, det A равенопределителю получающейся в итоге верхней треугольной матрицы U ,т. е. det A есть произведение диагональных элементов U .Другая возможная трактовка этого результата состоит в том, чтоесли A = LU — треугольное разложение матрицы A, то, как известноиз линейной алгебры,det A = det L · det U.Определитель нижней треугольной матрицы L равен 1, коль скоро наеё диагонали стоят все единицы. Следовательно, как и ранее, det A =3943. Численные методы линейной алгебрыdet U , а точнее — произведению всех диагональных элементов в верхнейтреугольной матрице U .Совершенно аналогичные выводы можно сделать и при использовании других матричных разложений.
Например, если нам удалось получить A = QR — разложение исходной матрицы в произведение ортогональной и правой треугольной, то, коль скоро det Q = 1, искомыйопределитель det A = det R и вычисляется по R как произведение еёдиагональных элементов.Рассмотрим теперь вычисление матрицы, обратной к данной матрице. Отметим, прежде всего, что в современных вычислительных технологиях это приходится делать не слишком часто. Один из примеров, когда подобное вычисление необходимо по существу, — нахождение дифференциала операции обращения матрицы A 7→ A−1 , равногоd(A−1 ) = −A−1 (dA) A−1 .(см., к примеру, [14]). Тогда коэффициенты чувствительности решения системы уравнений Ax = b по отношению к элементам матрицы иправой части (т.
е. производные решения по коэффициентам и правымчастям системы, см. §1.3) даются формулами∂xν= −zνi xj ,∂aij∂xν= zνi ,∂biν = 1, 2, . . . , n, где Z = (zij ) = A−1 — обратная к матрице A.Гораздо чаще встречается необходимость вычисления произведенияобратной матрицы A−1 на какой-то вектор b, и это произведение всегдаследует находить как решение системы уравнений Ax = b какими-либоиз методов для решения СЛАУ. Такой способ заведомо лучше, чем вычисление A−1 b через нахождение обратной A−1 , как по точности, таки по трудоёмкости.Матрица A−1 , обратная к данной матрице A, является решениемматричного уравненияAX = I.Но это уравнение распадается на n уравнений относительно векторных неизвестных, соответствующих отдельным столбцам неизвестнойматрицы X, и потому мы можем решать получающиеся уравнения порознь.3.13.
Вычисление определителей и обратных матриц395Из сказанного следует способ нахождения обратной матрицы: нужно решить n штук систем линейных уравненийAx = e(j) ,j = 1, 2, . . . , n,(3.134)где e(j) — j-ый столбец единичной матрицы I. Это можно сделать, кпримеру, любым из рассмотренных нами выше методов, причём прямые методы здесь особенно удобны в своей матричной трактовке. В самом деле, сначала мы можем выполнить один раз LU-разложение (илиQR-раложение) исходной матрицы A, а затем хранить его и использовать посредством схемы (3.65) (или (3.82)) для различных правыхчастей уравнений (3.134). Если матрица A — симмметричная положительно определённая, то очень удобным может быть разложение Холесского и последующее решение систем уравнений (3.134) с помощьюпредставления (3.73).В прямых методах решения СЛАУ прямой ход, т.
е. приведение исходной системы к треугольному виду, является наиболее трудоёмкойчастью всего алгоритма, которая требует обычно O(n3 ) арифметических операций. Обратный ход (обратная подстановка) — существенноболее лёгкая часть алгоритма, требующая всего O(n2 ) операций. Поэтой причине изложенный выше рецепт однократного LU-разложенияматрицы (или других разложений) позволяет сохранить общую трудоёмкость O(n3 ) для алгоритма вычисления обратной матрицы.Другой подход к обращению матриц — конструирование чисто матричных процедур, не опирающихся на методы решения систем линейных уравнений с векторными неизвестными.
Известен итерационныйметод Шульца для обращения матриц: задавшись специальным начальным приближением X (0) , выполняют итерацииX (k+1) ← X (k) 2I − AX (k) ,k = 0, 1, 2, . . . .(3.135)Метод Шульца — это не что иное как метод Ньютона для решения системы уравнений, применённый к X −1 −A = 0 (см. §4.5б).28 Его можнотакже рассматривать как матричную версию известной процедуры длявычисления обратной величины (см. [12], глава 3).28 Иногда этот метод называют также методом Хотеллинга, так как одновременно с Г.
Шульцем [96] его рассматривал американский экономист и статистикГ. Хотеллинг [91]. Кроме того, встречается (хотя и крайне редко) также названиеметод Бодевига.3963. Численные методы линейной алгебрыПредложение 3.13.1 Метод Шульца сходится тогда и только тогда, когда его начальное приближение X (0) удовлетворяет условиюρ(I − AX (0) ) < 1.Доказательство.
Расчётную формулу метода Шульца можно переписать в видеX (k+1) = 2X (k) − X (k) AX (k) .Умножим обе части этого равенства слева на (−A) и добавим к ним поединичной матрице I, получимI − AX (k+1) = I − 2AX (k) + AX (k) AX (k) ,что равносильноI − AX (k+1) = (I − AX (k) )2 ,k = 0, 1, 2, . . . .Отсюда, в частности, следует, чтоI − AX (k) = I − AX (0)2k,k = 0, 1, 2, . . . .2kЕсли X (k) → A−1 при k → ∞, то I − AX (0)→ 0 — последовательность степеней матрицы сходится к нулю. Тогда необходимоρ(I − AX (0) ) < 1 в силу Предложения 3.3.10.2kИ наоборот, если ρ(I − AX (0) ) < 1, то I − AX (0)→ 0 при k → ∞,и потому должна иметь место сходимость X (k) → A−1 .Из доказательства предложения следует, что метод Шульца имеетквадратичную сходимость.3.14Оценка погрешностиприближённого решенияВ этом параграфе мы рассмотрим практически важный вопрос обоценке погрешности приближённого решения систем линейных алгебраических уравнений.
Первый способ носит общий характер и можетприменяться в любых ситуациях, в частности, не обязательно в связис итерационными методами.3.14. Оценка погрешности приближённого решения397Пусть x̃ — приближённое решение системы уравнений Ax = b, тогдакак x∗ — её точное решение. Тогда, принимая во внимание, что I =A−1 A и Ax∗ = b,kx̃ − x∗ k = A−1Ax̃ − A−1Ax∗ = A−1 (Ax̃ − Ax∗ )≤ A−1 Ax̃ − b,(3.136)где матричная и векторная нормы, естественно, должны быть согласованы.
Величина (Ax̃ − b) — это невязка приближённого решения x̃,которую мы обычно можем вычислять непосредственно по x̃. Как следствие, погрешность решения можно узнать, найдя каким-либо образомили оценив сверху норму обратной матрицы kA−1 k.Иногда из практики можно получать какую-то информацию о значении kA−1 k. Например, если A — симметричная положительно определённая матрица и известна нижняя граница её спектра µ > 0, то изПредложения 3.2.3 следует, чтоkA−1 k2 ≤ µ−1 .Напомним, что аналогичную информацию о спектре матрицы СЛАУмы использовали при оптимизации скалярного предобуславливателя в§3.9г.
Такова ситуация с численным решением некоторых популярныхуравнений математической физики (уравнением Лапласа и его обобщениями, к примеру), для которых дискретные аналоги соответствующихдифференциальных операторов хорошо изучены и известны оценки ихсобственных значений.В общем случае нахождение kA−1 k или хотя бы разумных оценокдля kA−1 k в какой-то норме, которое было бы менее трудоёмким, чемрешение исходной СЛАУ, является нетривиальным делом. Краткий обзор существующих численных процедур для этой цели, которые называются «оценщиками обусловленности», а также дальнейшие ссылкина литературу можно найти в книге [13], §2.4.3.Для конкретных численных методов оценка погрешности приближённого решения иногда может быть выведена из свойств этих методов. Например, в стационарных одношаговых итерационных методах последовательность погрешностей приближений своими свойствами очень близка к геометрической прогрессии, и этим обстоятельствомможно с успехом воспользоваться.3983.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
