НОУ ИНТУИТ | Введение в вычислительную математику. Лекция 3: Численное решение систем линейных алгебраических уравнений

Учитесь и получайте официальные документы БЕСПЛАТНО. Вы можете поддержать наш проект.

Регистрация Вход

Твой путь к знаниям!

Опубликован: 25.10.2007 | Уровень: специалист | Доступ: платный | ВУЗ: Московский физико-технический институт

|

Вам нравится? Нравится 54 студентам

| Поделиться |

Поддержать программу

2.4. Прямые методы решения СЛАУ

Трудность численного решения рассматриваемых СЛАУ определяется видом матрицы $\mathbf{A}.$ Легко получается решение системы с диагональной матрицей, в этом случае система распадается на n линейных уравнений, каждое из которых содержит лишь одну неизвестную величину. Для диагональной системы очевидны явные формулы

$u_{k}= f_{k}/ a_{kk}, k = 1 \div n.$

В случае треугольной матрицы

$\mathbf{A}= \left( \begin{array}{cccc} a_{11} & a_{12} & \ldots & a_{1n} \\ 0 & a_{22} & \ldots & a_{2n} \\ \ldots & \ldots & \ldots & \ldots \\ 0 & 0 & \ldots & a_{nn} \end{array} \right)$

из последнего уравнения получаем $u_n = f_n /a_{nn}, (a_{ii} \ne 0 \mbox{, т.к. } \Delta =\det \mathbf{A} \ne 0).$

Решая систему линейных уравнений с треугольной матрицей "снизу вверх", для u_k имеем

$u_k= \frac{1}{a_{kk}}(f_k- a_{kn}u_n - a_{k,n- 1}u_{n - 1} - \ldots - a_{k,k + 1}u_{k + 1}), \\ \mbox{или }u_k= a_{kk}^{- 1}(f_k - \sum\limits_{j = k + 1}^n{a_{kj}u_j}),\quad k = n - 1, n - 2, \ldots , 1.$

Можно оценить количество арифметических действий, затрачиваемых на решение такой системы. Оно составляет O(n²) .

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

2.4.1. Метод исключения Гаусса

Рассматривается система уравнений

