Другие статьи


Операции с матрицами на C++.   Класс DMatrix



 

 

Applied Numerical Methods

( Прикладные численные методы )

 

B.Carnahan, H.A.Luther, J.O.Wilkes

 

1969

 

(фрагмент книги)

 

 

 

 

Предисловие переводчика

 

 

В книге B.Carnahan, H.A.Luther, J.O.WilkesApplied Numerical Methods” есть интересный фрагмент, содержащий оригинальную интерпретацию итерационного метода Качмажа, а именно – с использованием ортогонализации Грама-Шмидта. Несмотря на то, что впервые описание метода было опубликовано польским математиком С.Качмажем в 1937 году и несмотря на его высокую эффективность при прикладных расчетах, метод продолжает оставаться не очень широко известным и включается в курс высшей математики далеко не в каждом техническом ВУЗе.

Предлагаю Вам познакомиться с данным фрагментом книги в русском переводе. Первые два переведенных параграфа посвящены базовым сведениям по преобразованиям матриц; далее следует описание метода (с доказательством) и пример.

 

 

 

5.1  Введение

 

 

В главе 5 мы расскажем о методах решения совместных систем из n уравнений относительно n неизвестных  x1, … ,  xn:

 

f1(x1, … ,  xn) = 0,

f2(x1, … ,  xn) = 0,

…………………...

……………………

fn(x1, … ,  xn) = 0.

 

 

 

 

(5.1)

 

 

Общий случай, когда совокупность функций f1, … ,  fn не предполагает возможности упрощения, будет рассмотрен далее, в п.п. 5.8 и 5.9. Если же эти функции линейны по xi, система (5.1) может быть переписана в виде:

 

 

b11x1 + b12x2 + … + b1nxn = u1 ,

b11x1 + b12x2 + … + b1nxn = u1 ,

…………………………………

…………………………………

b11x1 + b12x2 + … + b1nxn = u1 ,

 

 

 

(5.2)

 

 

или более лаконично:

 

Bx = u ,                               (5.3)

 

где B – матрица коэффициентов,  u = [u1, u2, … , un]T – вектор правой части системы,  x = [x1, x2, … , xn]T – вектор решений.

Прямые методы решения системы (5.2), приводящие за конечное число шагов к точному решению (без учета незначительной ошибки, связанной с округлением), мы рассмотрим в п.п. 5.3, 5.4 и 5.5. Но прямые методы имеют практическое применение лишь в случае, если система содержит не очень большое количество уравнений (как правило, порядка 40 или менее). Итеративные методы для решения (5.2) рассмотрены в первом приближении в п.п. 5.6 и 5.7. Итеративные методы более эффективны, когда мы имеем дело с большим количеством совместных уравнений (как правило, порядка 100 и более); такие системы уравнений уже, как правило, обладают определенными специальными характеристиками.

 

 

 

5.2  Элементарные преобразования матриц

 

 

Прежде чем приступить к рассмотрению систем уравнений, вспомним о трех типах элементарных матриц:

1.     Элементарная матрица первого типа – это диагональная матрица Q размера n*n, полученная из единичной матрицы I заменой i-го диагонального элемента ненулевой константой q. Например, при n = 4 и i = 3:

Заметим, что det Q = q ,  а обратную матрицу  Q-1 = diag(1, 1, 1/q, 1) можно получить также из единичной матрицы с помощью замены i-ой единицы на 1/q.

 

2.     Элементарная матрица второго типа – это матрица R размера n*n, полученная из единичной матрицы I перестановкой каких-либо двух строк. Например, при n = 4, i = 1 и j = 3:

 

В этом случае det R = -1 и матрица является обратной к самой себе, т.е. RR = I.

 

3.     Элементарная матрица третьего типа – это матрица S размера n*n, полученная заменой одного из недиагональных элементов (i, j) (i j) матрицы I ненулевой константой. (Можно описать получение этой матрицы по-другому: берем матрицу I и прибавляем строку i, умноженную на s, к строке j.) Например, при n = 4, i = 3 и j = 1:

Заметим, что det Q = 1.

 

Умножение слева произвольной матрицы A размера n*p на одну из перечисленных элементарных матриц называется элементарным преобразованием матрицы A или, иначе говоря, элементарной операцией над строками матрицы A. В качестве примера рассмотрим произведения QA, RA и SA для n = 3, i = 2, j = 3 и p = 4.

 

1.

2.

3.

 

 

Мы видим, что умножение слева на элементарные матрицы приводит к следующим преобразованиям матрицы A:

