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

Материал из MachineLearning.

(Различия между версиями)
Перейти к: навигация, поиск
(категория)
 
(15 промежуточных версий не показаны.)
Строка 25: Строка 25:
:: 1. Предполагаем, что <tex>a_{11} \neq 0</tex> . Тогда первое уравнение системы {{eqref|2}} делим на коэффициент <tex>a_{11}</tex>, в результате получаем уравнение
:: 1. Предполагаем, что <tex>a_{11} \neq 0</tex> . Тогда первое уравнение системы {{eqref|2}} делим на коэффициент <tex>a_{11}</tex>, в результате получаем уравнение
<center><tex>x_{1}+\sum_{j=2}^n a_{ij}^{1}x_{j} = b_{i}^{1}</tex></center>.
<center><tex>x_{1}+\sum_{j=2}^n a_{ij}^{1}x_{j} = b_{i}^{1}</tex></center>.
-
::Затем из каждого из остальных уравнений вычитается первое уравнение, умноженное на соответствующий коэффициент <tex>a_{i1}</tex>. В результате уравнения преобразуются к виду:
+
::Затем из каждого из оставшихся уравнений вычитается первое, умноженное на соответствующий коэффициент <tex>a_{i1}</tex>. В результате система преобразуются к виду:
<center><tex>\left( \begin{array}{ccc} 1 & a_{12}^{1} & \ldots & a_{1n}^{1}\\ \\ 0 & a_{22}^{1} & \ldots & a_{2n}^{1}\\ \\ \vdots & \vdots & \ddots & \vdots\\ \\ 0 & a_{n2}^{1} & \ldots & a_{nn}^{1} \end{array}\right) \left( \begin{array}{c}x_1 \\ \\ x_2 \\ \\ \vdots \\ \\ \\ x_n \end{array}\right) = \left( \begin{array}{c}b_1^1 \\ \\ b_2^1 \\ \\ \vdots \\ \\ b_n^1 \end{array}\right)</tex></center>
<center><tex>\left( \begin{array}{ccc} 1 & a_{12}^{1} & \ldots & a_{1n}^{1}\\ \\ 0 & a_{22}^{1} & \ldots & a_{2n}^{1}\\ \\ \vdots & \vdots & \ddots & \vdots\\ \\ 0 & a_{n2}^{1} & \ldots & a_{nn}^{1} \end{array}\right) \left( \begin{array}{c}x_1 \\ \\ x_2 \\ \\ \vdots \\ \\ \\ x_n \end{array}\right) = \left( \begin{array}{c}b_1^1 \\ \\ b_2^1 \\ \\ \vdots \\ \\ b_n^1 \end{array}\right)</tex></center>
:: 2. В предположении, что <tex>a_{22}^1 \neq 0</tex>, делим второе уравнение на коэффициент <tex>a_{22}^1</tex> и исключаем неизвестное <tex>x_2</tex> из всех последующих уравнений и т.д.
:: 2. В предположении, что <tex>a_{22}^1 \neq 0</tex>, делим второе уравнение на коэффициент <tex>a_{22}^1</tex> и исключаем неизвестное <tex>x_2</tex> из всех последующих уравнений и т.д.
Строка 38: Строка 38:
== Анализ метода ==
== Анализ метода ==
 +
 +
Данный метод относится к классу прямых методов решения системы уравнений, а это значит, что за конечное число шагов можно получить точное решение, при условии, что входные данные ( матрица <tex>A</tex> и правая часть уравнения - <tex>b</tex>) заданы точно и вычисление ведется без округлений.
 +
Для получения решения требуется <tex>\frac{m^3}{3}+ m^2 - \frac{m}{3}</tex> умножений и делений, то есть порядка <tex>O(\frac{m^3}{3})</tex> операций.
 +
 +
Условия, при которых метод выдает точное решение, на практике не выполнимы - неизбежны как ошибки входных данных, так и ошибки округления.
 +
Тогда встает вопрос: насколько точное решение можно получить, используя метод Гаусса, насколько метод корректен?
 +
Определим устойчивость решения относительно входных параметров. Наряду с исходной системой {{eqref|1}} рассмотрим возмущенную систему:
 +
<center><tex>(A + \delta A)(x + \delta x)=b + \delta b</tex></center>
 +
Пусть введена некоторая норма <tex>|| . ||</tex>. <tex>CondA=|| A || * || A^{-1}||</tex> - называется числом обусловленности матрицы <tex>A</tex>.
 +
 +
Возможны 3 случая:
 +
 +
1) <tex>\delta A = 0 , \quad \delta b \neq 0</tex> :
 +
