Лекция № 4: Приведение матрицы к треугольному виду

Затвитить
эту лекцию
Опубликовано 12‑го октября 2008 г.
Обновлено 31‑го октября 2011 г.

Приведение матрицы к треугольному виду требуется, прежде всего, для решения систем линейных уравнений (СЛАУ). Однако, для простоты весь материал лекции будет рассмотрен на примере вычисления определителя матрицы. Ведь самый быстрый способ вычислить определитель — привести матрицу к треугольному виду, а затем перемножить её диагональные элементы. Если же вы хотите применить материал этой лекции к решению СЛАУ, то единственное что вам нужно сделать — добавить к операциям со строками матрицы такие же операции с элементами правой части СЛАУ (ну и не забыть найти значения неизвестных величин с помощью обратного хода алгоритма Гаусса).

Везде в этой лекции стро́ки и столбцы матриц нумеруются с нуля, чтобы более соответствовать языку программирования C++. Элементы матрицы обозначаются, как $a_{ij}$ , где i=0,...,N-1 — номер строки́, j=0,...,N-1 — номер столбца. Для сокращения записи стро́ки матрицы обозначаются, как $\vec a_i$ : с ними можно производить операции так же, как с векторами.

Алгоритм Гаусса

Матрицу можно привести к треугольному виду при помощи алгоритма Гаусса. Идея алгоритма заключается в том, что элементы $a_{ij}$ , i>j обнуляются путём вычитания строки́ номер j, умноженной на $a_{ij}/a_{jj}$ , из строки́ номер i. Формально алгоритм выглядит так:

$\begin{array}{l} \text{For }j=0,...,N-2\\ \left|\quad\begin{array}{l} \text{For }i=j+1,...,N-1\\ \left|\quad \vec a_i \leftarrow \vec a_i-\cfrac{a_{ij}}{a_{jj}}\cdot\vec a_j \right.\\ \end{array}\right.\\ \end{array}$

(1)

Алгоритм строит верхнетреугольную матрицу, изменяя исходную матрицу. Если исходная матрица должна сохраниться, то перед началом работы алгоритма нужно создать её копию.

Несмотря на то, что в алгоритме (1) показаны два вложенных цикла, на самом деле в нём три цикла, так как вычитание строк — операция над несколькими элементами. Соответственно, время работы алгоритма составляет $O\left(N^3\right)$ .

Алгоритм Гаусса не будет работать, если в какой-то момент окажется $a_{jj}=0$ в формуле (1). Необходимое и достаточное условие того, что в алгоритме не возникнет деления на ноль — все главные миноры матрицы должны быть отличны от нуля.

После приведения матрицы к верхнетреугольному виду её определитель можно вычислить перемножением диагональных элементов:

$\begin{array}{l} d\leftarrow a_{00}\\ \text{For }i=1,...,N-1\\ \left|\quad d\leftarrow d\cdot a_{ii} \right.\\ \end{array}$

(2)

Определитель треугольной матрицы, вычисленный по формуле (2), будет равен определителю исходной матрицы, так как определитель матрицы не меняется, если к одной строке матрицы прибавить другую строку, умноженную на произвольное число.

Пример работы алгоритма Гаусса

Вычислим определитель следующей матрицы по приведённому выше алгоритму:

$\begin{pmatrix}3&2&3&4\\ 4&4&3&2\\ 1&4&4&3\\ 2&3&1&1\end{pmatrix}$

(3)

Предположим, что компьютер при вычислениях хранит не более 5-ти значащих цифр: так более наглядно видны накапливающиеся ошибки округления. Умножим нулевую строку на 4/3 и вычтем её из первой строки́:

$\begin{pmatrix}3&2&3&4\\ 0&1.3333&-1&-3.3333\\ 1&4&4&3\\ 2&3&1&1\end{pmatrix}$

(4)

Аналогично вычтем нулевую строку из второй и третьей строк (умножив её перед этим на 1/3 и 2/3 соответственно):

$\begin{pmatrix}3&2&3&4\\ 0&1.3333&-1&-3.3333\\ 0&3.3333&3&1.6667\\ 0&1.6667&-1&-1.6667\end{pmatrix}$

