10. Вычисление обратной матрицы методом Гаусса-Жордана

Алиас: Гаус-Жордан.

Введение

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

На практике для получения обратной матрицы используются методы, основанные на решении системы линейных алгебраических уравнений.

В лабораторной работе рассматривается вычисление обратной матрицы методом Гаусса-Жордана. Гаусс предложил метод решения систем линейных алгебраических уравнений путем последовательного исключения переменных, а затем Жордан использовал его для нахождения обратной матрицы, добавив преобразование единичной матрицы.

Цель работы

• разработка программы, реализующей алгоритм вычисления обратной матрицы методом Гаусса-Жордана;
• использование разработанных функций для решения систем линейных уравнений.

Алгоритм вычисления обратной матрицы

Прежде чем приступать к программированию, разработаем численный пример получения обратной матрицы методом Гаусса-Жордана. Это позволит лучше понять алгоритм метода и использовать этот пример в качестве контрольного при разработке алгоритма программы. Предположим, имеется матрица

$$A= \begin{bmatrix} 2 & 5 & 7 \\ 3 & 9 & 15 \\ 5 & 16 & 20 \end{bmatrix}.$$

Метод Гаусса-Жордана предполагает нахождение обратной матрицы путем решения следующей системы уравнений

$$A \begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} a \\ b \\ c \end{bmatrix} \tag{1}$$

относительно $x$, $y$ и $z$ с помощью обычной процедуры последовательного исключения (здесь символы $x$, $y$, $z$ и $a$, $b$, с характеризуют скалярные величины).

Умножение обеих частей уравнения слева на $A^{-1}$ дает решение в виде

$$\begin{bmatrix} x \\ y \\ z \end{bmatrix} = A^{-1} \begin{bmatrix} a \\ b \\ c \end{bmatrix} \tag{2}$$

если $A^{-1}$ существует.

Следовательно, если решение системы $(1)$ относительно $x$, $y$ и $z$ выражено в явном виде через $a$, $b$ и $c$ и записано в матричной форме, то приводящая к решению матрица равна $A^{-1}$. Воспроизведем полностью этот процесс.

Допустим, что уравнения $(1)$ имеют следующий вид:

$$2x + 5y + 7z = a,\tag{3}$$ $$3x + 9y + 15z = b,\tag{4}$$ $$5x + 16y + 20z = c.\tag{5}$$

Деление $(3)$ на $2$ дает.

$$x + 2.5y + 3.5z = 0.5a.\tag{6}$$

Умножим обе части уравнения $(6)$ на три и вычтем их из $(4)$, а затем, умножив обе части уравнения $(6)$ на пять, вычтем их из $(5)$. Тогда

$$1.5y + 4.5z = b - 1.5a\tag{7}$$

и

$$3.5y + 2.5z = c - 2.5a.\tag{8}$$

Разделив $(7)$ на $1.5$ и умножив $(8)$ на $2$, получим

$$y + 3z = b/1.5 - a,\tag{9}$$ $$7y + 5z = 2c - 5a.\tag{10}$$

Вычитая $(9)$, умноженное на $7$, из $(10)$, мы получим

$$-16z = 2c - 5a - 7 (b/1.5 - a).$$

Следовательно,

$$z = (2a - 7b/1,5 + 2c)/(- 16) = (-3a + 7b - 3c)/24. \tag{11}$$

Подстановка $z$ в $(9)$ приводит к следующему результату:

$$y = (-15a - 5b + 9c)/24,\tag{12}$$

а подстановка $y$ и $z$ в $(6)$ дает

$$x = (60a - 12b - 12c)/24.\tag{13}$$

Последние три результата могут быть представлены в матричном виде следующим образом:

$$\begin{bmatrix} x \\ y \\ z \end{bmatrix} = \frac{1}{24}\begin{bmatrix} 60 & -12 & -12 \\ -15 & -5 & 9 \\ -3 & 7 & -3 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix} . \tag{14}$$

Сравнивая этот результат с $(2)$, мы приходим к заключению, что

