47

Лекция № 23

11. Алгоритмы упорядочения разреженных матриц

Оглавление

11.1. Особенности решения матричных уравнений. 2

11.2. Проблема упорядочения 9

11.3. Матрицы и графы 17

11.4. Методы упорядочения матриц и сопутствующие алгоритмы. 22

11.4.1.Ленточные методы 22

11.4.2. Профильные методы. 29

11.4.3. Универсальные разреженные методы. 39

11.4.4. Определение начального узла 42

11.5. Сравнение методов упорядочения матриц. 46

Общие сведения о матрицах.

А + В = В + А

АВ≠ВА;

АВ=С;

Транспонированная матрица А^Т получается отражением матрицы А относительно её диагонали, т.е. .

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

Единичная матрица обязательно квадратная. IA = AI.

Обратная матрица А^-1 есть такая матрица, для которой выполняется равенство А^-1А=1. Не всякая матрица имеет обратную. Действительно, произведение АВ может быть равным нулю, даже если А или В не были нулевым матрицами

.

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

Следующие высказывания эквивалентны:

матрица А – невырожденная;

обратная матрица А^-1 существует;

столбцы матрицы А линейно независимы;

строки матрицы А линейно независимы;

равенство Ах = 0 означает, что х = 0.

Задача решения линейных уравнений состоит в нахождении такого

Вектора х по данной матрице А и данному вектору b, для которого выполняется равенство Ах = b.

Если матрица невырожденная (тем самым А – квадратная), то эта задача всегда имеет единственное решение для любого b.

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

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

Система уравнений называется однородной, если правая часть равна нулю, например, Ах = 0.

11.1. Особенности решения матричных уравнений.

Алгоритм Гауссова исключения.

Система уравнений, которую надо решить есть Ах = b.

Матрица А записывается в виде A = LU, где L и U нижняя и верхняя треугольные матрицы. Тогда:

LUx = b , которое преобразуется в , после чего решаются два уравнения с треугольными матрицами.

В алгоритме Гауса процесс исключения состоит из n-1 шагов

А^(K+1) = L^(K)A^(K)

и начинается с матрицы А⁽¹⁾ =А

Элементы матрицы А^(K+1) находятся по формуле:

a_ij^(K+1)=a_ij^(K) - Кп_, где Кп = a_iK^(K) a_Ki^(K)/ a_KK^(K) ;

Конечным результатом исключения является верхняя треугольная матрица U=A^(U) , а весь процесс эквивалентен вычислению треугольного разложения:

A = LU, где L = L_n , L_n-1, .......L₂, L₁ , а U = А⁽ⁿ⁺¹⁾

Чтобы разложение удалось вычислить, необходимо, чтобы все знаменатели a_KK^(K) в формулах были отличны от нуля. Кроме того, для обеспечения приемлемой устойчивости вычислений желательно, чтобы поправочные коэффициенты Кп были не слишком велики. В результате гауссова исключения приходим к уравнению:

LUx = b , которое преобразуется в , после чего производится обратная подстановка. Повторим: U и L - верхняя и нижняя треугольные матрицы.

Ясно, что вычисления упростятся, если одна или более из величин формулы процесса факторизации

a_ij^(K+1)=a_ij^(K) - a_iK^(K) a_Ki^(K) / a_KK^(K), входящих в правую часть, будут равны нулю. Отметим этот немаловажный факт.

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

Исключение Гаусса с частичным выбором ведущего элемента

для решения системы Ах=b

Цикл по столбцам А

For K=1 to U-1 do

Поиск максимального по модулю элемента в столбце

i_max = индекс строки такой, чтоa_jmax,K= maxa_iK i K

Проверка на нуль ведущего элемента (на вырожденность)

BUFFER - машинно-зависимая константа, оценивающая величину ошибок округления

IFa_jmax,K BUFFER, то перейти на конец цикла К.

Перестановка строк с индексами К и IMAX

For j=1 to n переставить a_Kj и a_imax

Перестановка соответствующих правых частей

Переставить b_K и b_imax

Цикл по строкам

For i=K+1 to U do

Вычисление множителей и запоминание их в столбце К

m=a_iK - a_iK/a_KK

Цикл по элементам строки I

For j=K+1 to N do