(5)

Вычтем первую строку из второй и третьей (стро́ки у нас нумеруются с нуля!), умножив её перед этим на 3.3333/1.3333 и 1.6667/1.3333 соответственно:

$\begin{pmatrix}3&2&3&4\\ 0&1.3333&-1&-3.3333\\ 0&0&5.5&10.0001\\ 0&0&0.2501&2.5001\end{pmatrix}$

(6)

Наконец, вычтем из третьей строки́ вторую, умноженную на 0.2501/5.5 :

$\begin{pmatrix}3&2&3&4\\ 0&1.3333&-1&-3.3333\\ 0&0&5.5&10.0001\\ 0&0&0&2.0454\end{pmatrix}$

(7)

Перемножая диагональные элементы матрицы, получим значение определителя: 44.9977 . Нетрудно догадаться, что если бы мы вычисляли определитель точно (например, с использованием дробей), то он был бы равен .

Перестановка строк

Мы видим, что по мере работы алгоритма накапливается погрешность. Причём эта погрешность будет тем большей, чем ближе к нулю знаменатель в выражении (1). Более того, в какой-то момент этот знаменатель может оказаться нулевым, что приведёт к остановке алгоритма.

Замечу, что нулевой знаменатель в выражении (1) ещё не означает равенства нулю определителя всей матрицы; также этот факт не означает невозможности решения СЛАУ с данной матрицей. Кроме того, возможна такая ситуация, когда знаменатель в выражении (1) из-за ошибок округления не равен нулю в то время, когда он должен быть равен нулю. В этом случае алгоритм даст результат, сколь угодно сильно отличающийся от точного.

Для борьбы с описанными выше проблемами используется перестановка строк: перед каждым вычитанием строки́ из всех последующих (перед вторым циклом алгоритма (1)) находится строка номер k (среди строк i=j,...,N-1 ) с максимальным по модулю элементом $a_{jk}$ , и, если $k\ne j$ , производится перестановка строк k и j. После перестановки строка с максимальным по модулю первым ненулевым элементом оказывается на месте той строки́, которая вычитается из всех следующих. И в знаменатель выражения (1) попадает максимальный элемент среди всех, которые находятся снизу от $a_{jj}$ .

$\begin{array}{l} \text{For }j=0,...,N-2\\ \left|\quad\begin{array}{l} \text{Find line }k\ \left(k=j,...,N-1\right)\text{ having maximum }\left|a_{jk}\right|\\ \text{If }k\ne j\text{ then swap lines }k\text { and }j\\ \text{For }i=j+1,...,N-1\\ \left|\quad\vec a_i \leftarrow \vec a_i-\cfrac{a_{ij}}{a_{jj}}\cdot\vec a_j\right.\\ \text{End} \end{array}\right.\\ \text{End} \end{array}$

(8)

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

Если применить алгоритм (8) к матрице (3), то получим следующую матрицу (с учётом округления всех промежуточных результатов до пяти значащих цифр):

$\begin{pmatrix}4&4&3&2\\ 0&3&3.25&2.5\\ 0&0&1.8333&3.3333\\ 0&0&0&2.0454\end{pmatrix}$

(9)

За время работы алгоритма было произведено 2 перестановки, поэтому определитель равен 44.9980 , что немного точнее, чем результат, полученный ранее.

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

Программная реализация алгоритма Гаусса

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

Рисунок 1. Структура данных для хранения матрицы

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

Пусть мы программируем на C++, и числа матрицы имеют некоторый тип T, тогда указатели на стро́ки будут иметь тип T*, а сама матрица M будет иметь тип T**. Учитывая это, алгоритм (8) может быть реализован следующим образом:

template<class T> T GaussDeterminant(T **const M, int const N) 

