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

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

1) Решение систем линейных алгебраических уравне­ний (СЛАУ).

2) Вычисление определителей квадратной матрицы А.

3) Для данной квадратной матрицы А вычисление обрат­ной А^-1.

4) Определение собственных значений и собственных векторов квадратной матрицы А.

3.1. Нормы векторов и матриц

Определение 3.1. Нормой вектора называется неотрицательное число, которое обозначается символом || || и удовлетворяет следующим условиям:

1) || || > 0 при ¹ и || || = 0;

2) || || = | c |×|| || при любом числовом множителе;

3. || || ≤ || || + || ||. Это соотношение называют неравенством треугольника.

Норму вектора можно ввести различными способами. Наиболее часто для векторов n-мерного арифметического пространства

= (x₁, x₂, ..., x_n)^T используются следующие нормы:

1) || ||₁ = , i = 1, 2, ..., n (3.1)

(называется кубической);

2) || ||₂ = (3.2)

(называется октаэдрической);

3) || ||₃ = = (3.3)

(называется сферической, порождена скалярным произведением и определяет длину вектора ).

Норма (3.3) порождена скалярным произведением которое выражается формулой: .

Для скалярного произведения векторов справедливы соотношения:

= = .

Если A симметричная матрица, то = .

Определение 3.2. Если в пространстве векторов введена норма || ||, то согласованной с ней нормой в пространстве матриц называют норму

. (3.4)

Согласованные с нормами векторов (3.1)-(3.3) нормы матриц определяются формулами

, (3.5)

, (3.6)

. (3.7)

В формуле (3.7) – собственные значения матрицы A^TA, которая является симметрической.

Формула (3.7) следует из того, что для симметрической матрицы B можно доказать справедливость соотношения:

= , (3.8)

где l_i – собственные значения матрицы B.

Определение 3.3. Две нормы || ||_α и || ||_β называются эквивалентными, если существуют такие постоянные γ₁ и γ₂, что для всех векторов справедливы соотношения:

, и .

Нормы || ||₁, || ||₂ и || ||₃ эквивалентны между собой, так как выполнятся неравенства

|| ||₁ ≤ || ||₃ ≤ || ||₂ ≤ n×|| ||₁.

Определение 3.4. Будем говорить, что последовательность векторов сходится к вектору по данной норме || ||, если выполняется соотношение = 0.

Из эквивалентности норм || ||₁, || ||₂ и || ||₃ следует, что если последовательность векторов сходится по одной из этих норм, то она сходится и по остальным нормам.

Договоримся в дальнейшем под нормой || || подразумевать одну из указанных выше норм, а при необходимости конкретизировать, какую именно.

При этом под нормой матрицы будем понимать норму, согласованную с нормой матрицы.

3.2. Решение систем линейных алгебраических уравнений (СЛАУ)

Теоретические условия существования и единственно­сти решения систем линейных уравнений известны – главный определитель не должен быть равен нулю. Тог­да решение можно найти по правилу Крамера, или мето­дом исключения неизвестныхГаусса. Метод Гаусса и правило Крамера относятся к прямым методам решения систем линейных алгебраических уравнений. Они позво­ляют за конечное число действий получить точное реше­ние системы, при условии, что все действия выполняют­ся точно, без округления. Но на практике, при больших порядках системы, правило Крамера требует слишком много времени для вычисления определителей. Если оп­ределители вычислять формально по определению как сумму п!слагаемых, то число операций имеет порядок п!п.Правило Крамера используется чаще для теоретичес­ких исследований, а на практике почти не применяется.

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

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

а) Метод Гаусса для решения систем линейных уравнений

Прямой ходметода Гаусса преобразует систему (3.9) к треугольному виду исключением соответствующих неиз­вестных. Пусть a₁₁ ¹ 0 . Первый шаг заключается в исклю­чении переменной x₁с помощью первого уравнения из ос­тальных уравнений. Разделим первое уравнение на a₁₁:

; . (3.10)

Затем от второго уравнения вычтем первое уравне­ние, умноженное на a₂₁. В результате, на месте второго уравнения получим уравнение, не содержащее х₁. Чтобы исключить х₁ из третьего уравнения, вычтем из первое уравнение, умноженное на a₃₁. Аналогично ис­ключаем х₁из четвертого и последующих уравнений. Для исключения х₁из i-го уравнения (i = 2, 3, ... , п)приме­ним формулы

; . (3.11)

В результате этих вычислений получим систему вида

(3.12)

На втором шаге исключаем переменную х₂с помощью второго уравнения из третьего и последующих уравне­ний. Предположим, что ¹ 0. Разделим второе уравне­ние на :

; . (3.13)

Затем в системе (3.12) с помощью второй строки исключим

х₂ из i-го уравнения (i = 3, 4, ... , n), применяя формулы:

; . (3.14)

Система (3.12) преобразуется к следующему виду:

(3.15)

1. В общем случае на шаге т для т = 1, 2,..., n - 1,

делим сначала т-е уравнение на :

; , (3.16)

а затем исключаем переменную x_m c помощью m-го уравнения из i-го, где

i = m +1,..., n:

; . (3.17)

Здесь предполагается, что на каждом шаге выполня­ется условие ¹ 0.

В результате (n - 1)-го шага система (3.9) приобретает вид

(3.18)

2. Обратный ход метода Гаусса вычисляет неизвест­ные x_i в обратном порядке. Из последнего уравнения в (3.18) находим

х_п = (3.19)

Неизвестные x_i определяем по следующим формулам:

x_i = i = n - 1,n - 2,..., 1.(3.20)

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

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

Метод Гаусса с выбором главного элемента по стол­бцамотличается от алгоритма (3.16) – (3.20) только тем, что перед преобразованием (3.16) нужно выполнить по­иск максимального по модулю элемента в т-мстолбце и переставить строки системы уравнений так, чтобы мак­симальный элемент стал диагональным элементом матри­цы коэффициентов.

б) Алгоритм метода Гаусса с выбором главного элемента

По столбцам

1. Для т = 1, 2, ... , п - 1 выполним преобразования: Найдем максимальный по абсолютной величине эле­мент в т-м столбце. Пусть это будет элемент a_im.Если i ¹ т, то меняем местами i-ю и т-юстроки:

r = , = , = r , j = 1, ... , п;

r = b_i, b_i = b_т, b_т = r.

Элементы матрицы и вектора-столбца свободных членов после преобразования на т-м шаге обозначим , , причем = , = b_i.

Делим т-е уравнение на ведущий элемент :

; ,

затем исключаем переменную х_т с помощью m-го урав­нения из i-го, где

i = т + 1,..., п:

; .

2.Вычисляем неизвестные x_i в обратном порядке:

х_п = ; x_i = i = n - 1,n - 2,..., 1.

Приведенный вариант метода Гаусса дает решение и тогда, когда обычный метод Гаусса перестает работать из-за деления на 0.

При реализации метода Гаусса на каком-либо языке программирования удобно использовать исходные матри­цу a и вектор bдля хранения промежуточных результатов преобразований. Приведенные формулы преобразова­ний записываются как операторы присваивания, т. е. имена переменных a_ij и b_i записываются без верхних ин­дексов.

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

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

Отметим, что последний вари­ант с выбором главного элемента считается лучшим.

Общее число операций для решения системы уравне­ний методом Гаусса составляет приблизительно 2n³/3 + 2n², при этом большая часть, т. е. 2п³/3, операций при­ходится на прямой ход.

в) Итерационный метод

х = Вх + с. (3.22)

Тогда метод простой итерации определяется формулой:

x^k+1 = Bx^k + c, k = 0, 1, ... , (3.23)

где начальное приближение x⁰ задано.

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

||x^k+1 - x^k|| £ ε.

Сформулируем теоремы об условиях сходимости мето­да простых итераций.

Теорема 3.1 (достаточное условие сходимости).Если ||В|| < 1, то система (3.21) имеет единственное решение, а итерационный метод (3.23) сходится к решению со ско­ростью геометрической прогрессии.

Теорема 3.2. (необходимое и достаточное условие схо­димости).Пусть система (3.22) имеет единственное реше­ние. Итерационный процесс (3.23) сходится к решению системы (3.22) тогда и только тогда, когда все собствен­ные значения матрицы В по модулю меньше единицы.

Теоремы 3.1 и 3.2 не дают способов преобразования произвольной системы линейных уравнений к виду (3.22) так, чтобы метод итераций (3.23) при этом сходился к решению. Эти теоремы имеют важное теоретическое значение. На практике для проверки условия сходимости метода итераций более полезна теорема 3.1, а теорему 3.2 удается использовать тогда, когда собствен­ные значения матрицы В известны. Задача определения собственных значений матрицы в общем случае сложнее, чем задача решения системы линейных уравнений.

Для преобразования системы уравнений к виду, удоб­ному для итераций, можно умножить систему (3.21) на матрицу D = А^-1- ε, где ε – произвольная матрица. Тог­да система примет вид (3.22)

х= Вх + с, В = εА, с = Db (3.24)

и матрица В будет удовлетворять теореме 3.1, если выб­рать элементы матрицы ε достаточно малыми по модулю. Однако этот прием не выгоден, так как для вычисления обратной матрицы А^-1необходимо выполнить не меньше операций, чем при решении самой системы.

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

Далее будут рассмотрены разностные методы решения краевых задач для обыкновенных дифференциальных уравнений и уравнений в частных производных, которые приводят к решению систем линейных уравнений с боль­шим числом неизвестных. Учитывая специфику таких систем, часто удается построить эффективные итерацион­ные методы их решения.

г) Метод Зейделя