:<tex>\frac{||\delta x||}{||x||} &\leq& Cond A \frac{||\delta b ||}{|| b ||}</tex>
 +
2) <tex>\delta A \neq 0 , \quad \delta b = 0</tex> :
 +
:<tex><tex>\frac{||\delta x||}{||x||} \leq Cond A \frac{||\delta A ||}{|| A ||}</tex>
 +
3) <tex>\delta A \neq 0 , \quad \delta b \neq 0</tex> :
 +
:<tex><tex>\frac{||\delta x||}{||x||} \leq \frac{Cond A}{1-Cond A\frac{||\delta A ||}{|| A ||}}(\frac{||\delta A ||}{|| A ||}+ \frac{||\delta b ||}{|| b ||})</tex>
 +
 +
Число обусловленности матрицы <tex>A</tex> всегда <tex>\geq 1</tex>. Если оно велико ( <tex>Cond A \approx 10^k k \geq 2</tex> ) , то говорят, что матрица <tex>A</tex> плохо обусловлена. В этом случае малые возмущения правых частей системы {{eqref|1}}, вызванные либо неточностью задания исходных данных, либо вызванные погрешностями вычисления, существенно влияют на решение системы. Грубо говоря, если погрешность правых частей <tex>10^{-l}</tex> , то погрешность решения будет <tex>10^{-l+k}</tex> .
 +
 +
Проиллюстрируем полученные результаты на следующем числовом примере:
 +
Дана система
 +
 +
:<tex>\left\{\begin{array}{lcr}5x_1+7x_2 &=& 12 \\ \\ 7x_1+10x_2 &=& 17 \\ \end{array} \right.</tex>
 +
 +
Она имеет решение <tex>(1, \quad 1)</tex>.
 +
 +
Теперь рассмотрим возмущенную систему:
 +
 +
:<tex>\left\{\begin{array}{lcr}5x_1+7x_2 &=& 12.075 \\ \\ 7x_1+10x_2 &=& 16.905 \\ \end{array} \right.</tex>
 +
 +
Решением такой системы будет вектор <tex>(2.415, \quad 0)</tex>.
 +
 +
При совсем малом возмущении правой части получили несоизмеримо большое возмущение решения.
 +
Объяснить такую "ненадежность" решения можно тем, что матрица <tex>A</tex> почти вырожденная: прямые, соответствующие двум уравнениям, почти совпадают, что видно на графике:
 +
 +
[[Изображение:Слау.jpg|center|Геометрическое представление системы двух линейных алгебраических уравнений, которая является почти вырожденной. Прямые, соответствующие двум уравнениям, почти совпадают.]]
 +
 +
 +
 +
Такой результат можно было предвидеть в силу плохой обусловленностью матрицы <tex>A = \left( \begin{array}{ccc} 5 & 7\\ \\ 7 & 10 \end{array}\right)</tex> :
 +
<tex>Cond A = 17^2</tex> <ref>Здесь и далее в числовых примерах будем рассматривать метрику <tex>||.||_1</tex> : <tex>||x||_1=\sum_{i=1}^n | x_i |, \quad || A ||_1 = \max_{j} \left( \sum_{i=1}^n | a_{ij} |\right)</tex></ref>
 +
 +
Вычисление <tex>Cond A </tex> является достаточно сложным, сравнимо с решением всей системы, поэтому для оценки пограшности применяются более грубые, но простые в реализации методы.
== Способы оценки ошибок ==
== Способы оценки ошибок ==
 +
 +
:1) <u>'''Контрольная сумма:'''</u> обычно применяется для предупреждения случайных погрешностей в процессе вычисления без помощи компьютеров.
 +
 +
Составляем контрольный столбец <tex>a_{n+1} = \left( a_{1,n+1}, \ldots, a_{n, n+1}\right)</tex>, состоящий из контрольных элементов системы:
 +
:<tex>a_{i, n+1} = \sum_{j=1}^n a_{i, j} + b_i</tex>
 +
 +
При преобразовании уравнений над контрольными элементами производятся те же операции, что и над свободными членами уравнеий. В результате этого контрольный элемент каждого нового уравнения должен равняться сумме коэффициентов этого уравнения. Большое расхождение между ними указывает на погрешности в вычислениях или на неустойчивость алгоритма вычислений по отношению к вычислительной погрешности.
 +
 +
:2) <u>'''Относительная погрешность известного решения'''</u> позволяет без существенных дополнительных затрат получить суждение о погрешности решения.
 +
 +