$$A^{-1} = \frac{1}{24} \begin{bmatrix} 60 & -12 & -12 \\ -15 & -5 & 9 \\ -3 & 7 & -3 \end{bmatrix}.$$

Умножив обе части этого соотношения на $A$, можно убедиться в справедливости полученного результата. Эта процедура, естественно, очень утомительна для человека, но для компьютера она не представляется сложной: отдельные шаги расчетов по существу одни и те же как для больших, так и для маленьких матриц (для больших их просто больше).

Так как компьютер работает с числами (а не с переменными $a$, $b$, $c$), то для того чтобы решение системы $(1)$ относительно $x$, $y$ и $z$ выразить в явном виде через $a$, $b$ и $c$ и записать в матричной форме, умножим правую часть системы уравнений $(1)$ на единичную матрицу и будем выполнять над ней те же преобразования, что и над исходной матрицей.

Алгоритм преобразования состоит из двух шагов:

1. Прямой ход. В результате исходная матрица преобразуется в верхнюю треугольную матрицу.
2. Обратный ход. В результате исходная матрица преобразуется в единичную, а введенная единичная матрица – в матрицу, обратную исходной, то есть система уравнений $(1)$ преобразуется в систему уравнений $(2)$.

Выполним с помощью этого алгоритма рассмотренный выше пример:

$$\begin{bmatrix} 2 & 5 & 7 \\ 3 & 9 & 15 \\ 5 & 16 & 20 \end{bmatrix} \begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix}$$

Шаг 1. Прямой ход

Исключаем $x$ из нижележащих строк:

$$\begin{bmatrix} 1 & 2.5 & 3.5 \\ 0 & 1.5 & 4.5 \\ 0 & 3.5 & 2.5 \end{bmatrix} \quad \quad \begin{bmatrix} 0.5 & 0 & 0 \\ -1.5 & 1 & 0 \\ -2.5 & 0 & 1 \end{bmatrix}.$$

Исключаем $y$:

$$\begin{bmatrix} 1 & 2.5 & 3.5 \\ 0 & 1 & 3 \\ 0 & 0 & -16 \end{bmatrix} \quad \quad \begin{bmatrix} 0.5 & 0 & 0 \\ -1 & 2/3 & 0 \\ 2 & -14/3 & 2 \end{bmatrix}.$$

Устанавливаем коэффициент при $z$ равным $1$:

$$\begin{bmatrix} 1 & 2.5 & 3.5 \\ 0 & 1 & 3 \\ 0 & 0 & 1 \end{bmatrix} \quad \quad \begin{bmatrix} 0.5 & 0 & 0 \\ -1 & 2/3 & 0 \\ -1/8 & 7/24 & -1/8 \end{bmatrix}.$$

Шаг 2. Обратный ход

Исключаем $z$ из всех вышележащих строк:

$$\begin{bmatrix} 1 & 2.5 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad \quad \begin{bmatrix} 15/16 & -49/48 & 7/16 \\ -5/8 & -5/24 & 3/8 \\ -1/8 & 7/24 & -1/8 \end{bmatrix}.$$

Исключаем $y$ из всех вышележащих строк:

$$\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad \quad \begin{bmatrix} 5/2 & -1/2 & -1/2 \\ -5/8 & -5/24 & 3/8 \\ -1/8 & 7/24 & -1/8 \end{bmatrix}.$$

Представим результат в матричном виде:

$$\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x \\ y \\ z \end{bmatrix} = \frac{1}{24} \begin{bmatrix} 60 & -12 & -12 \\ -15 & -5 & 9 \\ -3 & 7 & -3 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix}.$$

Сравнив результат с $(14)$, убеждаемся, что все вычисления сделаны правильно.

Точность вычисления обратной матрицы

Основной источник погрешностей при выполнении вычислений на ЭВМ – это ошибки округления. Ошибки округления всегда возникают в последнем значащем разряде. Однако накопление ошибок в связи с выполнением большого количества вычислений может повлиять на последние два и более разрядов.