1. QA:   Умножение всех элементов строки на число.

2. RA:   Обмен местами двух строк.

3. SA:   Прибавление умноженных на число элементов одной строки к соответствующим элементам другой строки.

 

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

Умножение справа произвольной матрицы A размера p*n на одну из элементарных матриц называется элементарной операцией над столбцами. Три типа операций приводят к следующим результатам:

1. AQ:   Умножение всех элементов столбца на число.

2. AR:   Обмен местами двух столбцов.

3. AS:   Прибавление умноженных на число элементов одного столбца к соответствующим элементам другого столбца.

 

Если матрица A – произвольная, а матрица T – результат выполнения элементарных операций над строками или столбцами матрицы A, то матрицы T и A называются эквивалентными.

Из сказанного ранее следуют формулы для квадратной матрицы A:

det (QA) = det Q * det A = q det A

det (RA) = det R * det A = - det A

det (SA) = det S * det A = det A

 

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

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

 

 

 

Примечание переводчика.

 

В следующем параграфе используется представление об ортогональных матрицах и метод ортогонализации Грама-Шмидта.

Ортогональной называется квадратная матрица A прядка n, обладающая следующим свойством:  AAT = E, где AT – транспонированная матрица A (а следовательно, и ATA = E).

Строки и столбцы ортогональной матрицы образуют системы ортонормированных векторов, то есть   (ai, aj) = 0 при i j,   (ai, ai) = 1,  где ai – строки (столбцы) матрицы. Определитель ортогональной матрицы равен 1 или -1.

В приведенном ниже описании алгоритма Качмажа используется обозначение A*; обычно символом «*» обозначают матрицу сопряженного оператора, но в случае ортонормированного базиса такая матрица совпадает с транспонированной.

Процесс ортогонализации – это алгоритм построения для данной линейно независимой системы векторов евклидова или эрмитова пространства V ортогональной системы ненулевых векторов, порождающих то же самое подпространство в V. Наиболее известным является процесс ортогонализации Грама – Шмидта:

Пусть имеются линейно независимые векторы  a1, …, an.

Определим оператор проекции следующим образом:

projba = b,



где (a, b) — скалярное произведение векторов a и b. Этот оператор проецирует вектор a ортогонально на вектор b.

Классический процесс Грама — Шмидта выглядит так:

 \begin{array}{lclr}
\mathbf{b}_1 & = & \mathbf{a}_1 & (1) \\
\mathbf{b}_2 & = & \mathbf{a}_2-\mathbf{proj}_{\mathbf{b}_1}\,\mathbf{a}_2 & (2) \\
\mathbf{b}_3 & = & \mathbf{a}_3-\mathbf{proj}_{\mathbf{b}_1}\,\mathbf{a}_3-\mathbf{proj}_{\mathbf{b}_2}\,\mathbf{a}_3 & (3) \\
\mathbf{b}_4 & = & \mathbf{a}_4-\mathbf{proj}_{\mathbf{b}_1}\,\mathbf{a}_4-\mathbf{proj}_{\mathbf{b}_2}\,\mathbf{a}_4-\mathbf{proj}_{\mathbf{b}_3}\,\mathbf{a}_4 & (4) \\
& \vdots & & \\
\mathbf{b}_N & = & \mathbf{a}_N-\displaystyle\sum_{j=1}^{N-1}\mathbf{proj}_{\mathbf{b}_j}\,\mathbf{a}_N & (N)
\end{array}

 

На основе каждого вектора bj (j = 1, … , N) может быть получен нормированный вектор    ej = bj /  ||bj|| .

В результате выполнения процесса Грама — Шмидта получаем систему ортогональных векторов b1, …, bn, либо систему ортонормированных векторов e1, …, en .

Вычисление b1, …, bn носит название ортогонализации Грама — Шмидта, а e1, …, en  ортонормализации Грама — Шмидта.

 

 

 

5.5  Алгоритм Качмажа в конечной форме

 

 

Рассмотрим систему линейных уравнений:

 

 

b11x1 + b12x2 + … + b1nxn = u1 ,

b11x1 + b12x2 + … + b1nxn = u1 ,

…………………………………

…………………………………

b11x1 + b12x2 + … + b1nxn = u1 ,

 

 

 

(5.2)

 

 

или, кратко:                           Bx = u.

 

Решим систему (5.2) в предположении, что матрица B – невырожденная.

Способ решения заключается в предварительном преобразовании системы (5.2) в эквивалентную систему

Ax = v,          A* = A-1,                (5.9)