{

    T determinant(1);

    int exchanges(0);

    for(int l1(0); l1<N-1; ++l1)

    { //Перебираю все строки матрицы, кроме последней

        int maxN( l1 );

        T maxValue( fabs( M[l1][l1] ) );

        for(int l2(l1+1); l2<N; ++l2)

        { //Нахожу строку с максимальным по модулю элементом

            T const value( fabs( M[l2][l1] ) );

            if( value > maxValue ) { maxN=l2; maxValue=value; }

        }

        if( maxN > l1 )

        { //Нужен обмен

            T *const temp(M[l1]); M[l1]=M[maxN]; M[maxN]=temp;

            ++exchanges;

        } else { //Обмен не нужен, но нужна проверка на ноль

            if( maxValue == T(0) ) return maxValue;

        }

        T const value( M[l1][l1] );

        determinant *= value;

        for(int l2(l1+1); l2<N; ++l2)

        { //Вычитаю строку из всех последующих

            T const k( M[l2][l1]/value );

            M[l2][l1] = T(0);

            for(int c(l1+1); c<N; ++c) M[l2][c] -= M[l1][c]*k;

        }

    }

    determinant *= M[N-1][N-1];

    if(exchanges%2) return -determinant; else return determinant;

}

Листинг 1. Реализация алгоритма Гаусса с перестановкой строк

Модификация Барейса

Если матрица состоит из целых чисел, алгоритм Барейса позволяет привести её к треугольному виду с использованием только целочисленной арифметики!

Мы видели, что самый большой вклад в накапливаемую погрешность даёт деление на $a_{ii}$ в алгоритмах (1) и (8). Попытаемся избавиться от этого деления. Прежде всего заметим, что если сразу после вычитания одной строки́ из другой в формуле (1) умножить результат на $a_{ii}$ , то получится выражение, не содержащее деления:

$\vec a_i \leftarrow a_{jj}\cdot\vec a_i-a_{ij}\cdot\vec a_j.$

(10)

Уже с такой модификацией алгоритм можно использовать для решения СЛАУ, однако данный метод имеет следующие проблемы:

Модификация эквивалентна многократному умножению строк матрицы на её элементы. Если элементы матрицы больше единицы, то для матриц большого размера чи́сла к концу работы алгоритма вырастут очень значительно. Возможно переполнение используемого типа данных. Аналогично, если элементы меньше единицы, то к концу работы алгоритма числа станут очень малы.
Полученная треугольная матрица не пригодна для простого вычисления определителя исходной матрицы. Ибо, как несложно посчитать, определитель треугольной матрицы теперь равен определителю исходной матрицы, умноженному на $\left(a_{00}\right)^{N-1}\cdot\left(a_{11}\right)^{N-2}\cdot...\cdot\left(a_{N-2,N-2}\right)^1$ — диагональные элементы результирующей матрицы.

Барейс (Bareiss) доказал, что правая часть выражения (10) может быть поделена нацело на $a_{j-1,j-1}$ (а не на $a_{jj}$ , как в алгоритме Гаусса). При этом предполагается, что $a_{-1,-1}\equiv 1$ :

$\vec a_i \leftarrow \frac{a_{jj}\cdot\vec a_i-a_{ij}\cdot\vec a_j}{a_{j-1,j-1}}.$

(11)

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

Ещё одной интересной особенностью алгоритма является то, что в нём всего одно «нескомпенсированное» умножение для каждого диагонального элемента, то есть определитель D_T результирующей треугольной матрицы равен определителю D_R исходной матрицы, умноженному на $a_{00}\cdot a_{11}\cdot ...\cdot a_{N-2,N-2}$ (здесь имеются в виду диагональные элементы результирующей матрицы):

$D_T=D_R\times a_{00}\cdot a_{11}\cdot ...\cdot a_{N-2,N-2},$

(12)

но

$D_T=a_{00}\cdot a_{11}\cdot ...\cdot a_{N-1,N-1},$

(13)

поэтому

$D_R=\frac{a_{00}\cdot a_{11}\cdot ...\cdot a_{N-1,N-1}}{a_{00}\cdot a_{11}\cdot ...\cdot a_{N-2,N-2}}=a_{N-1,N-1}.$

(14)

То есть определитель исходной матрицы равен правому нижнему элементу матрицы, полученной в результате алгоритма Барейса!