Для уменьшения погрешностей округления нужно использовать более "длинные" типы данных для представления коэффициентов, например, double вместо float.

Влияние погрешностей округления на результат вычислений зависит от используемых методов расчета. Наиболее часто встречающиеся источники увеличения влияния погрешностей округления на результат – деление на малые величины и сложение больших величин с малыми величинами. Полностью предупредить подобные ошибки нельзя, однако в каждом конкретном случае нужно использовать методы расчета, которые сводят к минимуму влияние ошибок. Например, если нужно складывать большое количество малых величин с большой величиной, то нужно сложить сначала все малые величины, а затем их сумму, которая уже не будет такой малой, сложить с большой величиной. Если при сложении последовательно прибавлять малые величины к большой, то при выравнивании порядков малые величины будут превращаться в ноль и не повлияют на результат сложения. Общее правило заключается в том, что надо так организовать процесс вычислений, чтобы, по возможности, все операции выполнялись бы над величинами одного порядка.

Вообще говоря, ЭВМ не делает ошибок, но она выполняет ошибочные программы точно так же, как и правильные. Поэтому оценку точности результатов вычислений должен делать программист. В случае получения обратной матрицы правильность результата вычислений можно оценить путем умножения обратной матрицы на исходную. В результате должна получиться единичная матрица.

Задание

Создать функцию для вычисления обратной матрицы по методу Гаусса-Жордана. Размер матрицы передавать в функцию в качестве параметра. Для упрощения алгоритма следует присоединить единичную матрицу справа к исходной и выполнять все преобразования над объединенной матрицей размером $N*2N$. Обратная матрица получится на месте единичной в столбцах $N...2N$, а на месте исходной матрицы в столбцах $0...(N-1)$ должна получиться единичная матрица.

Включить в алгоритм проверку на существование обратной матрицы. Для этого в в прямом ходе перед делением выполнить проверку на ноль элементов главной диагонали исходной матрицы. Если элемент равен $0$, то нужно поменять местами текущую строку с одной из нижележащих строк, в которой элемент в соответствующем столбце не равен $0$.

Если таких строк нет, то выдать сообщение: "Обратная матрица не существует".

Применить функцию для решения системы линейных алгебраических уравнений.

Шаги выполнения задания

Часть 1

1. Разработать схему алгоритма для вычисления обратной матрицы методом Гаусса-Жордана.
2. Разработать программу, реализующую метод Гаусса-Жордана. Для ускорения процесса разработки программы на этом шаге следует сконцентрироваться на реализации алгоритма и не использовать динамическую память для хранения матриц и не оформлять реализацию алгоритма в виде функции.
3. Отладить программу, используя контрольный пример.
4. Оформить алгоритм в виде функции для вычисления обратной матрицы по методу Гаусса-Жордана. Размер матрицы передавать в функцию в качестве параметра.

Часть 2

Выполнить тестирование программы. Проверку правильности результатов вычислений выполнять путем умножения полученной обратной матрицы на исходную (в результате должна получиться единичная матрица). Для реализации контроля разработайте функцию, реализующую произведение прямоугольных матриц. В процессе выполнения программы на экран должны выводиться следующие данные: исходная матрица, промежуточные результаты после прямого и обратного ходов, обратная матрица и результат умножения прямой матрицы на обратную. Для печати матрицы использовать функцию из лабораторной работы 8.

Доработать алгоритм для случая, когда на главной диагонали появляется ноль. Включить в алгоритм проверку на существование обратной матрицы.

Применить функцию вычисления обратной матрицы для решения системы линейных алгебраических уравнений.

Исследовать влияние коэффициентов исходной матрицы на точность вычисления обратной матрицы. В качестве оценки погрешности использовать значения "нулей" единичной матрицы, полученной в результате перемножения исходной матрицы на обратную.

Проверка задания

Подготовленная программа для решения задания проверяется вручную преподавателем (визуальный контроль).

Методический материал

1. Требования к отчету