Теорема 3.3(достаточные условия сходимости).Пусть при всех i для коэффициентов системы уравнений (3.25) выполняются условия

. (3.27)

Тогда метод Зейделя сходится и выполняется неравен­ство

||x^k+1 - x^k||₁ £ q||x^k – x*||₁ (3.28)

где х* – точное решение системы (3.25).

Если матрица А удовлетворяет условию (3.27), гово­рят, что А – матрица с диагональным преобладанием.

Теорема 3.4(достаточные условия сходимости).Пусть матрица А системы (3.1) – вещественная, симмет­ричная и положительно определенная. Тогда метод Зей­деля сходится.

3.3.Погрешность решения и обусловленность системы уравнений

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

Ах = b₁,b₁ = b + η.

Пусть х¹ решение неточно заданной системы Ах = b₁, х – решение точной системы Ах = b.Обозначим погрешность решения через r = х¹ - х.Тогда можно запи­ть Ах¹= b₁в виде А (х + r) = b + η,и Аr= η.

Определение 3.5. Мерой обусловленности системы называется число

. (3.29)

Мера обусловленности системы равна верхней грани отношения относительной погрешности решения к относительной погрешности правой части. Из формулы (3.29) следует неравенство

(3.30)

Если мера обусловленности системы принимает боль­шое значение, это значит, что небольшая погрешность правой части может привести к большой погрешности решения, т. е. полученное приближенное решение ока­жется непригодным.

Учитывая, что r = А^-1η,можно получить формулу вы­числения меры обусловленности системы:

(3.31)

Определение 3.6. Мерой обусловленности матрицы А называется число

. (3.32)

Для вычисления меры обусловленности матрицы мож­но с помощью (3.31) получить формулу

. (3.33)

Учитывая (3.30), можно записать

(3.34)

Неравенство (3.34) связывает относительные погреш­ности правой части и решения системы через свойства матрицы системы.

Определение 3.7. Системы уравнений и матрицы назы­ваются плохо обусловленными,если их меры обусловлен­ности принимают большие значения, и хорошо обуслов­ленными,если меры обусловленности принимают малые значения.

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

3.4. Вычисление определителя и обратной матрицы

Теорема 3.5. Собственные значения симметричной матрицы с действительными элементами действительны, а собственные векторы, соответствующие различным соб­ственным значениям, взаимно ортогональны.

Теорема 3.6.Если λ_min и λ_mах – наименьшее и наиболь­шее собственные значения действительной симметричной матрицы А, то для любого вектора хсправедливо неравен­ство

λ_min(x, х)£(Ах, х)£ λ_max(х, х).(3.37)

Определение 3.9.Действительная симметричная мат­рица А называется положительно определенной, если для любого вектора х¹0 выполняется условие

(Ах, х)> 0.(3.38)

Теорема 3.7.Действительная симметричная матрица А является положительно определенной тогда и только тогда, когда все ее собственные значения положительны.

Теорема 3.8(критерий Сильвестра).Для того чтобы действительная симметричная матрица А = (а_ij) была по­ложительно определенной, необходимо и достаточно, что­бы все главные диагональные миноры ее определителя были положительны:

(3.39)

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

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

Алгоритм определения наибольшего по модулю соб­ственного значения и соответствующего собственного вектора матрицы с положительными элементами:

Алгоритм метода скалярных произведений:

1) зададим начальные приближения: х⁰ – ксобствен­ному вектору матрицы А и у⁰ = х⁰– к собственному век­тору транспонированной матрицы А^T;

k = 0;

2) вычислим (k + l)-e приближение к наибольшему собственному значению λпо формулам

x^{k +1}= Ах^k,y^{k +1} =A^Ty^k, λ_k = , k = k +1; (3.41)

3) если |λ_k+1- λ_k| ³ ε, переходим к п. 2, иначе – к п. 4;

4) конец.

7.Вычисление всех собственных значений положительно определенной симметричной матрицы

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

Пусть уже вычислены первые m собственных значений λ₁, λ₂, ..., λ_m и т соответствующих собственных векторов x₁, x₂, ..., x_m.

Алгоритм вычисления очередного (т + 1)-го собственного значения и

соответствующего собственного вектора:

0) выберем начальное приближение ; k = 0;

1) вычислим k-e приближение к собственному значе­нию λ_m+1:

; (3.42)

2) находим вектор из уравнения

; (3.43)

3) если m > 0, ортогонализируем вектор к первым т собственным векторам:

- ; (3.44)

4) нормируем полученный вектор:

; (3.45)

5) k = k+1.

Процесс 1–5 повторяется до тех пор, пока не будет выполнено условие сходимости итераций

(3.46)

где ε – заданная погрешность.

При вычислении первого собственного значения и со­ответствующего вектора п. 3 пропускается.

Этим алгоритмом можно вычислить все собственные значения и собственные векторы.