Задается некоторый ветор <tex>z</tex> с компонентами, имеющими по возможности тот же порядок и знак, что и компоненты искомого решения <ref>Часто из-за отсутствия доплнительной информации берут <tex>z=\left(1, \ldots, 1\right)</tex></ref>. Вычисляется вектор <tex>c=Az</tex>, и на ряду с исходной системой уравнения решается система <tex>Az=c</tex>.
 +
 +
Пусть <tex>x'</tex> и <tex>z'</tex> - реально получаемые решения этих систем. Суждение о погрешности <tex>x-x'</tex> искомого решения можно получить, основываясь на гипотезе: относительные погрешности при решении методом исключения систем с одной и той же матрицей и различными правыми частями, которыми являются соответственно величины <tex>\frac{||x-x'||}{||x'||}</tex> и <tex>\frac{||z-z'||}{||z'||}</tex>, отличаются не в очень большое число раз.
 +
 +
:3) <u>'''Изменение масштабов'''</u> - прием, применяющийся для получения представления о реальной величине погрешности, возникающей за счет округлений при вычислениях.
 +
 +
Наряду с исходной системой тем же методом решается система
 +
:<tex>\left( \alpha A\right)x' = \beta b</tex>, где <tex>\alpha</tex> и <tex>\beta</tex> - числа
 +
 +
Если бы не было погрешности округления, то выполнялось бы равенство для решений исходной и масштабированной систем: <tex>x=\alpha \beta^{-1}x'</tex>. Поэтому при <tex>\alpha</tex> и <tex>\beta</tex>, не являющихся степенями двойки, сравнение векторов <tex>x</tex> и <tex>\alpha \beta^{-1}x'</tex> дает представление о величине вычислительной погрешности <ref>Например, можно взять <tex>\alpha=\sqrt{2}, \quad \beta=\sqrt{3}</tex></ref>
== Улучшение метода исключения Гаусса==
== Улучшение метода исключения Гаусса==
 +
 +
Рассмотренные ниже модификации метода Гаусса позволяют уменьшить погрешность результата.
 +
=== Выбор главного элемента ===
=== Выбор главного элемента ===
 +
Основное увеличение ошибки в методе происходит во время прямого хода, когда ведущая <tex>k</tex>-я строка умножается на коэффициенты <tex>\frac{a_{i, k}}{a_{k, k}}, \quad i=k+1, \ldots, n</tex>.Если коэффициенты <tex> >1 </tex>, то ошибки, полученные на предыдущих шагах накапливаются. Чтобы этого избежать, применяется модификация метода Гаусса с выбором главного элемента. На каждом шаге к обычной схеме добавляется выбор максимального элемента по столбцу следующим образом:
 +
 +
Пусть по ходу исключения неизвестных получена система уравнений:
 +
 +
:<tex>x_i+\sum_{j=i+1}^{n}a_{ij}^i x^j = b_i^i, \quad \quad i=1, \ldots, k-1</tex>,
 +
:<tex>\sum_{j=k}^{n}a_{ij}^{k-1}x_j = b_i^{k-1}, \quad \quad i=k, \ldots, n</tex>.
 +
 +
Найдем такое <tex>l</tex>, что <tex>|a_{l,k}^{k-1}|=\max_{j=k,...,n} |a_{j,k}^{k-1}|</tex> и поменяем местами <tex>k</tex>-е и <tex>l</tex>-е уровнения.
 +
 +
Такое преобразование во многих случаях существенно уменьшает чувствительность решения к погрешностям округления при вычислениях.
 +
=== Итеративное улучшение результата ===
=== Итеративное улучшение результата ===
 +
 +
Если есть подозрение, что полученное решение <tex>x^{1}</tex> сильно искажено, то можно улучшить результат следующим образом. Величина <tex>b^1=b-Ax^{0}</tex> называется невязкой. Погрешность <tex>r^1=x-x^{1}</tex> удовлетворяет системе уравнений
 +
:<tex>Ar^1=Ax-Ax^1=b^1</tex>.
 +
Решая эту систему, получаем приближение <tex>r^{(1)}</tex> к <tex>r^1</tex> и полагаем
 +
:<tex>x^2=x^1+r^{(1)}</tex>.
 +
Если точность данного приближения неудовлетворительна, то повторяем эту операцию.
 +
 +
Процесс можно продолжать до тех пор, пока все компоненты <tex>r^{(i)}</tex> не станут достаточно малыми. При этом нельзя останавливать вычисления только потому, что все компоненты вектора невязки стали достаточно малыми: это может быть результатом плохой обусловленности матрицы коэффициентов.
 +
 +
==Числовой пример==
 +
 +
Рассмотрим для примера матрицу Вандермонда размером 7х7 и 2 различные правые части:
 +