Формальная запись алгоритма Барейса:

$\begin{array}{l} a_{-1,-1}\equiv 1\\ \text{For }j=0,...,N-2\\ \left|\quad\begin{array}{l} \text{For }i=j+1,...,N-1\\ \left|\quad\vec a_i \leftarrow \cfrac{a_{jj}\cdot\vec a_i-a_{ij}\cdot\vec a_j}{a_{j-1,j-1}}\right. \end{array}\right.\\ \end{array}$

(15)

Алгоритм Барейса так же, как и алгоритм Гаусса, может быть улучшен путём добавления перестановки строк.

Пример работы алгоритма Барейса

Рассмотрим работу алгоритма на примере матрицы (3). Умножаем первую, вторую, и третью стро́ки на 3, и отнимаем от них нулевую строку, умноженную на 4, 1, и 2 соответственно. Деление на этом шаге не нужно, так как $a_{-1,-1}\equiv 1$ . Получим:

$\begin{pmatrix}3&2&3&4\\ 0&4&-3&-10\\ 0&10&9&5\\ 0&5&-3&-5\end{pmatrix}$

(16)

Теперь вычитаем из второй и третьей строк, умноженных на 4, первую строку, умноженную на 10 и 5 соответственно (стро́ки нумеруются с нуля). После этого делим стро́ки на 3:

$\begin{pmatrix}3&2&3&4\\ 0&4&-3&-10\\ 0&0&22&40\\ 0&0&1&10\end{pmatrix}$

(17)

Наконец, умножаем последнюю строку на 22, вычитаем из неё предпоследнюю, и делим результат на 4:

$\begin{pmatrix}3&2&3&4\\ 0&4&-3&-10\\ 0&0&22&40\\ 0&0&0&45\end{pmatrix}$

(18)

Обратите внимание на нижний правый элемент: это определитель исходной матрицы. И нам даже не потребовались дроби или числа с плавающей точкой для приведения матрицы к треугольному виду и вычисления этого определителя.

Программная реализация алгоритма Барейса

Реализация алгоритма Барейса на языке программирования C++ почти не отличается от соответствующей реализации алгоритма Гаусса. Отличия показаны восклицательными знаками:

template<class T> T BareissDeterminant(T **const M, int const N)

{

    T denom(1); //!

    int exchanges(0);

    for(int l1(0); l1<N-1; ++l1)

    { //Перебираю все строки матрицы, кроме последней

        int maxN( l1 );

        T maxValue( fabs( M[l1][l1] ) );

        for(int l2(l1+1); l2<N; ++l2)

        { //Нахожу строку с максимальным по модулю элементом

            T const value( fabs( M[l2][l1] ) );

            if( value > maxValue ) { maxN=l2; maxValue=value; }

        }

        if( maxN > l1 )

        { //Нужен обмен

            T *const temp(M[l1]); M[l1]=M[maxN]; M[maxN]=temp;

            ++exchanges;

        } else { //Обмен не нужен, но нужна проверка на ноль

            if( maxValue == T(0) ) return maxValue;

        }

        T const value1( M[l1][l1] ); //!

          //!

        for(int l2(l1+1); l2<N; ++l2)

        { //Вычитаю строку из всех последующих

            T const value2( M[l2][l1] ); //!

            M[l2][l1] = T(0);

            for(int c(l1+1); c<N; ++c) //!

                M[l2][c]=(M[l2][c]*value1-M[l1][c]*value2)/denom;

        }

        denom = value1; //!

    }

    //!

    if(exchanges%2) return -M[N-1][N-1]; else return M[N-1][N-1];

}

Листинг 2. Реализация алгоритма Барейса с перестановкой строк

Проверка алгоритмов на матрице Гильберта

Матрицей Гильберта называется матрица вида:

$H_{ij}=\frac{1}{i+j+1},\quad i,j=0,...,N-1;$

(19)

или

$H^N=\begin{pmatrix}1&1/2&\ldots &1/N\\ 1/2&1/3&\ldots &1/(N+1)\\ \vdots &\vdots &\ddots &\vdots\\ 1/N&1/(N+1)&\ldots &1/(2N-1)\end{pmatrix}.$

