Московский физико-технический институт
Опубликован: 25.10.2007 | Доступ: свободный | Студентов: 3790 / 1086 | Оценка: 4.50 / 4.33 | Длительность: 24:00:00
ISBN: 978-5-9556-0065-9
Специальности: Программист, Математик
Лекция 3:

Численное решение систем линейных алгебраических уравнений

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).

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

$ 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(n2) .

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

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 и исключим u1 из всех уравнений, начиная со второго, для чего ко второму уравнению прибавим первое, умноженное на -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), исключаем u2 из последних (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. Продолжая алгоритм, т.е. исключая ui (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)

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

\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*}

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

$ u_2  = \frac{1 + 4 \cdot 10^{- 7}}{1 + 2 \cdot 10^{- 7}} $.
После вычислений с той же точностью имеем: u2 = 1,000000, u1 = 2,000000, что является правильным решением (с заданным количеством значащих цифр).

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

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

|{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.

Решив последнюю систему, получаем \mathbf{\varepsilon}^1, после чего уточняем решение:

{\mathbf{u}}^2 = {\mathbf{u}}^1 + {\mathbf{\varepsilon }}^1.

Эту процедуру можно продолжить.