<center><tex>\left( \begin{array}{ccc} 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ \\ 64 & 32 & 16 & 8 & 4 & 2 & 1\\ \\ 729 & 243 & 81 & 27 & 9 & 3 & 1 \\ \\ 4096 & 1024 & 256 & 64 & 16 & 4 & 1 \\ \\ 15625 & 3125 & 625 & 125 & 25 & 5 & 1 \\ \\ 46656 & 7776 & 1296 & 216 & 36 & 6 & 1 \\ \\ 117649 & 16807 & 2401 & 343 & 49 & 7 & 1\end{array}\right) \left( \begin{array}{c}x_1 \\ \\ x_2 \\ \\ x_3 \\ \\ x_4 \\ \\ x_5 \\ \\ x_6 \\ \\ x_7 \end{array}\right) = \left( \begin{array}{ccc}7 & 4 \\ \\ 127 & 85\\ \\ 1093 & 820\\ \\ 5461 & 4369 \\ \\ 19531 & 16276\\ \\ 55987 & 47989\\ \\ 137257 & 120100 \end{array}\right)</tex></center>
 +
 +
Данные системы были решены двумя способами. Тип данных - float. B итоге получили следующие результаты:
 +
 +
{| class="standard"
 +
|
 +
{| class="standard"
 +
! Способ решения
 +
|-
 +
! Номер столбца в <tex>b</tex>
 +
|-
 +
! Номер итерации
 +
|-
 +
! x1
 +
|-
 +
! x2
 +
|-
 +
! x3
 +
|-
 +
! x4
 +
|-
 +
! x5
 +
|-
 +
! x6
 +
|-
 +
! x7
 +
|-
 +
! Абс. погрешность
 +
|-
 +
! Отн. погрешность
 +
|-
 +
! Невязка
 +
|}
 +
|
 +
{| class="standard" width=40%
 +
! colspan=4| Обычный метод
 +
|-
 +
! colspan=2 | 1
 +
! colspan=2 | 2
 +
|-
 +
! 1 || 2 || 1 || 2
 +
|-
 +
|0.999991 || 1 || 0.999996 || 1
 +
|-
 +
|1.00019 || 1 || 7.4774e-005 || 2,33e-008
 +
|-
 +
|0.998404 || 1 || 0.999375 || 1
 +
|-
 +
|1.00667 || 1 || 0.00263727 || 1,12e-006
 +
|-
 +
|0.985328 || 1 || 0.994149 || 1
 +
|-
 +
|1.01588 || 1 || 0.00637817 || 3,27e-006
 +
|-
 +
|0.993538 || 1 || 0.99739 || 1
 +
|-
 +
| 0,045479 || 2,9826e-006 || 0,01818 ||8,8362e-006
 +
|-
 +
| 0,006497 || 4,2608e-007 || 0,0045451 ||2,209e-006
 +
|-
 +
| 0,040152 || 4,344e-005 || 0,083938 ||2,8654e-006
 +
|}
 +
|
 +
{| class="standard" width=40%
 +
! colspan=4| С выбором ведущего элемента по строке
 +
|-
 +
! colspan=2 | 1
 +
! colspan=2 | 2
 +
|-
 +
! 1 || 2 || 1 || 2
 +
|-
 +
| 1 || 1 || 1 || 1
 +
|-
 +
|1 || 1|| -3.57628e-005 || 1,836e-007
 +
|-
 +
|1.00001 || 1|| 1.00031 || 1
 +
|-
 +
|0.999942 || 1|| -0.00133276 || 7,16e-006
 +
|-
 +
|1.00005 || 1|| 1.00302 || 0,99998
 +
|-
 +
|1.00009 || 1|| -0.0033505 || 1,8e-005
 +
|-
 +
|0.99991 || 1|| 1.00139 || 0,99999
 +
|-
 +
|0,000298 || 4,3835e-007 || 0,009439 || 5,0683e-005
 +
|-
 +
| 4,2571e-005 || 6,2622e-008 || 0,0023542 || 1,2671e-005
 +
|-
 +
| 0,010622 || 9,8016e-007 || 0,29402 ||1,4768e-006
 +
|}
 +
|
 +
{| class="standard"
 +
! colspan=2| Реальные результаты
 +
|-
 +
! 1
 +
! 2
 +
|-
 +
! - || -
 +
|-
 +
|1 || 1
 +
|-
 +
|1 || 0
 +
|-
 +
|1 || 1
 +
|-
 +
|1 || 0
 +
|-
 +
|1 || 1
 +
|-
 +
|1 || 0
 +
|-
 +
|1 || 1
 +
|-
 +
| 0 || 0
 +
|-
 +
| 0 || 0
 +
|-
 +
| 0 || 0
 +
|}
 +
|}
 +