a_ij=a_ij - m _* a_Kj

Вычисление правой части строки I

b_i=b_i - m _* b_K

Конец цикла по i.

Конец цикла по K.

Обратная подстановка

Цикл по строкам в обратном порядке

For i=N to 1 do

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

Проверка на нуль ведущего элемента

IF(a_ij> BUFFER) then x_i=r / a_ij

Проверка на нуль правой части

Присвоить нуль компоненте решения и печать сообщения

Else IFr  BUFFER then x_i=0

или печать сообщения об ошибке

else печать несовместимая система, А - вырождена.

Конец цикла по i.

Матрицы L и U хранятся на месте матрицы А. Каждый множитель помещен на месте обнуляемого элемента А.

В обычной записи:

a₁₁X₁+a₁₂X₂+a₁₃X₃+...+a_1nX_n=b₁

a₂₁X₁+a₂₂X₂+a₂₃X₃+...+a_2nX_n=b₂

a₃₁X₁+a₃₂X₂+a₃₃X₃+...+a_3nX_n=b₃

......................................................

......................................................

a_n1X₁+a_n2X₂+a_n3X₃+...+a_nnX_n=b_n

......,

a₂₁X₁+a₂₂X₂+a₂₃X₃+...+a_2nX_n=b₂

a₁₁ X₁+a₁₂ X₂+a₁₃ X₃+...+a_1n X_n=b₁

Масштабирование и анализ машинного нуля.

Алгоритм Гаусса требует анализа ситуаций, когда некоторые получаемые числа равны 0.

С математической точки зрения было бы правильно положить BUFFER=0, однако на практике этого сделать нельзя. Числа которые при использовании такой арифметики должны быть нулями, почти наверняка не равны нулю при использовании арифметики с плавающей точкой из-за ошибок округления. Чтобы определить, является ли число нулем, необходимо рассмотреть два вида масштабов:

масштаб представления чисел в исходной задаче и

масштаб представления чисел в машине.

Машинный масштаб определяется точностью (числом оперируемых цифр) самой машины. Поэтому для компьютера, который выполняет действия с 10 десятичными цифрами, ведущий элемент 1,210 ^(-10) может рассматриваться как результат влияния ошибок округления и по праву может быть положен равным нулю.

Однако это было бы абсолютно неправильным, если все элементы матрицы А имели бы порядок 10^(-8). Таким образом целесообразно масштабировать систему так, чтобы все элементы были порядка 1 и тогда было бы правильным считать, что малые числа порядка 1,210^(-10) представляют собой помехи из-за ошибок округления.

Существует два естественных способа масштабирования системы Ax=b.

a) Мы можем умножить уравнения на некоторый множитель (например, на 10⁸, если элементы а_ij имеют порядок 10^(-8)) для выравнивания коэффициентов.

b) Мы можем умножить также столбцы матрицы А: это соответствует изменению единиц измерения неизвестных (например, переходу от сантиметров к метру).

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

Например:

К счастью, большинство задач могут быть масштабированы.

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

Наиболее простой способ масштабирования состоит в таком умножении каждой строки (уравнения), чтобы наибольший элемент строки был равен 1.

В предположении, что это сделано, перейдем к выбору величины BUFFER.

Если машина имеет t разрядов с основанием b, то в качестве значения BUFFER может быть всего С*b^(-t), где С - некоторое небольшое число (3  4) и должно постепенно расти с увеличением порядка матрицы.

Над каждым элементом а_ij выполняется самое большое 2n² арифметических операций, поэтому С=n² - самое большое значение, которое следует рассматривать. Вообще C=n достаточно надежный выбор.

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

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

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

3. На практике обычно масштабируют строки (делением каждого уравнения на его наибольший коэффициент). Это достаточно надежный способ для обычно встречающихся задач.

4. Для получения LU - разложения исключением Гаусса требуется порядка операций сложения и умножения.

Отметим дополнительно следующее:

б) формулами процесса факторизации пользуются,

если a_iK^(K) 0 и a_Kj^(K) 0.

в) новый элемент возникает, если a_ij^(K) =0 , а a_ij^(K+1) 0

Последнее явление называют заполнением.

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

Поясним эффект заполнения, возникающий в разложении Гаусса.

Пример 1.

Дана система уравнений.