x ( k+1) достигают предельных значений (конечно, в пределах определенного количества десятичных разрядов). Следовательно, использование рекуррентного процесса (4.48) позволяет за конечное число итераций получить значения корней исходной системы (4.1) с заданной погрешностью. Метод Зейделя иногда дает более быструю сходимость по сравнению с методом простой итерации. Условия сходимости этих методов, вообще говоря, не совпадают. Существуют такие матрицы коэффициентов при неизвестных, для которых метод Зейделя сходится, а метод простой итерации – нет. Встречаются и противоположные ситуации. Для формулировки условий сходимости метода Зейделя целесообразно ввести понятие нормальной системы линейных уравнений. Система линейных уравнений вида (4.1) называется нормальной, если выполняются два следующих условия. 1. Матрица А коэффициентов при неизвестных симметрична, т.е. a ij =a ji 2. Квадратичная форма U, соответствующая матрице А, U = ∑ = n i 1 ∑ = n j j i ij x x a 1 (4.49) 81 положительно определена (т.е. принимает при любых значениях аргументов x i неотрицательные значения, причем обращается в нуль, только когда все аргументы x i равны нулю). Если матрица коэффициентов А = [ a ij ] неособенная (т.е. det A ≠ 0), то система линейных уравнений общего вида (4.1) может быть приведена к нормальному виду. Для этого обе части матричного уравнения (4.3) надо умножить слева на транспонированную матрицу А ′′′′ = [ a ji ]: А ′′′′ А x = А′′′′ b. (4.50) Произведение квадратных матриц А ′′′′А представляет собой тоже матрицу А N размером n × n. В курсе линейной алгебры доказывается, что матрица А N симметрическая и соответствует положительно определенной квадратичной форме. Произведение А ′′′′b = b N является n- мерным вектором. Таким образом, нормализация исходной системы линейных уравнений (4.1) дает эквивалентную систему вида А N x = b N (4.51) После получения нормальной системы (4.51) можно ее преобразовать к виду (4.35). В линейной алгебре доказана следующая теорема: если система линейных уравнений вида (4.1) − нормальная, то процесс Зейделя для эквивалентной системы (4.35) сходится к точному решению системы при любом выборе начального приближения. Пример 5. Рассмотрим следующую систему трех линейных уравнений: 9 x 1 − 2 x 2 + x 3 = 8, 2 x 1 − 7 x 2 + x 3 = −4, (4.52) x 1 + 3 x 2 + 8 x 3 = 12. Легко видеть, что для уравнений системы (4.52) выполнены условия (4.44), поэтому возможно решить данную систему итерационным методом. Используем метод Зейделя. Приведем систему (4.52) к виду (4.35), необходимому для проведения итераций:
82 x1 = 0,88889 + 0,22222 x2 − 0,11111 x3 , x2 = 0,57142 + 0,28571 x1 + 0,14286 x3 , (4.53) x3 = 1,5 − 0,125 x1 − 0,375 x2 Так как выполнение условий (4.44) обеспечивает сходимость при любом начальном приближении, в качестве такового выберем нулевое: ) 0 ( 1 x=0; ) 0 ( 2 x=0; ) 0 ( 3 x=0. Подставим эти значения в правую часть первого уравнения системы (4.53) и легко получим первое приближение первого корня ) 1 ( 1 x=0,88889. При вычислении первого приближения второго корня ) 1 ( 2 x используем уже полученное значение ) 1 ( 1 x: ) 1 ( 2 x= 0,57142 + 0,28571 ) 1 ( 1 x + 0,14286 ) 0 ( 3 x = 0,82539. Аналогично вычислим первое приближение третьего корня ) 1 ( 3 x= 1,5 − 0,125 ) 1 ( 1 x− 0,375 ) 1 ( 2 x = 1,07937. Продолжим итерации методом Зейделя, а результаты запишем в таблицу 4.2. В том, что точным решением системы (4.52) является x1 = x2 = x3 = 1, легко убедиться подстановкой полученных значений в исходную систему. Таблица 4.2 Номер итерации x1 x2 x3 0 0 0 0 1 0,88889 0,82539 1,07937 2 0,95238 0,99773 1,0068 3 0,99874 1,00061 0,99993 4 1,00014 1,00003 0,99997 5 1,00001 0,999999 0,999999 6 1,000000 1,000000 1,000000 83 Из таблицы 4.2 видно, что значения всех трех корней быстро сходятся к единице. Точность до шестого знака после запятой достигается на шестой итерации.
84 Глава 5. ВЫЧИСЛЕНИЕ ДЕТЕРМИНАНТОВ Каждая квадратная матрица A характеризуется определенным числом det A, которое называется детерминантом или определителем этой матрицы. Порядок детерминанта совпадает с порядком соответствующей матрицы. Как следует из предыдущей главы, при решении системы линейных уравнений полезно (или даже необходимо) знать численное значение детерминанта матрицы коэффициентов при неизвестных. Значения детерминантов матриц требуются и в других прикладных задачах физики. Поэтому целесообразно в этой главе рассмотреть вопрос о вычислении детерминантов. В линейной алгебре величина детерминанта определяется индуктивно. Квадратная матрица порядка n = 1 состоит из единственного элемента a11 . Детерминант такой матрицы принимается равным этому элементу. Детерминант квадратной матрицы порядка n = 2 определяется следующей формулой: det A = a11 ⋅ a22 − a12 ⋅ a21 (5.1) Детерминанты порядка n > 2 определяются через детерминанты более низких порядков с помощью разложения по минорам. Как известно из линейной алгебры, каждому элементу aij квадратной матрицы A порядка n ставится в соответствие минор Mij. Минор Mijпредставляет собой детерминант матрицы Aij, которая получается из исходной A вычеркиванием i-й строки и j-го столбца. Ясно, что каждая матрица Aij и соответствующий минор Mij имеют порядок n − 1. В линейной алгебре доказана теорема, что детерминант n-го порядка может быть вычислен по следующей формуле: det A = ∑ = + − njji1 ) 1 ( aij⋅ Mij, (5.2) где номер строки может принимать любое значение i = 1, …, n. Формула (5.2) называется разложением детерминанта по минорам i-й строки. Справедлива и другая формула вычисления детерминанта n-го порядка разложением по минорам j-го столбца: 85 det A = ∑ = + − niji1 ) 1 ( aij⋅ Mij, (5.3) где номер j может быть выбран любым из множества: 1, …, n. Таким образом, для вычисления одного детерминанта n-го порядка необходимо вычислить n детерминантов ( n – 1)-го порядка. Для вычисления каждого минора можно вновь применить формулу (5.2). В результате, чтобы вычислить один детерминант n-го порядка этим методом приходится рассчитать n!/2 детерминантов 2-го порядка. Кроме того, необходимо выполнить умножения миноров на коэффициенты aij исходной матрицы и провести суммирования по формулам (5.2) или (5.3). Ясно, что для высокого порядка вычисление детерминантов разложением по минорам является длительной операцией даже при использовании персональных компьютеров. В линейной алгебре также используется метод выражения детерминантов непосредственно через коэффициенты матрицы aij без промежуточного вычисления миноров. Этот метод можно представить состоящим из следующих этапов. Вначале составляется произведение коэффициентов aij вида a1 α a2 β ...anω , где индексы α, β, ..., ω взяты из множества целых чисел (1, 2, ..., n) и все различны. В качестве первого набора можно взять α = 1, β=2, ..., ω = n. Затем перебираются всевозможные перестановки индексов α, β, ..., ω в произведениях a1 α a2 β ...anω , начиная с исходной. Если перестановка нечетная, то произведение a1 α a2 β ...anω умножается на −1. После перебора всех перестановок индексов α, β, ..., ω все полученные произведения суммируются. Заметим, что количество перестановок равно n!. Справедливость такого представления для детерминанта 2-го порядка видна непосредственно из формулы (5.1). Детерминант 3-го порядка, построенный последним методом, дает следующее выражение: det A = a11 a22 a33 − a12 a21 a33 + a13 a21 a32 − a11 a23 a32 + + a12 a23 a31 − a13 a22 a31 , (5.4) что совпадает с результатом, полученным по методу разложения по минорам. Пример 1. Необходимо вычислить детерминант следующей матрицы размера 4 ×4: 86 det A = − − − − − 1 2 0 1 6 3 3 2 1 1 1 1 2 0 4 5 (5.5) Разложим данный детерминант 4-го порядка на миноры по первой строке. Применив формулу (5.2) и пользуясь тем, что элемент a 13 ⋅ = 0, получим разложение: det A = 5 1 2 0 6 3 3 1 1 1 − − − −(−4) 1 2 1 6 3 2 1 1 1 − − − − −2 2 0 1 3 3 2 1 1 1 − . (5.6) Далее перед нами два разных пути. Можно каждый из трех полученных детерминантов 3-го порядка вычислить по формуле (5.4). Иначе можно детерминанты 3-го порядка предварительно разложить и свести к вычислению нескольких детерминантов 2-го порядка. Если проводить разложение по минорам нижних строк, то (из-за того, что a 42 = 0) для вычисления первого и третьего детерминантов в (5.6) потребуется вычислить только по два детерминанта 2-го порядка. Предлагаем читателям провести необходимые выкладки и убедиться, что det A = −6. Описанный алгоритм может быть запрограммирован и реализован на персональном компьютере. Нетрудно подсчитать, что для вычисления детерминанта n-го порядка необходимо выполнить n ⋅n! умножений и n! сложений (или вычитаний). Следовательно, метод разложения по минорам практически не пригоден для вычисления детерминантов высоких порядков. Более эффективный метод (требующий гораздо меньшего количества операций) базируется на следующей теореме линейной алгебры. Если к элементам любой строки (или любого столбца) детерминанта прибавить соответствующие элементы другой строки ( или столбца), умноженные на произвольный множитель, то величина детерминанта не изменится. В частности, подобными операциями исходная квадратная матрица приводится к треугольному виду. Способ
87 приведения квадратной матрицы к треугольному виду аналогичен процедуре прямого хода в методе Гаусса при решении системы линейных уравнений. Процедура начинается с выбора ведущего элемента. Это может быть любой коэффициент ai1 , не равный нулю. Переставляя строки (или столбцы) детерминанта, всегда можно ведущим элементом сделать коэффициент a11 . При этом следует помнить, что при перестановке двух строк (или двух столбцов) знак детерминанта изменяется на противоположный. Затем коэффициент a11 выносится из первого столбца в качестве общего множителя детерминанта. Вычитаем из элементов aij каждого j- го столбца ( j ≥ 2) соответствующие элементы ai1 первого столбца, умноженные на a1 j. Преобразование элементов можно записать в виде: ) 1 ( ija = ija− 1 ia11 1 aaj, i = 1, ..., n, j = 2, ..., n. (5.7) В преобразованном детерминанте элемент ) 1 ( 11 a = 1. Остальные элементы первой строки равны нулю. Разложение этого детерминанта по элементам первой строки, согласно (5.2), показывает, что он равен детерминанту порядка ( n − 1), составленному из элементов ) 1 ( ija ( i, j = 2, ..., n ). Таким образом, исходный детерминант n-го порядка выразится в виде произведения: det A = a11 ∆ n− 1 , (5.8) где ∆ n− 1 = ) 1 ( ) 1 ( 3 ) 1 ( 2 ) 1 ( 3 ) 1 ( 33 ) 1 ( 32 ) 1 ( 2 ) 1 ( 23 ) 1 ( 22 ... ... ... nnnnnnaaaaaaaaa (5.9) новый детерминант ( n − 1)-го порядка. Следующий шаг заключается в применении к новому детерминанту ∆ n− 1 той же процедуры. Выбирается новый ведущий элемент, не равный нулю. Переставляя при необходимости строки или столбцы нового 88 детерминанта, делаем ведущим элемент ) 1 ( 22 a. Выносим его из первого столбца нового детерминанта в качестве общего множителя. Процедура продолжается, причем на каждом шаге данного алгоритма уменьшается порядок детерминанта. Преобразования элементов приводят к тому, что первая строка каждого нового детерминанта начинается с единицы. Остальные элементы этой строки – нули. В результате, если на каждом шаге удается найти ненулевой ведущий элемент, искомый детерминант n-го порядка представится в виде произведения этих ведущих элементов: det A = ) 1 ( ) 1 ( 22 11 − ⋅ ⋅ ⋅ nnnaaa (5.10) Этот способ вычисления детерминанта называется методом Гаусса. Если на каком-либо шаге не удается отыскать ненулевой ведущий элемент, то это значит, что детерминант исходной матрицы равен нулю. В терминологии матричной алгебры можно сказать, что проведенное преобразование исходной квадратной матрицы к треугольной позволяет выразить ее детерминант в виде произведения диагональных элементов полученной треугольной матрицы. Можно подсчитать, что количество арифметических операций для вычисления детерминанта n-го порядка методом Гаусса не превышает n3 . Следовательно, метод Гаусса предпочтительнее рассмотренных выше уже при n > 5. Пример 2. Вычислим детерминант (5.5), приведенный в примере 1 данной главы, но на этот раз используем метод Гаусса, изложенный выше. Первый ведущий элемент a11 = 5. Преобразуем 2-ю строку данной матрицы по формуле (5.7): ) 1 ( 2 ja= ja2 − 21 a11 1 aaj, j = 2, 3, 4 (номер столбца). ) 1 ( 22 a=1–(–1) ⋅(–4)/5=1/5, ) 1 ( 23 a=1–(–1) ⋅0/5=1, ) 1 ( 24 a=–1 – (–1) 2/5=–3/5. Аналогично преобразуются 3-я и 4-я строки данной матрицы: 89 ) 1 ( 3 j a = j a 3 − 31 a 11 1 a a j , ) 1 ( 4 j a = j a 4 − 41 a 11 1 a a j , j = 2, 3, 4. Получаем новый детерминант ∆ 3 = ) 1 ( 44 ) 1 ( 43 ) 1 ( 42 ) 1 ( 34 ) 1 ( 33 ) 1 ( 32 ) 1 ( 24 ) 1 ( 23 ) 1 ( 22 a a a a a a a a a = 5 / 7 2 5 / 4 5 / 34 3 5 / 23 5 / 3 1 5 / 1 − − − Новый ведущий элемент ) 1 ( 22 a = 1/5. Предпоследняя и последняя строки детерминанта ∆ 3 преобразуются по тому же алгоритму (см. первую формулу в (4.15)): ) 1 ( 22 ) 1 ( 2 ) 1 ( 32 ) 1 ( 3 ) 2 ( 3 a a a a a j j j − = , ) 1 ( 22 ) 1 ( 2 ) 1 ( 42 ) 1 ( 4 ) 2 ( 4 a a a a a j j j − = , j = 3, 4. Следующий детерминант получает вид: ∆ 2 = ) 2 ( 44 ) 2 ( 43 ) 2 ( 34 ) 2 ( 33 a a a a = 1 2 7 20 − − Теперь вычислим последний сомножитель произведения (5.10): ) 1 ( − n nn a = ) 3 ( 33 ) 2 ( 34 ) 2 ( 43 ) 2 ( 44 ) 3 ( 44 a a a a a − = = 1–(–2) ⋅7/(– 20) = 3/10. Искомый детерминант равен det A = a 11 ) 1 ( 22 a ) 2 ( 33 a ) 3 ( 44 a =5 ⋅ 5 1 ⋅(–20)⋅ 10 3 = –6. Конечно, детерминант ∆ 2 можно было бы сосчитать по формуле (5.1) и представить результат в виде произведения det A = a 11 ⋅ ) 1 ( 22 a ⋅∆ 2 Искомое значение det A = –6 достигается еще быстрее.
90 Глава 6. ОБРАЩЕНИЕ МАТРИЦ Как известно из линейной алгебры, обратной матрицей A−−−− 1 по отношению к данной A называется матрица, которая будучи умноженной на A как справа, так и слева, дает единичную матрицу E: A⋅⋅⋅⋅ A−−−− 1 = A−−−− 1⋅ A = E. (6.1) Единичную матрицу E удобно кратко записывать с помощью символов Кронекера: E = [ δ ij ], (6.2) где δ ij = ≠ = ) ( 0 ) ( 1 jiji (6.3) Матрицы, обратные данным, используются во многих прикладных задачах, в т.ч. и физических. Следовательно, существует потребность в методах, позволяющих вычислять элементы обратной матрицы, выражая их через элементы исходной. В курсе линейной алгебры доказано, что обратная матрица A−−−− 1существует, если исходная матрица A является неособенной, т.е. ее детерминант не равен нулю (det A ≠ 0). Методы вычисления детерминантов приведены в предыдущей главе. Одна из теорем линейной алгебры утверждает, что обратная матрица A−−−− 1 может быть представлена как транспонированная матрица алгебраических дополнений исходной A, причем каждый элемент этой транспонированной матрицы должен быть поделен на детерминант det A. Алгебраическое дополнение Aij любого элемента матрицы aijвыражается через его минор Mij: Aij = ( −1) i+ j⋅ Mij (6.4) Таким образом, обратная матрица A−−−− 1 выражается в следующем виде: 91 A−−−− 1 = ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ nnnnnnAAAAAAAAA... ... ... 1 1 2 22 12 1 21 11 , (6.5) где использовано обозначение ∆ = det A. В принципе обратную матрицу можно вычислить по формуле (6.5). Операция транспонирования не представляет трудности. Но для вычисления каждого алгебраического дополнения необходимо вычислить минор − детерминант порядка ( n − 1). Следовательно, для расчета обратной матрицы размером ( n× n) необходимо вычислить один детерминант n-го порядка и n2 детерминантов ( n−1)-го порядка. При больших значениях n вычисление элементов обратной матрицы по формуле (6.5) требует значительного времени, даже при использовании современных персональных компьютеров, поэтому описанный метод обращения матриц практически не используется. За годы развития прикладной математики было разработано много различных способов обращения матриц: разбиение на клетки, окаймление, представление произведением треугольных матриц (LU- факторизация), схема Халецкого, и т.д. Эффективность некоторых методов повышается при специфическом виде обращаемой матрицы. Для плохо обусловленных матриц, т.е. тех, у которых число обусловленности (4.32) много больше единицы, разработаны специальные методы обращения, обеспечивающие достаточную точность вычислений. В эпоху распространения персональных компьютеров одним из широко используемых методов обращения хорошо обусловленных матриц является метод Гаусса. Пусть нам дана хорошо обусловленная квадратная матрица A порядка n вида (4.2). Обозначим через xij элементы обратной матрицы A−−−− Поделитесь с Вашими друзьями: |