== Программа, реализующая метод на C++ ==
== Программа, реализующая метод на C++ ==
-
== Рекомендации программисту ==
+
 
 +
[[Медиа:Gauss_Elimination.zip| Gauss_Elimination.zip [30КБ] ]] - В архиве содержится исходный код, exe-файл и пример файла с входными данными
 +
 
 +
==Рекомендации пользователю==
 +
*Метод Гаусса удобно применять для систем маленькой и средней размерности (до порядка <tex>10^4</tex>). Для больших же размерностей или разреженных матриц более эффективными представляются итерационные методы.
 +
 
 +
*Рекомендуется использовать '''метод Гаусса с выбором главного элемента''' по столбцу, как более устойчивый к ошибкам, но при этом не требующий больших дополнительных затрат. В этом плане метод выбора главного элемента по строке и столбцу представляется менее эффективным, так как требует гораздо больше вычислительных затрат, но дает небольшую прибавку в точности.
 +
 
 +
*Итерационное улучшение результата стоит провести 2-3 раза, если погрешность уменьшается очень медленно - возможно, матрица плохо обусловлена, тогда метод в любом случае даст не лучшие результаты - лучше попробовать применить итерационные методы.
 +
 
 +
*Важно, чтобы данные считывались из файла правильно.
 +
При записи данных помните, что, '''целая и дробная части числа отделяются точкой''', а не запятой. В последнем случае данные будут считаны неправильно, но не будет выведено никакого сообщения об ошибке.
 +
 
 +
*Есть возможность менять точность вычислений (float и double), но в исходниках.
 +
 
 +
==Сноски==
 +
<references/>
 +
 
== Список литературы ==
== Список литературы ==
* ''Н.С.Бахвалов, Н.П.Жидков, Г.М.Кобельков'' Численные методы
* ''Н.С.Бахвалов, Н.П.Жидков, Г.М.Кобельков'' Численные методы
 +
* ''Д. Мак-Кракен, У.Дорн'' Численные методы и программирование на фортране
 +
* ''J.C.Nash'' Compact numerical methods for computers. Linear algebra anf function minimization. Second edition.
 +
* ''Nicholas J.Higham'' How Accurate is Gaussian Elimination?
 +
== Внешние ссылки ==
== Внешние ссылки ==
 +