(20)

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

Определители обратной матрицы Гильберта образуют следующую последовательность чисел (определитель обратной матрицы равен обратной величине определителя):

$\begin{aligned} \left|H^1\right|^{-1}&=1,\\ \left|H^2\right|^{-1}&=12,\\ \left|H^3\right|^{-1}&=2160,\\ \left|H^4\right|^{-1}&=6048000,\\ \left|H^5\right|^{-1}&=266716800000,\\ \left|H^6\right|^{-1}&=186313420339200000,\\ \left|H^7\right|^{-1}&=2067909047925770649600000,\\ \left|H^8\right|^{-1}&=365356847125734485878112256000000,\\ \left|H^9\right|^{-1}&=1028781784378569697887052962909388800000000,\\ \left|H^{10}\right|^{-1}&={}\\ {}=46206&89394791469131629562883903627872698368000000000. \end{aligned}$

(21)

$\begin{aligned} \left|H^1\right|^{-1}&=1,\\ \left|H^2\right|^{-1}&=12,\\ \left|H^3\right|^{-1}&=2160,\\ \left|H^4\right|^{-1}&=6048000,\\ \left|H^5\right|^{-1}&=266716800000,\\ \left|H^6\right|^{-1}&=186313420339200000,\\ \left|H^7\right|^{-1}&=2067909047925770649600000,\\ \left|H^8\right|^{-1}&=365356847125734485878112256000000,\\ \left|H^9\right|^{-1}&=1028781784378569697887052962909388800000000,\\ \left|H^{10}\right|^{-1}&=4620689394791469131629562883903627872698368000000000. \end{aligned}$

(21)

Попробуем найти $\left|H^{10}\right|^{-1}$ с помощью описанных алгоритмов, используя числа двойной точности:

	Без перестановки строк	С перестановкой строк
Алгоритм Гаусса	$4.620287...\cdot 10^{52}$	$4.620207...\cdot 10^{52}$
Алгоритм Барейса	$4.620204...\cdot 10^{52}$	$4.620375...\cdot 10^{52}$

Таблица 1. Обратная величина определителя матрицы Гильберта,
вычисленного различными способами

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

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

Параллельная реализация

Рассмотрим пока возможность параллельной реализации описанного алгоритма с теоретической точки зрения (практической реализацией мы займёмся в одной из следующих лекций).

Если посмотреть листинг 2, то можно заметить, что «вычитание» строки́ из всех последующих — независимые операции, и обработка элементов каждой строки́ — тоже независимые операции. Но имеется зависимость обработки элементов одной строки́ от обработки элементов другой строки́ через переменную value2.

Можно либо давать каждому вычислителю на обработку целое число строк и свой экземпляр переменной value2 (что может привести к дисбалансу нагрузки, так как число строк может не поделиться нацело на число вычислителей), либо вообще не использовать переменную value2 (а «первые» элементы строк обнулять в конце, или вообще не обнулять, — на значение определителя это не повлияет), и разбить всё множество операций «вычитания» поровну между процессорами. В последнем случае дисбаланс нагрузки если и возникнет, то это будет не одна строка, как при разбиении по строкам, а один элемент, что гораздо менее заметно. Заведомо нулевые элементы векторов не обрабатываются, поэтому общее число операций на этом этапе равно $4\cdot \left(N-l_1-1\right)^2$ штук.

А вот нахождение строки́ с максимальным по модулю элементом — зависимые операции (они оперируют общими переменными maxN и maxValue). Можно, конечно, применить в этом случае алгоритм сдваивания, но гораздо проще оставить эту часть алгоритма последовательной, ибо нахождение максимального элемента — малая часть операций (примерно $4\cdot \left(N-l_1-1\right)$ штук) по сравнению с вычитанием строк.

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

Найдём параметр $\alpha$ для закона Амдаля. По определению, $\alpha$ равен отношению количества последовательных операций к общему числу операций:

$\alpha=\frac{N_1}{N_1+N_p}.$

(22)

