Системы уравнений | Линейная алгебра

Все курсы > Линейная алгебра > Занятие 4

Начнем изучать системы линейных уравнений.

Содержание занятия

Обратная матрица

Обратимся к следующей системе.

$$ \begin{bmatrix} 1 & 1 & 3 \\ 1 & 2 & 4 \\ 1 & 1 & 2 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 15 \\ 21 \\ 13 \end{bmatrix}$$

Перепишем ее в виде

$$ A \mathbf x = \mathbf b $$

Вспомним, что если существует обратная матрица (inverse matrix) $A^{-1}$, то при умножении на матрицу $A$ мы получаем единичную матрицу $I$, т.е. $A^{-1}A = I$. Умножим обе стороны уравнения на $A^{-1}$.

$$ A^{-1}A \mathbf x = A^{-1} \mathbf b $$

$$ I \mathbf x = A^{-1} \mathbf b $$

$$ \mathbf x = A^{-1} \mathbf b $$

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

$$ A \mathbf x = \mathbf b \rightarrow A^{-1} \mathbf b = \mathbf x $$

С точки зрения линейных преобразований матрица A перемещает вектор $ \mathbf x$ на вектор $ \mathbf b$. Обратная же матрица «возращает» вектор $ \mathbf b$ на вектор $ \mathbf x$.

Решим систему уравнений с помощью метода .inv() модуля linalg библиотеки Numpy.

# создадим матрицу коэффициентов при неизвестных

A = np.array([[1, 1, 3],

[1, 2, 4],

[1, 1, 2]])

# и вектор-столбец свободных коэффициентов

b = np.array([[15],

[21],

[13]])

# решим систему уравнений

x = np.linalg.inv(A).dot(b)

array([[5.],

[4.],

[2.]])

Метод Гаусса

Приведение матрицы к ступенчатому виду

Рассмотрим метод Гаусса (Gaussian Elimination) для решения систем линейных уравнений.

Метод заключается в том, чтобы с помощью ряда элементарных операций привести расширенную матрицу (augmented matrix), т.е. матрицу коэффициентов при неизвестных ($A$) + вектор свободных коэффициентов ($ \mathbf b $), к ступенчатому виду (row-echelon form, ref) путем последовательного исключения коэффициентов (elimination).

Получив ступенчатую матрицу, мы сможем последовательно найти все переменные системы, начиная с последней (back-substitution).

Приведем перечень трех возможных операций:

перестановка местами двух строк матрицы
умножение на число, отличное от нуля
прибавление одной строки к другой, умноженной на константу (вычитание предполагает умножение на $-1$ и сложение)

Приведем пример решения⧉ системы уравнения вручную и пример алгоритма⧉ на Питоне, выполняющего эти операции.

100

101

102

103

104

105

106

107

108

109

110

# function to get matrix content

def gaussianElimination(mat):

N = mat.shape[0]

# reduce matrix to r.e.f.

singular_flag = forwardElim(mat, N)

# if matrix is singular

if (singular_flag != -1):

print("Singular Matrix.")

# if the RHS of equation corresponding to

# zero row is 0, * system has infinitely

# many solutions, else inconsistent

if (mat[singular_flag][N]):

print("Inconsistent System.")

else:

print("May have infinitely many solutions.")

return

# get solution to system and print it using

# backward substitution

backSub(mat, N)

# function for elementary operation of swapping two rows

def swap_row(mat, i, j, N):

for k in range(N + 1):

temp = mat[i][k]

mat[i][k] = mat[j][k]

mat[j][k] = temp

# function to reduce matrix to r.e.f.

def forwardElim(mat, N):

for k in range(N):

# initialize maximum value and index for pivot

i_max = k

v_max = mat[i_max][k]

# find greater amplitude for pivot if any

for i in range(k + 1, N):

if (abs(mat[i][k]) > v_max):

v_max = mat[i][k]

i_max = i

# if a principal diagonal element is zero,

# it denotes that matrix is singular, and

# will lead to a division-by-zero later

if not mat[k][i_max]:

return k # matrix is singular

# swap the greatest value row with current row

if (i_max != k):

swap_row(mat, k, i_max, N)

for i in range(k + 1, N):

# factor f to set current row kth element to 0,

# and subsequently remaining kth column to 0

f = mat[i][k]/mat[k][k]

# subtract fth multiple of corresponding kth

# row element

for j in range(k + 1, N + 1):

mat[i][j] -= mat[k][j]*f

# filling lower triangular matrix with zeros

mat[i][k] = 0

# print(mat); // for matrix state

return -1

# function to calculate the values of the unknowns

def backSub(mat, N):

x = [None for _ in range(N)] # an array to store solution

# start calculating from last equation up to the first

for i in range(N-1, -1, -1):

