B.Carnahan, H.A.Luther, J.O.Wilkes. Applied Numerical Methods. Метод Качмажа

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

Applied Numerical Methods

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

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

1969

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

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

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

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

5.1 Введение

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

f₁(x₁, … , x_n) = 0,

f₂(x₁, … , x_n) = 0,

…………………...

……………………

f_n(x₁, … , x_n) = 0.

(5.1)

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

b₁₁x₁ + b₁₂x₂ + … + b_1nx_n = u₁ ,

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

b₁₁x₁ + b₁₂x₂ + … + b_1nx_n = u₁ ,

(5.2)

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

Bx = u , (5.3)

где B – матрица коэффициентов, u = [u₁, u₂, … , u_n]^T – вектор правой части системы, x = [x₁, x₂, … , x_n]^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.

Мы видим, что умножение слева на элементарные матрицы приводит к следующим преобразованиям матрицы 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, обладающая следующим свойством: AA^T = E, где A^T – транспонированная матрица A (а следовательно, и A^TA = E).

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

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

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

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

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

proj_ba = 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}$

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

В результате выполнения процесса Грама — Шмидта получаем систему ортогональных векторов b₁, …, b_n, либо систему ортонормированных векторов e₁, …, e_n .

Вычисление b₁, …, b_n носит название ортогонализации Грама — Шмидта, а e₁, …, e_n — ортонормализации Грама — Шмидта.

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

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

b₁₁x₁ + b₁₂x₂ + … + b_1nx_n = u₁ ,

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

b₁₁x₁ + b₁₂x₂ + … + b_1nx_n = u₁ ,

(5.2)

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

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

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

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

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

r_j = r_j_-1 – [(α_j, r_j_-1) - ν_j] α_j , 1 ≤ j ≤ n (5.10)

где A = (a_ij) и α_j = [a_j₁, a_j₂, … , a_jn]^T.

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

Br_n = u. (5.11)

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

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

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

(α_j, r_j) = ν_j ,

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

(α_i, r_j) = (α_i, r_j-1), (i ≠ j)

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

Вернемся к решению системы (5.2). Пусть B* = [β₁, β₂, … , β_n], где β_i = [b_i₁, b_i₂, … , b_in]^T. Система, эквивалентная (5.2) и имеющая свойства (5.9), строится на основе линейно независимых векторов β_i традиционным способом, а именно – с помощью процедуры ортогонализации Грама-Шмидта. Пусть

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

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

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

или

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

(β_j, x) = u_j.

На практике расчет может проводиться последовательным движением по массиву [B|u] с одновременным вычислением соответствующих строк массива [A|v]. Если необходимо уточнить вектор u после того, как матрица A уже была построена, нужно записать n чисел и (n² - 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} и {r_j} можно рассчитывать последовательно по шагам, поэтому метод можно считать n-шаговым.