Число последовательных операций:

$N_1 = \sum_{l_1=0}^{N-2} 4\cdot \left(N-l_1-1\right) = 4\cdot\frac{N^2-N}{2}.$

(23)

Число параллельных операций:

$N_p = \sum_{l_1=0}^{N-2} 4\cdot \left(N-l_1-1\right)^2 = 4\cdot\frac{2N^3-3N^2+N}{6}.$

(24)

Подставляя выражения (23) и (24) в (22), получим:

$\alpha\left(N\right)=\frac{3}{2N+2}.$

(25)

Полученное выражение можно подставить в выражение (7) предыдущей лекции, и использовать его, как оценку сверху для возможного ускорения алгоритма. Видно, что чем больше N, тем меньше будет $\alpha$ , и тем ближе кривая ускорения будет к прямой максимального ускорения S_p=p . Эта закономерность часто встречается на практике: чем больше исходных данных даётся алгоритму, тем обычно с большей эффективностью этот алгоритм можно распараллелить.

Что почитать по теме

Chee K. Yap, “Fundamental Problems in Algorithmic Algebra” — лекции по численным алгоритмам алгебры, с доказательствами всех утверждений. Алгоритм Барейса описан в лекции 10 (.ps-файлы более новые — читайте их).

Затвитить
эту лекцию

13 отзывов на «приведение матрицы к треугольному виду»

Отзыв № 1
Автор: Евгения. Дата: 10-го января 2009 г. Время: 00:29.

Отличный конспект. Единственное, еще бы метод Халецкого сюда добавить — было бы вообще СУПЕР!!!!!!!!!

Отзыв № 2
Автор: Александр. Дата: 16-го декабря 2009 г. Время: 18:05.

а я так ни чего и не понял!=(

Отзыв № 3
Автор: Михаил. Дата: 14-го февраля 2010 г. Время: 00:12.

Здравствуйте, если я решаю систему линейных однородных уравнений, то тогда мне надо, чтобы определитель как раз равнялся 0. Можно ли с помощью этого метода определить, равен ли определитель 0 (с некоторой степенью точности)?

Отзыв № 4 (ответ на отзыв № 3).
Автор: Антон. Дата: 14-го февраля 2010 г. Время: 00:34.

На самом деле это очень сложный вопрос. Что́ считать равенством или не равенством, если вычисления выполняются приближённо?

Если у вас все коэффициенты в исходной матрице — целые числа, то вам очень повезло — берите алгоритм Барейса с перестановкой строк (листинг 2) и, если функция вернула 0, то у вас определитель в точности равен нулю.

Если исходные числа рациональные — то выполняйте все вычисления в дробях, и полу́чите точный результат.

Если же ваши числа — произвольные числа с плавающей точкой, то ситуация гораздо хуже. Ведь вам для решения нужно привести матрицу не к треугольной, а к трапецевидной форме: то есть при обнаружении нулевого подстолбца нужно смещаться вправо на один столбец. Каждое такое «смещение вправо» уменьшает ранг матрицы на единицу. Ранг — это целое число, которое вам нужно знать точно, ибо ранг определяет размерность подпространства решений СЛАУ. Малейшие погрешности при вычислениях могут привести к какому угодно изменению ранга матрицы. Поэтому совершенно непонятно, что делать в этом случае. Что́ считать числом, близким к нулю?

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

Отзыв № 5
Автор: Стас. Дата: 21-го марта 2010 г. Время: 02:03.

А какой(или чей) метод нахождения определителя матрицы самый точный?

Отзыв № 6 (ответ на отзыв № 5).
Автор: Антон. Дата: 21-го марта 2010 г. Время: 10:50.

Всё зависит от того, в каком виде у вас представлены данные.

Если данные у вас представлены в целых числах, и вычисляете вы в целых числах с достаточным диапазоном, то алгоритм Барейса даст точный ответ (не важно, с перестановкой строк или без, лишь бы деления на ноль не было). Если данные представлены в дробях, и вычисления вы ведёте в дробях, то все методы дадут точный ответ. Если же у вас числа с плавающей точкой ограниченной точности (и исходные данные того же типа), то самый точный ответ будет у алгоритма Барейса с перестановкой строк.

Отзыв № 7
Автор: Jason. Дата: 9-го апреля 2010 г. Время: 15:13.

ОГРОМНЕШОЕ СПАСИБО, особенно за реализацию алгоритма Гаусса
3 дня промучался над ним а ощибка была всего в недостоющей 1-й строки проги

ещё раз спасибо

Отзыв № 8
Автор: ziwert. Дата: 11-го июля 2010 г. Время: 18:01.

Спасибо большое за лекцию, очень помогла при решении практики студенту-математику)