# start with the RHS of the equation

x[i] = mat[i][N]

# initialize j to i+1 since matrix is upper

# triangular

for j in range(i + 1, N):

# subtract all the lhs values

# except the coefficient of the variable

# whose value is being calculated

x[i] -= mat[i][j] * x[j]

# divide the RHS by the coefficient of the

# unknown being calculated

x[i] = (x[i] / mat[i][i])

print("Solution for the system:")

for i in range(N):

print("{:.2f}".format(x[i]))

print()

print('Row-echelon form:')

print(mat)

Применим этот алгоритм для решения приведенной в начале системы. Одновременно выведем получившуюся расширенную матрицу в ступенчатом виде.

# расширенная матрица

A = np.array([[1.0, 1.0, 3.0, 15.0],

[1.0, 2.0, 4.0, 21.0],

[1.0, 1.0, 2.0, 13.0]])

gaussianElimination(A)

Solution for the system:

5.00

4.00

2.00

Row-echelon form:

[[ 1. 1. 3. 15.]

[ 0. 1. 1. 6.]

[ 0. 0. -1. -2.]]

Важно отметить, что, не зная обратной матрицы, мы нашли решение лишь для конкретного $\mathbf b$.

Упрошенный ступенчатый вид матрицы

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

$$ \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 5 \\ 4 \\ 2 \end{bmatrix} $$

В этом можно убедиться, выполнив:

$ R_3 \times (-1)$
$ -3R_3 + R_1 $
$ -1R_3 + R_2 $
$ -1R_2 + R_1 $

Если матрица приведена к ступенчатому виду и кроме этого:

Первый ненулевой коэффициент каждой ненулевой строки равен единице (его называют разрешающим или ведущим коэффициентом, по-английски pivot, leading coefficient)
Столбец, в котором есть разрешающий элемент (единица), содержит только нулевые значения

то говорят, что матрица приведена к упрощенному ступенчатому виду (reduced row echelon form, rref). Ниже для сравнения приведены две матрицы, одна в ступенчатом виде (ref), вторая в упрощенном ступенчатом виде (rref).

сравнение ступенчатого и упрощенного ступенчатого вида матрицы

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

$$ rref(A) = I \quad \iff \quad A \text{ is non-singular} $$

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

Метод Гаусса-Жордана

Приведение матрицы к виду выше (когда слева остается единичная матрица) можно использовать для нахождения обратной матрицы методом Гаусса-Жордана (Gauss–Jordan elimination). Суть метода заключается в том, чтобы (пример взять отсюда⧉):

Во-первых, расширить исходную матрицу с помощью единичной матрицы.

метод Гаусса-Жордана: расширение исходной матрицы

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

метод Гаусса-Жордана: элементарные преобразования

Получившаяся справа матрица $B$ будет обратной матрице $A$ ($B = A^{-1}$).

Теперь с обратной матрицей мы можем найти $ \mathbf x$ вне зависимости от значений $ \mathbf b$.

Идея такого преобразования заключается в следующем. Приведенные выше операции с матрицей можно выразить через умножение трансформационной матрицы $Е$ на матрицу $A$. Например, возьмем вектор $ \mathbf b$.

Для того чтобы превратить 8 в 4 можно применить следующую трансформационную матрицу.

Эта матрица аналогична операции $ -2R_1 + R_2 $. В целом все операции по приведению матрицы $A$ к виду $I$ можно записать в виде матрицы $Е$. Тогда получаем, что $EA = I$. При этом по определению $E = A^{-1}$.

Найдем обратную матрицу с помощью написанного вручную алгоритма, а затем проверим решение через метод np.linalg.inv().

# создадим расширенную матрицу

A = np.array([[1, 1, 3],

[1, 2, 4],

[1, 1, 2]])

n = len(A)

I = np.identity(n)

M = np.concatenate((A, I), axis = 1)

# iterate over matrix rows

for i in range(0, n):

# select pivot value

pivot = M[i][i]

# extract row

row = M[i]

# get 1 along the diagonal

M[i] = row / pivot

# iterate over all rows except pivot to get augmented matrix into reduced row echelon form

for j in [k for k in range(0, n) if k != i]:

# subtract current row from remaining rows

M[j] = M[j] - M[i] * M[j][i]

# extract inverse matrix

inverse = M[:, n:]

inverse

array([[ 0., -1., 2.],

[-2., 1., 1.],

[ 1., -0., -1.]])

1	np.linalg.inv(A)

array([[ 0., -1., 2.],

[-2., 1., 1.],

[ 1., -0., -1.]])

1	np.dot(A, np.linalg.inv(A))

array([[1., 0., 0.],

[0., 1., 0.],

[0., 0., 1.]])

Подведем итог

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

Перейдем к изучению определителя матрицы.