то есть системы (5.9) и (5.2) имеют один и тот же вектор решений x, A ортогональная матрица. Затем, произвольно определив значение начального вектора r0, мы вводим рекуррентную формулу:

rj = rj-1 – [(αj, rj-1) - νj] αj ,             1 ≤ j n                (5.10)

где A = (aij) и αj = [aj1, aj2, … , ajn]T.

Далее докажем, что  Arn = v  и, следовательно,

Brn = u.                                     (5.11)

 

Покажем, что выражение (5.10) определяет такой вектор rn, что Arn = v. Уравнения (5.9) можно записать в виде:

 

(αi, x) = νi ,             1 ≤ i n.

 

Умножая (5.10) слева на αj и учитывая ортогональность м-цы A, получим:

 

(αj, rj) = νj ,

 

то есть rj удовлетворяет j-му уравнению системы (5.9). С другой стороны, умножение (5.10) на αi при  i j дает нам

 

(αi, rj) = (αi, rj-1),       (i ≠ j)

 

так как (αi, αj) = 0 при i j. Отсюда по индукции делаем вывод, что Arn = v ,  то есть единственное решение (5.9) найдено.

 

Вернемся к решению системы (5.2). Пусть B* = [β1, β2, … , βn], где βi = [bi1, bi2, … , bin]T. Система, эквивалентная (5.2) и имеющая свойства (5.9), строится на основе линейно независимых векторов βi традиционным способом, а именно – с помощью процедуры ортогонализации Грама-Шмидта. Пусть

 

 

Отсюда с помощью индукции нетрудно вывести, что, при i j, (αi, αi) = 1,  и (αi, αj) = 0.  Следовательно, если A* = [α1, α2, … , αn], то A* = A-1, что проверяется непосредственным умножением.

Далее, положим:

 

 

Определенные таким образом αi и νi  дают нам систему (5.9), имеющую тот же вектор решений, что и (5.2). Ибо, если (αi, x) - νi = 0, то из (5.12) имеем:

 

или

 

С помощью (5.13), получаем:

 

(βj, x) = uj.

 

На практике расчет может проводиться последовательным движением по массиву [B|u] с одновременным вычислением соответствующих строк массива [A|v]. Если необходимо уточнить вектор u после того, как матрица A уже была построена, нужно записать n чисел  и (n2 - n) / 2 чисел (βj, αi) ,  1 ≤ i < j,  2 ≤ j n. Построение матрицы [A|v] на основе [B|u] лучше проводить сначала с помощью вычисления (5.12)  (а не (5.13)). Первая строка [A|v] формируется по первой строке [B|u], вторая – по второй строке [B|u] и уже построенной строке [A|v] и т.д. Каждая операция может быть рассмотрена как умножение слева на элементарную матрицу первого или третьего типа. Таким образом, существует невырожденная матрица φ, такая, что:

 

φ[B|u] = [A|v].

 

Вычисление детерминанта этой матрицы приводит к формуле:

 

 

Знание этого факта может быть полезным в том случае, когда матрица B – плохо обусловленная, то есть содержит строки или столбцы со столь слабо выраженным свойством линейной независимости, что ошибки округления или усечения могут приводить к заметному отклонению рассчитанных решений от их истинных значений. Вспомним, что |det(A)| = 1, из чего следует, что (5.14) приводит к разрешению этой проблемы. Заметим также, что последовательности {αj}, {νj} и {rj} можно рассчитывать последовательно по шагам, поэтому метод можно считать n-шаговым.

Также после ортогонализации вектор решений x можно найти как A* v.

 

Пример. В качестве простой иллюстрации метода Качмажа рассмотрим следующий пример (он будет также рассматриваться в п.п. 5.6 и 5.7):

 

4x1 + 2x2 + x3 = 11,

x1 + 2x2 = 3,

2x1 + x2 + 4x3 = 16.

 

Составим матрицу [B|u]:

Первую строку [A|v] получим из первой строки [B|u] делением на

:   .

 

Вторая строка перед нормализацией выглядит так:

 

 

Фактор нормализации:   и, в результате, вторая строка [A|v] становится равной

 

 

Третья строка перед нормализацией:

 

 

Итого имеем матрицу [A|v]:

 

 

Выбрав в качестве начального вектор r0 = [1, 1, 1]T, получим

 

Затем

 

И, наконец, получаем решение системы:

 

 

 

 

 

ЧИТАТЬ  ФРАГМЕНТ  КНИГИ  НА  АНГЛИЙСКОМ