Отзыв № 9
Автор: Чайник!!!. Дата: 6-го октября 2010 г. Время: 16:35.

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

Отзыв № 10
Автор: Ержан. Дата: 13-го ноября 2010 г. Время: 17:49.

Здравствуйте! Спасибо за столь хороший конспект! Очень помогло. Но не до конца( Учусь программировать в СИ, надо написать программу, которая будет решать СЛАУ методом Гаусса. В принципе почти сделал, но никак не могу написать условие, когда матрица не имеет решения или когда бесконечное множество решений. Если кто может — помогите пожалуйста. Скидываю код в СИ:

#include
#include
#define N 20
#define M (N+1)

void glavelem( int k, double mas[][N + 1], int n, int otv[] )
{
int i, j, i_max = k, j_max = k;
double temp;
//Ищем максимальный по модулю элемент
for ( i = k; i < n; i++ )
for ( j = k; j < n; j++ )
if ( fabs( mas[i_max] [j_max] ) < fabs( mas[i] [j] ) )
{
i_max = i;
j_max = j;
}
//Переставляем строки
for ( j = k; j < n + 1; j++ )
{
temp = mas[k] [j];
mas[k] [j] = mas[i_max] [j];
mas[i_max] [j] = temp;
}
//Переставляем столбцы
for ( i = 0; i N )
{
printf("Слишком большое кол-во переменных.\n");
return 0;
}
}
while ( N < n );

printf( "Задайте %d элементов матрицы:\n", n*(n+1));
for ( i = 0; i < n; i++ )
for ( j = 0; j < n + 1; j++ )
scanf( "%lf", &mas[i][j] );

/*——————————————————*/
//Сначала все корни по порядку
for ( i = 0; i < n + 1; i++ )
otv[i] = i;
//Прямой ход метода Гаусса
for ( k = 0; k = k; j— )
mas[k] [j] /= mas[k] [k];
for ( i = k + 1; i = k; j— )
mas[i] [j] -= mas[k] [j] * mas[i] [k];
}
//Обратный ход
for ( i = 0; i = 0; i— )
for ( j = i + 1; j < n; j++ )
x[i] -= x[j] * mas[i] [j];
//Вывод результата
printf( "Одно решение:\n" );
for ( i = 0; i < n; i++ )
for ( j = 0; j < n; j++ )
if ( i == otv[j] )
{ //Расставляем корни по порядку
printf( "x%d = %f\n", y++, x[j] );
break;
}
return ( 0 );
}

П.С заранее спасибо тем, кто откликнется. Работу сдавать завтра. По возможности напишите недостающие строки кода:
1)если матрица не имеет решения
2)если матрица имеет бесконечное множество решений

Отзыв № 11
Автор: Ержан. Дата: 13-го ноября 2010 г. Время: 18:32.

не знаю почему, но часть кода проглатывает… но если вдруг кто то захочет помочь — пишите на ящик: apelsin_ast@mail.ru
Надеюсь на Вашу помощь

Отзыв № 12
Автор: Yaroslav. Дата: 22-го декабря 2010 г. Время: 03:29.

<< "(практической реализацией мы займёмся в одной из следующих лекций)."

Почему нет обещанной параллельной реализации приведения к треугольному виду?
Это наиболее интересует.

Отзыв № 13
Автор: Наталья. Дата: 24-го мая 2011 г. Время: 18:34.

Спасибо, отличный материал..очень помогли!

Image Processing

Лицензия

Реклама