*[http://unn.ru/rus/books/met_files/metgauss.doc Лабораторная работа по методу Гаусса в НГУ им. Лобачевского]
 +
*[http://ru.wikipedia.org/wiki/Метод_Гаусса Метод Гаусса на Википедии]
*[http://ru.wikipedia.org/wiki/%D0%9E%D0%B1%D1%80%D0%B0%D1%82%D0%BD%D0%B0%D1%8F_%D0%BC%D0%B0%D1%82%D1%80%D0%B8%D1%86%D0%B0 Обратная матрица]
*[http://ru.wikipedia.org/wiki/%D0%9E%D0%B1%D1%80%D0%B0%D1%82%D0%BD%D0%B0%D1%8F_%D0%BC%D0%B0%D1%82%D1%80%D0%B8%D1%86%D0%B0 Обратная матрица]
Строка 55: Строка 291:
{{stub}}
{{stub}}
-
[[Категория:Решение СЛАУ]]
+
[[Категория:Линейная алгебра]]
[[Категория:Учебные задачи]]
[[Категория:Учебные задачи]]

Текущая версия

Содержание

Постановка задачи

Дана система линейных алгебраических уравнений (СЛАУ), состоящая из n уравнений с n неизвестными :

(1)

\left\{\begin{array}{lcr} a_{11}x_1+\ldots+a_{1n}x_n &=& b_1 \\ \\ \cdots & & \\ \\ a_{n1}x_1+\ldots+a_{nn}x_n &=& b_n \\ \end{array} \right. \quad \Longleftrightarrow \quad A\vec{x}=\vec{b}, \quad A=\left( \begin{array}{ccc} a_{11} & \ldots & a_{1n}\\ \\ \cdots & & \\ \\ a_{n1} & \ldots & a_{nn} \end{array}\right),\quad \vec{b}=\left( \begin{array}{c}b_1 \\ \\ \vdots \\ \\ b_n \end{array} \right).

Предполагается, что существует единственное решение системы, то есть detA \neq 0.

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

Описание метода

Процесс решения системы линейных уравнений

(2)
\sum_{j=1}^n a_{ij}x_{j} = b_{i}, \qquad i=1,\ldots,n,

по методу Гаусса состоит из 2х этапов:

  • Прямой ход
    Система (2) приводится к треугольному виду
1. Предполагаем, что a_{11} \neq 0 . Тогда первое уравнение системы (2) делим на коэффициент a_{11}, в результате получаем уравнение
x_{1}+\sum_{j=2}^n a_{ij}^{1}x_{j} = b_{i}^{1}
.
Затем из каждого из оставшихся уравнений вычитается первое, умноженное на соответствующий коэффициент a_{i1}. В результате система преобразуются к виду:
\left( \begin{array}{ccc} 1 & a_{12}^{1} & \ldots & a_{1n}^{1}\\ \\ 0 & a_{22}^{1} & \ldots & a_{2n}^{1}\\ \\ \vdots & \vdots & \ddots & \vdots\\ \\ 0 & a_{n2}^{1} & \ldots & a_{nn}^{1} \end{array}\right) \left( \begin{array}{c}x_1 \\ \\ x_2 \\ \\ \vdots \\ \\ \\ x_n \end{array}\right) = \left( \begin{array}{c}b_1^1 \\ \\ b_2^1 \\ \\ \vdots \\ \\ b_n^1 \end{array}\right)
2. В предположении, что a_{22}^1 \neq 0, делим второе уравнение на коэффициент a_{22}^1 и исключаем неизвестное x_2 из всех последующих уравнений и т.д.
3. Получаем систему уравнений с треугольной матрицей:
(3)
\left( \begin{array}{ccc} 1 & a_{12}^n & \ldots & a_{1n}^n\\ \\ 0 & 1 & \ldots & a_{2n}^n\\ \\ \vdots & \vdots & \ddots & \vdots\\ \\ 0 & 0 & \ldots & 1 \end{array}\right) \left( \begin{array}{c}x_1 \\ \\ x_2 \\ \\ \vdots \\ \\ \\ x_n \end{array}\right) = \left( \begin{array}{c}b_1^n \\ \\ b_2^n \\ \\ \vdots \\ \\ b_n^n \end{array}\right)
  • Обратный ход
    Непосредственное определение неизвестных
1. Из n-го уравнения системы (3) определяем x_{n}:\: \quad x_n=b_n^n
2. Из (n-1)-го - определяем x_{n-1}:\: \quad x_{n-1}=b_{n-1}^n-a_{(n-1)n}^n x_n и т.д.

Анализ метода

Данный метод относится к классу прямых методов решения системы уравнений, а это значит, что за конечное число шагов можно получить точное решение, при условии, что входные данные ( матрица A и правая часть уравнения - b) заданы точно и вычисление ведется без округлений. Для получения решения требуется \frac{m^3}{3}+ m^2 - \frac{m}{3} умножений и делений, то есть порядка O(\frac{m^3}{3}) операций.

Условия, при которых метод выдает точное решение, на практике не выполнимы - неизбежны как ошибки входных данных, так и ошибки округления. Тогда встает вопрос: насколько точное решение можно получить, используя метод Гаусса, насколько метод корректен? Определим устойчивость решения относительно входных параметров. Наряду с исходной системой (1) рассмотрим возмущенную систему:

(A + \delta A)(x + \delta x)=b + \delta b

Пусть введена некоторая норма || . ||. CondA=|| A || * || A^{-1}|| - называется числом обусловленности матрицы A.

Возможны 3 случая:

1) \delta A = 0 , \quad \delta b \neq 0 :

\frac{||\delta x||}{||x||} &\leq& Cond A \frac{||\delta b ||}{|| b ||}

2) \delta A \neq 0 , \quad \delta b = 0 :

<tex>\frac{||\delta x||}{||x||} \leq Cond A \frac{||\delta A ||}{|| A ||}

3) \delta A \neq 0 , \quad \delta b \neq 0 :

<tex>\frac{||\delta x||}{||x||} \leq \frac{Cond A}{1-Cond A\frac{||\delta A ||}{|| A ||}}(\frac{||\delta A ||}{|| A ||}+ \frac{||\delta b ||}{|| b ||})

Число обусловленности матрицы A всегда \geq 1. Если оно велико ( Cond A \approx 10^k k \geq 2 ) , то говорят, что матрица A плохо обусловлена. В этом случае малые возмущения правых частей системы (1), вызванные либо неточностью задания исходных данных, либо вызванные погрешностями вычисления, существенно влияют на решение системы. Грубо говоря, если погрешность правых частей 10^{-l} , то погрешность решения будет 10^{-l+k} .

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

\left\{\begin{array}{lcr}5x_1+7x_2 &=& 12 \\ \\ 7x_1+10x_2 &=& 17 \\ \end{array} \right.

Она имеет решение (1, \quad 1).

Теперь рассмотрим возмущенную систему:

\left\{\begin{array}{lcr}5x_1+7x_2 &=& 12.075 \\ \\ 7x_1+10x_2 &=& 16.905 \\ \end{array} \right.

Решением такой системы будет вектор (2.415, \quad 0).

При совсем малом возмущении правой части получили несоизмеримо большое возмущение решения. Объяснить такую "ненадежность" решения можно тем, что матрица A почти вырожденная: прямые, соответствующие двум уравнениям, почти совпадают, что видно на графике:

Геометрическое представление системы двух линейных алгебраических уравнений, которая является почти вырожденной. Прямые, соответствующие двум уравнениям, почти совпадают.


Такой результат можно было предвидеть в силу плохой обусловленностью матрицы A = \left( \begin{array}{ccc} 5 & 7\\ \\ 7 & 10 \end{array}\right) : Cond A = 17^2 [1]

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

Способы оценки ошибок

1) Контрольная сумма: обычно применяется для предупреждения случайных погрешностей в процессе вычисления без помощи компьютеров.

Составляем контрольный столбец a_{n+1} = \left( a_{1,n+1}, \ldots, a_{n, n+1}\right), состоящий из контрольных элементов системы:

a_{i, n+1} = \sum_{j=1}^n a_{i, j} + b_i

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

2) Относительная погрешность известного решения позволяет без существенных дополнительных затрат получить суждение о погрешности решения.

Задается некоторый ветор z с компонентами, имеющими по возможности тот же порядок и знак, что и компоненты искомого решения [1]. Вычисляется вектор c=Az, и на ряду с исходной системой уравнения решается система Az=c.

Пусть x' и z' - реально получаемые решения этих систем. Суждение о погрешности x-x' искомого решения можно получить, основываясь на гипотезе: относительные погрешности при решении методом исключения систем с одной и той же матрицей и различными правыми частями, которыми являются соответственно величины \frac{||x-x'||}{||x'||} и \frac{||z-z'||}{||z'||}, отличаются не в очень большое число раз.

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

Наряду с исходной системой тем же методом решается система

\left( \alpha A\right)x' = \beta b, где \alpha и \beta - числа

Если бы не было погрешности округления, то выполнялось бы равенство для решений исходной и масштабированной систем: x=\alpha \beta^{-1}x'. Поэтому при \alpha и \beta, не являющихся степенями двойки, сравнение векторов x и \alpha \beta^{-1}x' дает представление о величине вычислительной погрешности [1]

Улучшение метода исключения Гаусса

Рассмотренные ниже модификации метода Гаусса позволяют уменьшить погрешность результата.

Выбор главного элемента

Основное увеличение ошибки в методе происходит во время прямого хода, когда ведущая k-я строка умножается на коэффициенты \frac{a_{i, k}}{a_{k, k}}, \quad i=k+1, \ldots, n.Если коэффициенты >1 , то ошибки, полученные на предыдущих шагах накапливаются. Чтобы этого избежать, применяется модификация метода Гаусса с выбором главного элемента. На каждом шаге к обычной схеме добавляется выбор максимального элемента по столбцу следующим образом:

Пусть по ходу исключения неизвестных получена система уравнений:

x_i+\sum_{j=i+1}^{n}a_{ij}^i x^j = b_i^i, \quad \quad i=1, \ldots, k-1,
\sum_{j=k}^{n}a_{ij}^{k-1}x_j = b_i^{k-1}, \quad \quad i=k, \ldots, n.

Найдем такое l, что |a_{l,k}^{k-1}|=\max_{j=k,...,n} |a_{j,k}^{k-1}| и поменяем местами k-е и l-е уровнения.

Такое преобразование во многих случаях существенно уменьшает чувствительность решения к погрешностям округления при вычислениях.

Итеративное улучшение результата

Если есть подозрение, что полученное решение x^{1} сильно искажено, то можно улучшить результат следующим образом. Величина b^1=b-Ax^{0} называется невязкой. Погрешность r^1=x-x^{1} удовлетворяет системе уравнений

Ar^1=Ax-Ax^1=b^1.

Решая эту систему, получаем приближение r^{(1)} к r^1 и полагаем

x^2=x^1+r^{(1)}.

Если точность данного приближения неудовлетворительна, то повторяем эту операцию.

Процесс можно продолжать до тех пор, пока все компоненты r^{(i)} не станут достаточно малыми. При этом нельзя останавливать вычисления только потому, что все компоненты вектора невязки стали достаточно малыми: это может быть результатом плохой обусловленности матрицы коэффициентов.

Числовой пример

Рассмотрим для примера матрицу Вандермонда размером 7х7 и 2 различные правые части:

\left( \begin{array}{ccc} 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ \\ 64 & 32 & 16 & 8 & 4 & 2 & 1\\ \\ 729 & 243 & 81 & 27 & 9 & 3 & 1 \\ \\ 4096 & 1024 & 256 & 64 & 16 & 4 & 1 \\ \\ 15625 & 3125 & 625 & 125 & 25 & 5 & 1 \\ \\ 46656 & 7776 & 1296 & 216 & 36 & 6 & 1 \\ \\ 117649 & 16807 & 2401 & 343 & 49 & 7 & 1\end{array}\right) \left( \begin{array}{c}x_1 \\ \\ x_2 \\ \\ x_3 \\ \\ x_4 \\ \\ x_5 \\ \\ x_6 \\ \\ x_7 \end{array}\right) = \left( \begin{array}{ccc}7 & 4 \\ \\ 127 & 85\\ \\ 1093 & 820\\ \\ 5461 & 4369 \\ \\ 19531 & 16276\\ \\ 55987 & 47989\\ \\ 137257 & 120100 \end{array}\right)

Данные системы были решены двумя способами. Тип данных - float. B итоге получили следующие результаты:

Способ решения
Номер столбца в b
Номер итерации
x1
x2
x3
x4
x5
x6
x7
Абс. погрешность
Отн. погрешность
Невязка
Обычный метод
1 2
1 2 1 2
0.999991 1 0.999996 1
1.00019 1 7.4774e-005 2,33e-008
0.998404 1 0.999375 1
1.00667 1 0.00263727 1,12e-006
0.985328 1 0.994149 1
1.01588 1 0.00637817 3,27e-006
0.993538 1 0.99739 1
0,045479 2,9826e-006 0,01818 8,8362e-006
0,006497 4,2608e-007 0,0045451 2,209e-006
0,040152 4,344e-005 0,083938 2,8654e-006
С выбором ведущего элемента по строке
1 2
1 2 1 2
1 1 1 1
1 1 -3.57628e-005 1,836e-007
1.00001 1 1.00031 1
0.999942 1 -0.00133276 7,16e-006
1.00005 1 1.00302 0,99998
1.00009 1 -0.0033505 1,8e-005
0.99991 1 1.00139 0,99999
0,000298 4,3835e-007 0,009439 5,0683e-005
4,2571e-005 6,2622e-008 0,0023542 1,2671e-005
0,010622 9,8016e-007 0,29402 1,4768e-006
Реальные результаты
1 2
- -
1 1
1 0
1 1
1 0
1 1
1 0
1 1
0 0
0 0
0 0

Программа, реализующая метод на C++

Gauss_Elimination.zip [30КБ] - В архиве содержится исходный код, exe-файл и пример файла с входными данными

Рекомендации пользователю

  • Метод Гаусса удобно применять для систем маленькой и средней размерности (до порядка 10^4). Для больших же размерностей или разреженных матриц более эффективными представляются итерационные методы.
  • Рекомендуется использовать метод Гаусса с выбором главного элемента по столбцу, как более устойчивый к ошибкам, но при этом не требующий больших дополнительных затрат. В этом плане метод выбора главного элемента по строке и столбцу представляется менее эффективным, так как требует гораздо больше вычислительных затрат, но дает небольшую прибавку в точности.
  • Итерационное улучшение результата стоит провести 2-3 раза, если погрешность уменьшается очень медленно - возможно, матрица плохо обусловлена, тогда метод в любом случае даст не лучшие результаты - лучше попробовать применить итерационные методы.
  • Важно, чтобы данные считывались из файла правильно.

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

  • Есть возможность менять точность вычислений (float и double), но в исходниках.

Сноски


Список литературы

  • Н.С.Бахвалов, Н.П.Жидков, Г.М.Кобельков Численные методы
  • Д. Мак-Кракен, У.Дорн Численные методы и программирование на фортране
  • J.C.Nash Compact numerical methods for computers. Linear algebra anf function minimization. Second edition.
  • Nicholas J.Higham How Accurate is Gaussian Elimination?

Внешние ссылки

См. также

Источник — «http://machinelearning.ru/wiki/index.php?title=%D0%9C%D0%B5%D1%82%D0%BE%D0%B4%D1%8B_%D0%B8%D1%81%D0%BA%D0%BB%D1%8E%D1%87%D0%B5%D0%BD%D0%B8%D1%8F_%D0%93%D0%B0%D1%83%D1%81%D1%81%D0%B0»
Личные инструменты