$\left\{ \begin{array}{ccc} {a_{11}u_1 + a_{12}u_2 + & \ldots & + a_{1n}u_n = f_1 ,} \\ {a_{21}u_1 + a_{22}u_2 + & \ldots & + a_{2n}u_n = f_2 ,} \\ & \ldots & \\ {a_{n1}u_1 + a_{n2}u_2 + & \ldots & + a_{nn}u_n = f_n .} \\ \end{array} \right.$

( 2.10)

Прямой ход метода Гаусса состоит в следующем. Положим, что $a_{11} \ne 0$ и исключим u₁ из всех уравнений, начиная со второго, для чего ко второму уравнению прибавим первое, умноженное на $-a_{21}/a_{11} = - \eta _{21}$ , к третьему прибавим первое, умноженное на $-a_{31}/a_{11} = - \eta _{31}$ и т.д. После этих преобразований получим эквивалентную систему:

$\left\{ \begin{array}{l} {a_{11}u_1 + a_{12}u_2 + \ldots + a_{1n}u_n = f_1 ,} \\ {a_{22}^1{u_2} + \ldots + a_{2n}^1{u_n} = f_2 ^1 ,} \\ \ldots \\ {a_{n2}^1{u_2} + \ldots + a_{nn}^1{u_n} = f_n^1 ,} \\ \end{array} \right.$

( 2.11)

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

$a_{ij}^1 = a_{ij}- \eta_{i1}a_{1j};\quad f_i^1 = f_i - \eta_{i1}f_1;\quad i,j = 2, \ldots , n.$

Теперь положим $a_{22}^1 \ne 0.$ Аналогично, вычислив множители второго шага $- a_{i2}^1/a_{22}^1 = -\eta_{i2}\quad (i = 3, \ldots , n)$ , исключаем u₂ из последних (n - 2) уравнений системы (2.17). В результате преобразований получим новую эквивалентную систему уравнений

$\left\{ \begin{array}{l} {a_{11}u_1 + a_{12}u_2 + a_{13}u_3 + \ldots + a_{1n}u_n = f_1} \\ {a_{22}^1 u_2 + a_{23}^1 u_3 + \ldots + a_{2n}^1 u_n = f_2^1} \\ {a_{33}^2 u_3 + \ldots + a_{3n}^2 u_n = f_3^2} \\ \ldots \\ {a_{n3}^2 u_3 + \ldots + a_{nn}^2 u_n = f_n^2} \\ \end{array} \right.$

в которой $a_{ij}^2 = a_{ij}^1 - \eta_{i2}a_{2j}^1;\quad f_i^2 = f_i^1 - \eta_{i2}f_2^1;\quad i,j = 3, \ldots , n.$ Продолжая алгоритм, т.е. исключая u_i (i = k + 1, ..., n), приходим на n - 1 шаге к системе с треугольной матрицей

$\left\{ \begin{array}{l} a_{11}u_1 + a_{12}u_2 + a_{13}u_3 + \ldots + a_{1n}u_n = f_1 \\ a_{22}^1 u_2 + a_{23}^1 u_3 + \ldots + a_{2n}^1 u_n = f_2^1 \\ a_{33}^2 u_3 + \ldots + a_{3n}^2 u_n = f_3^2 \\ \ldots \\ a_{nn}^{(n - 1)}u_n = f_n^{(n - 1)}. \\ \end{array} \right.$

( 2.12)

Обратный ход метода Гаусса позволяет определить решение системы линейных уравнений. Из последнего уравнения системы находим u_n ; подставляем это значение в предпоследнее уравнение, получим u_n-1. Поступая так и далее, последовательно находим u_n-2, u_n-3, ..., u₁. Вычисления компонент вектора решения проводятся по формулам

$\begin{gather*} u_n = f_n^{(n - 1)}/a_{nn}^{(n - 1)}, \\ \ldots \\ u_k = \frac{1}{a_{kk}^{(k - 1)}}(f_k^{(k - 1)}- a_{k,k + 1}^{(k - 1)}u_{k + 1}- \ldots - a_{kn}^{(k - 1)}u_n ),\quad k = n - 1, n - 2, \ldots , 1, \\ \ldots \\ u_2 = \frac{1}{a_{22}^1}(f_2^1 - a_{23}^1 u_3 - \ldots - a_{2n}^1 u_n ), \\ u_1 = \frac{1}{a_{11}}(f_1 - a_{12}u_2 - \ldots - a_{1n}u_n ). \end{gather*}$

Этот алгоритм прост и легко реализуем при условии, что $a_{11} \ne 0$ , $a_{22} \ne 0$ и т.д. Количество арифметических действий прямого хода $\approx 2/3n^{3}$ , обратного $\approx n^{2}.$ Это уже приемлемая для современных компьютеров величина.

Рассмотрим метод Гаусса с позиции операций с матрицами. Пусть ${\mathbf{A}}_1$ — матрица системы после исключения первого неизвестного

${\mathbf{A_1}}_1 = \left( \begin{array}{ccccc} a_{11} & a_{12} & a_{13} & \ldots & a_{1n}\\ 0 & a^1_{22} & a^1_{23} & \ldots & a^1_{2n}\\ 0 & a^1_{32} & a^1_{33} & \ldots & a^1_{3n}\\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 0 & a^1_{n2} & a^1_{n3} & \ldots & a^1_nn \\ \end{array} \right),\quad {\mathbf{f}}_1 = {\{f_1, f_2^1, \ldots , f_n^1\}}^T.$

Введем новую матрицу

${\mathbf{N}}_1 = \left(\begin{array}{ccccc} 1 & 0 & 0 & \ldots & 0 \\ - \eta_{21} & 1 & 0 & \ldots & 0 \\ - \eta_{31} & 0 & 1 & \ldots & 0 \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ - \eta_{n1} & 0 & 0 & \ldots & 1 \\ \end{array} \right).$

Очевидно, ${\mathbf{A}}_1={\mathbf{N}}_1\mathbf{A}$ , ${\mathbf{f}}_1={\mathbf{N}}_1\mathbf{f}.$ Аналогично, после второго шага система приводится к виду ${\mathbf{A}}_2\mathbf{u}={\mathbf{f}}_2$ , где ${\mathbf{A}}_2={\mathbf{N}}_2{\mathbf{A}}_1$ , ${\mathbf{f}}_2={\mathbf{N}}_2{\mathbf{f}}_1$ ,

${\mathbf{A}}_1 = \left( \begin{array}{ccccc} a_{11} & a_{12} & a_{13} & \ldots & a_{1n} \\ 0 & a_{22}^1 & a_{23}^1 & \ldots & a_{2n}^1 \\ 0 & 0 & a_{33}^2 & \ldots & a_{3n}^2 \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 0 & 0 & a_{n3}^2 & \ldots & a_{nn}^2 \\ \end{array} \right), \quad {\mathbf{N}}_2 = \left( \begin{array}{ccccc} 1 & 0 & 0 & \ldots & 0 \\ 0 & 1 & 0 & \ldots & 0 \\ 0 & - \eta_{32} & 1 & \ldots & 0 \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 0 & - \eta_{n2} & 0 & \ldots & 1 \\ \end{array} \right), \\ {{\mathbf{f}}_2} = \left\{{f_1, f_2^1, f_3^2, \ldots , f_n^2}\right\}^T.$

После n - 1 шага получим ${\mathbf{A}}_{n - 1}\mathbf{u} = {\mathbf{f}}_{n - 1}, {\mathbf{A}}_{n - 1} = {\mathbf{N}}_{n - 1} \cdot {\mathbf{A}}_{n - 2}, {\mathbf{f}}_{n - 1} = {\mathbf{N}}_{n - 1}{\mathbf{f}}_{n - 2}$ ,

${\mathbf{A}}_{n - 1}= \left( \begin{array}{ccccc} a_{11} & a_{12} & a_{13} & \ldots & a_{1n} \\ 0 & a_{22}^1 & a_{23}^1 & \ldots & a_{2n}^1 \\ 0 & 0 & a_{33}^2 & \ldots & a_{3n}^2 \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 0 & 0 & 0 & \ldots & a_{nn}^{(n - 1)} \\ \end{array} \right), \\ {\mathbf{N}}_{n - 1}= \left( \begin{array}{ccccc} 1 & 0 & \ldots & 0 & 0 \\ 0 & 1 & \ldots & 0 & 0 \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 0 & 0 & \ldots & 1 & 0 \\ 0 & 0 & \ldots & - \eta_{n,n}- 1} & 1 \\ \end{array} \right),\quad {\mathbf{f}}_{n - 1}= \{f_1, f_2^1, f_3^2, \ldots , f_n^{n - 1}\}^T.$

В итоге получаются матрица и вектор ${\mathbf{A}}_{n - 1} = {\mathbf{N}}_{n - 1} \ldots {\mathbf{N}}_2{\mathbf{N}}_1\mathbf{A}$ , ${\mathbf{f}}_{(n - 1)} = {\mathbf{N}}_{n - 1} \ldots {\mathbf{N}}_2{\mathbf{N}}_1\mathbf{f}$ , откуда $\mathbf{A}={\mathbf{N}}_1^{- 1}{\mathbf{N}}_2^{- 1} \ldots {\mathbf{N}}_{n - 1}^{- 1} \cdot {\mathbf{A}}_{n - 1}.$ При этом

${\mathbf{N}}_1^{- 1}= \left( \begin{array}{ccccc} 1 & 0 & 0 & \ldots & 0 \\ \eta_{21} & 1 & 0 & \ldots & 0 \\ \eta_{31} & 0 & 1 & \ldots & 0 \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ \eta_{n1} & 0 & 0 & \ldots & 1 \\ \end{array} \right),\quad {\mathbf{N}}_2^{- 1}= \left( \begin{array}{ccccc} 1 & 0 & 0 & \ldots & 0 \\ 0 & 1 & 0 & \ldots & 0 \\ 0 & \eta_{32} & 1 & \ldots & 0 \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 0 & \eta_{n2} & 0 & \ldots & 1 \\ \end{array} \right), \\ {\mathbf{N}}_{n - 1}^{- 1}= \left( \begin{array}{ccccc} 1 & 0 & \ldots & 0 & 0 \\ 0 & 1 & \ldots & 0 & 0 \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 0 & 0 & \ldots & 1 & 0 \\ 0 & 0 & \ldots & \eta_{n,n}- 1} & 1 \\ \end{array} \right).$

После введения обозначений $\mathbf{U} = {\mathbf{A}}_{n - 1}$ , $\mathbf{L} = {\mathbf{N}}_1^{- 1}{\mathbf{N}}_2^{- 1} \ldots {\mathbf{N}}_{n - 1}^{- 1}$ , где

$\mathbf{L} = \left( \begin{array}{ccccc} 1 & 0 & 0 & \ldots & 0 \\ \eta_{21} & 1 & 0 & \ldots & 0 \\ \eta_{31} & \eta_{32} & 1 & \ldots & 0 \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ \eta_{n1} & \eta_{n2} & \eta_{n3} & \ldots & 1 \\ \end{array} \right),$

получим $\mathbf{A} = \mathbf{LU}.$

Это представление матрицы $\mathbf{A}$ называется LU-разложением (на произведение нижней и верхней треугольных матриц $\mathbf{L}$ и $\mathbf{U}$ ). Прямой ход метода Гаусса можно рассматривать как один из вариантов представления матрицы в виде произведения двух треугольных матриц, или LU-разложения. Его можно провести и другими способами.

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

$\begin{gather*} -10^{- 7}u_1 + u_2 = 1, \\ u_1 + 2u_2 = 4. \end{gather*}$

Исключая u₁ из первого уравнения и подставляя во второе, получим u₂ = (10⁷ + 4)/(10⁷ + 2). После вычислений с семью значащими цифрами получаем u₁ = 0,000000, u₂ = 1,000000, что неверно (см. второе уравнение). Теперь исключим u₁ из второго уравнения и подставим в первое. При этом получим

$$ u_2 = \frac{1 + 4 \cdot 10^{- 7}}{1 + 2 \cdot 10^{- 7}} $.$

После вычислений с той же точностью имеем: u₂ = 1,000000, u₁ = 2,000000, что является правильным решением (с заданным количеством значащих цифр).

В реальных вычислениях используются методы с выбором главного (или ведущего ) элемента. Выбор главного элемента по столбцам реализуется следующим образом: перед исключением u₁ отыскивается $\max\limits_i|{a_{i1}}|.$ Пусть максимум достигается при i = k. В этом случае меняются местами первое и k уравнения (или в матрице меняются местами две строки) и реализуется процедура исключения.

Затем отыскивается $\max\limits_i|{a_{i2}^1}|$ , и процедура поиска главного элемента в столбцах повторяется. Так же реализуется выбор главного элемента по строкам: перед исключением u₁ отыскивается $\max\limits_j|{a_{kj}}|.$ Если максимум достигается при i = k, то у u₁ и u_k меняются номера, то есть максимальный элемент из коэффициентов первого уравнения окажется на месте a₁₁, и т.д. Наиболее эффективным является метод Гаусса с выбором главного элемента по всей матрице. Во многих методах важным является условие диагонального преобладания

$|{a_{ii}}| \ge \sum\limits_{\substack{j = 1 \\ j \ne i}}^n{|{a_{ij}}|}$

для i = 1, ..., n, при выполнении которого проблемы, появляющиеся в методе Гаусса, не возникают. Если для всех строк матрицы выполняются строгие неравенства, то говорят о строгом диагональном преобладании.

Полученное решение можно улучшить следующим образом. Пусть ${\mathbf{r}}^1 = \mathbf{f} - {\mathbf{Au}}^1$ есть невязка, допущенная при решении рассматриваемой системы ( ${\mathbf{u}}^1$ — полученное численное решение) за счет ошибки округлений. Очевидно, что погрешность ${\mathbf{\varepsilon}}^1 = \mathbf{u} - {\mathbf{u}}^1$ удовлетворяет СЛАУ ${\mathbf{A\varepsilon}}^1 = {\mathbf{r}}^1$ , так как ${\mathbf{A\varepsilon}}^1 = \mathbf{Au} - {\mathbf{Au}}^1 = \mathbf{f}- {\mathbf{Au}}^1.$