В.Б. Андреев ЧИСЛЕННЫЕ МЕТОДЫ Часть I
2
Глава I Вычислительные методы линейной алгебры3 5 С вычислительной точки зрения в линейной алгебре имеются, если понимать их достаточно широко, две основные задачи: 1 ? решение систем линейных уравнений, 2 ? вычисление собственных значений и собственных векторов матрицы. Основное внимание в лекциях будет уделено решению первой задачи, да и то при весьма ограничительных предположениях. Вторая задача более трудная ее мы коснемся менее подробно. В силу теоремы Кронекера - Капелли система линейных алгебраических уравнений Ax = b разрешима тогда и только тогда, когда ранг матрицы A равен рангу расширенной матрицы [Ab]. Это заведомо так, если матрица A квадратная и невырожденная, т.е. det A = 0 . В этом случае система не только разрешима при любых b, но и имеет единственное решение. (Разрешима однозначно). Именно этот случай мы и будем изучать. Методы решения систем линейных алгебраических уравнений делятся на две груп- пы. К первой группе принадлежат так называемые прямые методы алгоритмы, позволяющие получить решение за конечное число арифметических действий. Сюда относятся известное правило Крамера нахождения решения при помощи определите- лей, метод исключения Гаусса, метод прогонки метод решения систем с трехдиа- гональными матрицами. Существуют и другие методы, из которых отметим методХолецкого (метод квадратных корней), применяемый к системам с симметричными положительно определенными матрицами, метод вращений и метод отражений. Вторую группу составляют приближенные методы, в частности, итерационные. В итерационных методах решение системы получается как предел при стремлении числа итераций n к бесконечности. При конечных n, как правило, получаются лишь приближенные решения. Прямые и итерационные методы имеют свою область применения: если размер- ность системы не слишком велика, то часто предпочтительнее использовать прямые методы. Итерационные методы выгодны для систем большого порядка. Особенно в случае матриц специального вида. В настоящем курсе основное внимание будет уделено прямым методам, а итераци- онных методов коснемся лишь кратко. Более подробно итерационные методы будут изложены во второй части курса "Численные методы", которая читается на четвертом курсе. 6
џ 1 Метод исключения Гаусса и треугольное (LU) разложение матрицы 1.1 Метод исключения Гаусса Система Ax = b в развернутой форме имеет вид n j=1 a ij x j = b i , i = 1, . . . , n. (1.1) Как известно, метод Гаусса или метод последовательного исключения неизвестных состоит в том, что неизвестные x j , j = 1, . . . , n ? 1 последовательно исключаются из соответствующих уравнений системы (1.1) , в результате чего она преобразуется к эквивалентной системе с треугольной матрицей a (0) 11 x 1 + a (0) 12 x 2 + a (0) 13 x 3 + . . . + a (0) 1n x n = b (0) 1 , a (1) 22 x 2 + a (1) 23 x 3 + . . . + a (1) 2n x n = b (1) 2 , a (i?1) ii x i + . . . + a (i?1) in x n = b (i?1) i , a (n?1) nn x n = b (n?1) n , (1.2) коэффициенты a (k) ij которой и компоненты ее правой части b (k) i вычисляются по фор- мулам a (k) ij = a (k?1) ij ? l ik a (k?1) kj , i, j = k + 1, . . . , n; k = 1, . . . , n ? 1; (1.3) 7
8 џ 1. МЕТОД ГАУССА И ТРЕУГОЛЬНОЕ РАЗЛОЖЕНИЕ МАТРИЦЫ b (k) i = b (k?1) i ? l ik b (k?1) k , i, j = k + 1, . . . , n; k = 1, . . . , n ? 1; (1.4) а l ik = a (k?1) ik /a (k?1) kk , i = k + 1, . . . , n; k = 1, . . . , n ? 1; (1.5) причем a (0) ij = a ij , b (0) i = b i Вычисления по формулам (1.3)-(1.5)называются прямым ходом метода Гаусса. После этого неизвестные x k последовательно, начиная с x n , находятся из (1.2) по формулам x i = b (i?1) i ? n j=i+1 a (i?1) ij x j a (i?1) ii , i = n, . . . , 1. (1.6) Вычисления по этим формулам называют обратным ходом метода Гаусса. Замечание 1.1. В формулах (1.2) при пребразованиях системы (1.1) первое уравне- ние осталось без изменения. С равным успехом может быть использован и другой вариант исключения, когда первое уравнение (1.1) делится на a 11 , а вместо (1.2) получается система с единичными коэффициентами при x j в j-ом уравнении. Замечание 1.2. Вычисления по формулам (1.5), (1.6) , а, следовательно, и по фор- мулам (1.3), (1.4) возможны лишь тогда, когда все числа a (i?1) ii = 0, i = 1, . . . , n. (1.7) Необходимые и достаточные условия выполнения (1.7) устанавливаются в доказы- ваемой чуть позже теореме 1.2 1.2 LU разложение матрицы. Покажем, что метод Гаусса эквивалентен разложению матрицы A системы (1.1) в произведение нижней L и верхней U треугольных матриц с последующим решением вспомогательных систем с этими матрицами. В самом деле, из (1.4) находим, что b (k?1) i = b (k?2) i ? l i k?1 b (k?2) k?1 Подставляя это соотношение в (1.4), получим b (k) i = b (k?2) i ? l i k?1 b (k?2) k?1 ? l ik b (k?1) k Точно так же подставляя сюда выражения для b (k?2) i , а затем для b (k?3) i и т.д., будем иметь b (k) i = b (0) i ? l i1 b (0) 1 ? l i2 b (1) 2 ? · · · ? l ik b (k?1) k 1.2. LU РАЗЛОЖЕНИЕ МАТРИЦЫ. 9 Полагая здесь k = i ? 1 и выражая b (0) i = b i через b (j) i , получим b i = i?1 j=1 l ij b (j?1) j + b (i?1) i (1.8) Обозначим столбец правой части системы (1.2) через y = [y 1 . . . y n ] T , полагая y i = b (i?1) i (1.9) В этих обозначениях (1.8) перепишется так b i = i?1 j=1 l ij y j + y i , i = 1, . . . , n. (1.10) Обозначим через L нижнюю треугольную матрицу с коэффициентами l ij , вычисляе- мым по формулам (1.5), и единичной главной диагональю L = ? ? ? ? ? ? 1 0 0 . . . 0 l 21 1 0 . . . 0 l 31 l 32 1 . . . 0 l n1 l n2 l n3 . . . 1 ? ? ? ? ? ? (1.11) Тогда (1.10) можно записать в матричном виде b = Ly. (1.12) Если верхнюю треугольную матрицу системы (1.2) обозначить через U и переписать (1.2) в матричном виде, то будем иметь Ux = y. (1.13) Действуя теперь на левую и правую часть (1.13) невырожденной матрицей L и при- нимая во внимание (1.12), получим LUx = Ly = b ? A = LU. (1.14) Итак, мы показали, что реализация вычислений по формулам (1.3) и (1.5) прямого хода метода Гаусса эквивалентна разложению матрицы A системы (1.1) в произве- дение нижней треугольной матрицы с единичной главной диагональю L и верхней треугольной матрицы U. При этом элементы матрицы L вычисляются по формулам (1.5), а элементы матрицы U суть u kj = a (k?1) kj (1.15)
10 џ 1. МЕТОД ГАУССА И ТРЕУГОЛЬНОЕ РАЗЛОЖЕНИЕ МАТРИЦЫ и вычисляются по формулам (1.3). После разложения матрицы A в произведение двух треугольных для отыскания решения системы (1.1) нужно решить две системы с треугольными матрицами системы (1.12) и (1.13). Решение системы (1.12) заменяет преобразование вектора правой части системы (1.1) по формулам (1.4) прямого хода метода Гаусса. Решение же x системы (1.13) с учетом обозначений (1.9) и (1.15) определяется формулами (1.6) обратного хода метода Гаусса. Замечание 1.3. Соотношения (1.3) содержат формулы для u kj (1.15) и промежуточ- ные значения, которые тоже нужно запоминать и хранить. Мы сейчас преобразуем эти формулы к такому виду, при котором хранение промежуточных значений не требуется. Пусть матрица L имеет вид (1.11), т.е. ее элементы l ik = 0 при k > i, (1.16) а U = ? ? ? ? ? ? u 11 u 12 u 13 . . . u 1n u 22 u 23 . . . u 2n u 33 . . . u 3n u nn ? ? ? ? ? ? , т.е. u kj = 0 при k > j. (1.17) Поскольку LU = A, то по правилу умножения матриц находим, что a ij = n k=1 l ik u kj (1.18) Преобразуем эту формулу двумя способами. В силу (1.11), (1.16) n k=1 l ik u kj = i?1 k=1 l ik u kj + l =1 ii u ij + n k=i+1 0 l ik u kj = = i?1 k=1 l ik u kj + u ij , а в силу (1.17) n k=1 l ik u kj = j?1 k=1 l ik u kj + l ij u jj + n k=j+1 l ik 0 u kj = = j?1 k=1 l ik u kj + l ij u jj 1.2. LU РАЗЛОЖЕНИЕ МАТРИЦЫ. 11 Отсюда и из (1.18) имеем u ij = a ij ? i?1 k=1 l ik u kj i = 1, . . . , n; j = i, . . . , n; l ij = 1 u jj a ij ? j?1 k=1 l ik u kj j = 1, . . . , n; i = j + 1, . . . , n. (1.19) Очевидно, что реализация формул (1.19) возможна только тогда, когда все u jj = a (j?1) jj в силу (1.15) отличны от нуля (ср. с (1.7)). Замечание 1.4. Формулы (1.19) устроены так, что нельзя сначала вычислить все u ij , а затем все l ij или наоборот. Можно предложить следующий порядок вычислений по формулам (1.19): u 1j = a 1j , j = 1, 2, . . . , n; l i1 = a i1 /u 11 , i = 2, 3, . . . , n; u 2j = a 2j ? l 21 u 1j , j = 2, 3, . . . , n; l i2 = (a i2 ? l i1 u 12 )/u 22 , i = 3, 4, . . . , n; и т.д., т.е. чередовать вычисление строк матрицы U и столбцов матрицы L. После построения матриц L и U решение систем (1.12) и (1.13) с треугольными матрицами находятся по формулам y i = b i ? i?1 k=1 l ik y k , i = 1, 2, . . . , n, (1.20) (вычисления ведутся сверху вниз). x k = 1 u kk y k ? n j=k+1 u kj x j , k = n, n ? 1, . . . , 1 (1.21) (вычисления ведутся снизу вверх). Одной из важнейших характеристик любого численного метода является его тру- доемкость. Под трудоемкостью метода, предназначенного для решения системы (1.1), обычно понимают число арифметических действий, необходимых для нахождения искомого решения. Часто в трудоемкость метода включают лишь действия умножения и деления, как наиболее трудоемкие операции с точки зрения работы компьютера. Так будем поступать и мы.
12 џ 1. МЕТОД ГАУССА И ТРЕУГОЛЬНОЕ РАЗЛОЖЕНИЕ МАТРИЦЫ Легко видеть, что для вычислений по формулам (1.19) (для получения треуголь- ного разложения) требуется Q = n i=1 n j=i (i ? 1) + n j=1 n i=j+1 j = n i=1 [(i ? 1)(n ? i + 1) + i(n ? i)] = = 2 n i=1 [(n + 1)i ? i 2 ] ? n(n + 1) = n(n + 1) 2 ? n(n + 1)(2n + 1) 3 ? n(n + 1) = = n(n 2 ? 1) 3 = n 3 3 + O(n) ? n 3 3 (1.22) действий умножения и деления. Для вычислений по формулам (1.20) и (1.21) имеем соответственно ? q = n i=1 (i ? 1) = n(n ? 1) 2 и q = n k=1 (n ? k + 1) = n(n + 1) 2 , т.е. общее число действий для решения систем (1.12) и (1.13) по формулам (1.20), (1.21) есть q = ? q + q = n 2 (1.23) Замечание 1.5. Из формул (1.22) и (1.23) следует, что при больших n основной объем работы, которую нужно выполнить для решения системы (1.1) описанным методом, падает на преобразование коэффициентов матрицы системы, т.е. на постро- ение треугольного разложения, в то время как преобразование вектора правой части (решение системы (1.12)) и на отыскание самого решения трудозатраты сравниетельно невелики. В связи с этим при больших n решение нескольких систем с различными правыми частями и одной и той же матрицей оказывается по трудоемкости практи- чески таким же как и решение одной системы. Выясним теперь условия, при которых вычисления по формулам (1.19)-(1.21) воз- можны, т.е. все u jj отличны от нуля. Теорема 1.1. Пусть A невырожденная матрица, L нижняя треугольная мат- рица с единичной главной диагональю, а U невырожденная верхняя треугольная матрица. Тогда, если A = LU, то это представление единственно. Для доказательства теоремы 1.1 нам потребуется Лемма 1.1. Произведение нижних (верхних)треугольных матриц есть нижняя (верхняя) треугольная матрица. Обратная к невырожденной нижней (верхней) тре- угольной матрице есть нижняя (верхняя) треугольная матрица. Упражнение 1.1. Доказать лемму 1.1.
1.2. LU РАЗЛОЖЕНИЕ МАТРИЦЫ. 13 Доказательство теоремы 1.1. Пусть A = L 1 U 1 = L 2 U 2 . Тогда L 2 = L 1 U 1 U ?1 2 и L ?1 1 L 2 = U 1 U ?1 2 Слева стоит произведение нижних треугольных матриц, а справа верхних. Поэтому произведение есть диагональная матрица D, т.е. L ?1 1 L 2 = D . Отсюда находим, что L 2 = L 1 D . Поскольку главные диагонали L 1 и L 2 единичные, то главная диагональ L 1 D совпадает с главной диагональю D, следовательно, D = I. Отсюда L 1 = L 2 и U 1 = U 2 . Теорема доказана. Теорема 1.2. Пусть A квадратная невырожденная матрица, L нижняя тре- угольная матрица с единичной главной диагональю, а U невырожденная верхняя треугольная матрица. Разложение A = LU существует тогда и только тогда, когда все угловые миноры матрицы A отличны от нуля. Напомним, что угловыми минорами матрицы A называются величины ? 1 = a 11 , ? 2 = det a 11 a 12 a 21 a 22 , . . . , ? n = det[A]. Доказательство. 1 ? . (Необходимость) Пусть разложение A = LU существует. Тогда по теореме 1.1 оно единственно. Представим матрицы A, L и U в блочном виде A = A m A 12 A 21 A 22 , L = L m 0 L 21 L 22 , U = U m U 12 0 U 22 , где A m , L m , U m и A 22 , L 22 , U 22 квадратные матрицы размерностей m Ч m и (n ? m) Ч (n ? m) соответственно, а m < n произвольное число. Разложение A = LU в блочном представлении имеет вид A m A 12 A 21 A 22 = L m 0 L 21 L 22 U m U 12 0 U 22 = L m U m L m U 12 L 21 U m L 21 U 12 + L 22 U 22 (1.24) Отсюда следует, что A m = L m U m (1.25) Поскольку матрица U треугольная и невырожденная, то все ее диагональные элемен- ты отличны от нуля. Поэтому невырождена и треугольная матрица U m . Тем самым ? m = det[A m ] = det[L m ] det[U m ] = u 11 . . . u mm = 0 при m = 1, . . . , n. 2 ? . (Достаточность) Пусть теперь ? 1 ? 2 . . . ? n = 0 . Для доказательства суще- ствования треугольного разложения воспользуемся методом полной математической индукции по порядку системы n. При n = 1 матрица A = a 11 = ? 1 = 0 , матрица L = 1 14 џ 1. МЕТОД ГАУССА И ТРЕУГОЛЬНОЕ РАЗЛОЖЕНИЕ МАТРИЦЫ и поэтому U = u 11 = a 11 = det U = 0 . Существование искомого разложения при n = 1 доказано. Пусть A k матрица порядка k и разложение A k = L k U k существует с det U k = 0 при k = 1, . . . , m ? 1. Докажем, что существует и A m = L m U m , причем det U m = 0 Пусть a ·m = [a 1m . . . a m?1 m ] T столбец, a m· = [a m1 . . . a m m?1 ] строка и разложение A m будем искать в виде A m = A m?1 a ·m a m· a mm = L m?1 0 l m· 1 U m?1 u ·m 0 u mm = = L m?1 U m?1 L m?1 u ·m l m· U m?1 l m· u ·m + u mm Отсюда следует, что неизвестный столбец u ·m , неизвестная строка l m· и u mm опреде- ляются следующими соотношениями: L m?1 u ·m = a ·m l m· U m?1 = a m· ? U T m?1 l T m· = a T m· , u mm = a mm ? l m· u ·m (1.26) Поскольку L m?1 и U m?1 невырожденные, из первого и второго соотношений (1.26) можно найти u ·m и l m· , соответственно, после чего третье соотношение дает u mm Существование разложения (1.24) для m доказано. Осталось доказать, что det U m = 0 Но с учетом (1.25) 0 = ? m = det A m = det U m , что и требовалось доказать. Теорема полностью доказана.
џ 2 Ленточные методы 2.1 Метод прогонки Рассмотрим систему линейных алгебраических уравнений Ax = b, (2.1) матрица которой является трехдиагональной. Запишем эту систему в развернутом виде. Пусть b 1 x 1 + c 1 x 2 = d 1 , a 2 x 1 + b 2 x 2 + c 2 x 3 = d 2 , a i x i?1 + b i x i + c i x i+1 = d i , a n x n?1 + b n x n = d n (2.2) Алгоритм метода прогонки метода решения системы (2.2) состоит в следующем (см. курс "Введение в численные методы", но, может быть, с другими обозначениями!) а) Нахождение прогоночных коэффициентов (прямая прогонка) по формулам ? i = ?c i /? i , i = 1, 2, . . . , n ? 1, ? i = b i + a i ? i?1 , i = 2, . . . , n, ? 1 = b 1 , ? i = (d i ? a i ? i?1 )/? i , i = 2, . . . , n, ? 1 = d 1 /b 1 (2.3) б) Нахождение самого решения (обратная прогонка) x i = ? i x i+1 + ? i , i = n ? 1, . . . , 1, x n = ? n (2.4) 15
16 џ 2. ЛЕНТОЧНЫЕ МЕТОДЫ Из (2.3),(2.4) следует, что общее число умножений и делений при вычислении коэф- фициентов ? i и ? i Q = 2(n ? 1, ) (2.5) а при вычислении коэффициентов ? i и решения x i q = 3(n ? 1). (2.6) Сравнение (2.5),(2.6) с (1.22), (1.23) свидетельствуют о том, что прогонка существенно менее трудоемка по сравнению с общим методом Гаусса. Связано это с тем, что мы явным образом воспользовались тем, что значительная часть элементов матрицы A равна нулю. 2.2 Ленточные матрицы Определение 2.1. Матрица A называется ленточной с полушириной ленты p, если ее элементы a ij = 0 при |i ? j| > p, но существует по крайней мере один элемент a ij = 0 при |i ? j| = p Пример 2.1. Для диагональной матрицы a ij = 0 при |i ? j| > 0 и, следовательно, ее полуширина равна нулю. Ее лента состоит из одной диагонали и ширина равна 1. Условно диагональную матрицу можно изобразить как на рис. 1. ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Рис. 1. Пример 2.2. Полная матрица имеет 2n?1 диагоналей. Это и есть ширина ее ленты, а полуширина будет p = n ? 1. ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Рис. 2. 2.3. ЛЕНТОЧНЫЙ ВАРИАНТ ТРЕУГОЛЬНОГО РАЗЛОЖЕНИЯ 17 Пример 2.3. На рис. 3 изображена матрица с полушириной p = 1. ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Рис. 3. Матрицы такой структуры называются трехдиагональными. Ширина ленты 3. Пример 2.4. Матрица, изображенная на рис.4, имеет полуширину p = 1 и ширину ленты 2. Матрицы такой структуры называются трапециевидными. Изображенная матрица также называется правой ленточной. ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Рис. 4. Пример 2.5. У матрицы на рис. 5 ? ? ? ? ? ? ? ? ? 0 ? 0 ? 0 ? ? 0 ? 0 ? ? 0 ? 0 ? ? 0 ? 0 ? 0 ? ? ? ? ? ? ? ? ? Рис. 5. полуширина p = 2, ширина 5 и всего три ненулевых диагонали. Определение 2.2. Рисунок, на котором (звездочками) отмечены позиции, где толь- ко и могут располагаться ненулевые элементы матрицы A, называется ее портретом. 2.3 Ленточный вариант треугольного разложения Модифицируем алгоритм исключения Гаусса на случай ленточных матриц, т.е. зара- нее отбросим те вычисления, которые заведомо приводят к нулевым элементам. Это
18 џ 2. ЛЕНТОЧНЫЕ МЕТОДЫ позволит нам сэкономить в трудозатратах на решение системы. Обратимся сразу к варианту, основанному на треугольном разложении матрицы A. Нам потребуется Лемма 2.1. Если полуширина матрицы A равна p, то в треугольном разложении A = LU полуширина L (U) не больше p. Доказательство. Пусть a ij = 0 при |i ? j| > p. (2.7) Докажем, что l ij = 0 при i ? j > p. (2.8) Для доказательства применим метод полной математической индукции по номерам столбцов матрицы L. При j = 1 из (1.19) находим, что l i1 = a i1 /a 11 , i = 2, . . . , n. Отсюда с учетом (2.7) приходим к (2.8) с j = 1, т.е. l i1 = 0 при i ? 1 > p. Пусть теперь утверждение (2.8) верно для столбцов матрицы L с номерами k = 1, 2, . . . , j ? 1, т.е. l ik = 0 при i ? k > p, k = 1, 2, . . . , j ? 1. (2.9) Докажем справедливость (2.8) для j-ого столбца. Пусть i ? j > p. (2.10) Тогда в силу (1.19) и (2.7) l ij = 1 u jj 0 a ij ? j?1 k=1 l ik u kj = ? 1 u jj j?1 k=1 l ik u kj (2.11) Оценим разность индексов i ? k у первых сомножителей под знаком суммы в (2.11). С учетом (2.10) и (2.11) будем иметь i ? k > p + j ? k p + 1. Но тогда в силу (2.9) первые сомножители в сумме (2.11) обращаются в нуль и соотношение (2.8) установлено. Лемма доказана. Преобразуем формулы (1.19)-(1.21) на случай ленточной матрицы A. Сначала вы- ясним, для каких значений индексов i и j нужно проводить вычисления по формулам (1.19). Так как в силу леммы 2.1 u ij = 0 при j ? i > p, (2.12) 2.3. ЛЕНТОЧНЫЙ ВАРИАНТ ТРЕУГОЛЬНОГО РАЗЛОЖЕНИЯ 19 то ненулевые элементы u ij могут быть лишь при j?i p , т.е. при j p+i . Аналогично l ij = 0 при i ? j > p (2.13) и, следовательно, ненулевые элементы l ij могут быть только при i ? j p , т.е. при i p + j . Отсюда u ij = 0 i = 1, . . . , n, j = i, . . . , min[n, p + i], l ij = 0 j = 1, . . . , n, i = j + 1, . . . , min[n, p + j]. Теперь преобразуем суммы в (1.19). В силу (2.12) ненулевые слагаемые в суммах (1.19) могут быть только при j ? k p , т.е. при k j ? p. а в силу (2.13) только при i ? k p , т.е. при k i ? p. Объединяя эти неравенства и принимая во внимание, что k натуральное, будем иметь k max[1, i ? p, j ? p]. Поскольку в формулах для u ij индексы i, j подчинены ограничению j i , то в этих формулах k max[1, j ? p]. В формулах же для l ij наоборот i > j и поэтому в них k max[1, i ? p]. С учетом сказанного, для ленточной матрицы A с полушириной p формулы (1.19) принимают вид u ij = a ij ? i?1 k=max[1,j?p] l ik u
Поделитесь с Вашими друзьями: |