Описанные выше методы дают возможность точно рассчитать собственные значения матриц, но при исследовании матриц большой размерности становится ясно, что данные методы неприменимы, так как характеристическое уравнение данных матриц будет обладать слишком большой вычислительной сложностью, поэтому в данном случае прибегают к итерационным методам. Эти методы находят приближение к собственным значениям, не используя характеристическое уравнение. Они обладают достаточно высокой точностью и упрощают вычислительный процесс, поэтому рассмотрим их подробнее.
1.3 Итерационные методы нахождения собственных значений
1.3.1 Степенной метод
Степенной метод обычно используется для приближенного вычисления крайних собственных значений, в данном случае рассмотрим нахождение максимального по модулю числа. Рассмотрим данный метод подробнее.
Если матрица имеет размеры , то у неё обязательно имеются собственных чисел, при этом данные собственные числа не обязательно должны быть различны, и, предположим, что матрица имеет собственных векторов и полную систему из собственных векторов. Тогда это даёт возможность разложить любой вектор в линейную комбинацию собственных векторов (1.14). Далее, умножая уравнение (1.1) на исходную матрицу, получим выражение (1.15), проведя данное умножение k раз получаем равенство (1.16) и выполнив его преобразования получим (1.17):
(1.14)
где - собственные вектора.
(1.15)
(1.16)
(1.17)
Таким образом, воспользовавшись всеми введёнными утверждениями и формулами можно сделать вывод, что данным метод решения частичной задачи нахождения собственного значения реализуется с помощью алгоритма:
1) выбирается или случайно генерируется вектор начального приближения;
2) по формулам (1.18) и (1.19) строится приближение вектора:
(1.18)
(1.19)
где не противоречит условию
3) Итерационно строим приближения, пока не будет достигнуто условие сходимости (1.20):
(1.20)
где - заданная точность;
4) при достижении условия сходимости полученное приближение будет являться собственным вектором, используя который, можно найти собственное значение по формуле (1.21):
(1.21)
поскольку .
1.3.2. Метод скалярный произведений
Для нахождения первого собственного значения действительной матрицы можно указать несколько иной итерационный процесс, являющийся иногда более выгодным. Метод основан на образовании скалярных произведений.
(1.22)
где - матрица, транспонированная с матрицей ,
- выбранный каким-либо образом начальный вектор.
Пусть - действительная матрица и - ее собственные значения, которые предполагаются различными, причем
(1.23)
Возьмем некоторый ненулевой вектор и с помощью матрицы построим последовательность итерации
(1.24)
Для вектора образуем также с помощью транспонированной матрицы вторую последовательность итерации
(1.25)
где
В пространстве выберем два собственных базиса и соответственно для матриц и , удовлетворяющих условиям биортонормировки: .в базисе - через , т. е.
(1.26)
отсюда
(1.27)
(1.28)
Составим скалярное произведение
(1.29)
таким образом,
(1.30)
1.3.3 Метод вращений Якоби
Данный метод является одним из наиболее часто используемых при решения полной задачи нахождения собственного значения при условии, что используемая матрица симметрична.
Дадим определение матрице вращения. Ортогональная матрица
(1.31)
порождает преобразование поворота на угол в двумерной плоскости. Матрица , отличающаяся от единичной только лишь четырьмя элементами (1.32) называется матрицей вращения:
(1.32)
где - заданные целые числа.
В данном случае - ортогональная матрица. Она порождает преобразование поворота на угол в двумерной плоскости, натянутой на векторы канонического базиса с номерами.
Далее опишем основную идею метода вращений Якоби. Как уже говорилось матрица - симметричная матрица. Образуем матрицу , столбцами которой будут ортонормированные собственные векторы . для данной матрицы справедлива формула (1.33)
(1.33)
где -- диагональная матрица с элементами на диагонали.
В методе вращения Якоби матрица строится как предел последовательности ортогональных матриц так, что:
(1.34)
причем при каждом матрица конструируется как произведение матриц вращения.
Образуем по исходной матрице матрицу , и подбираем параметры матрицы вращения, то есть значения , так, чтобы матрица была максимально близка к диагональной. Опуская простые вычисления, представим формулу для вычисления суммы квадратов внедиагональных элементов матрицы :
(1.35)
Далее, определим числа k,l из условия (1.36):
(1.36)
и затем угол ц из условия (1.37)
(1.37)
(1.38)
Сумма квадратов внедиагональных элементов матрицы принимает наименьшее значение, если использовать данное условие при выборе параметров для матрицы вращения.
Опишем и сам алгоритм метода вращений Якоби. Пусть = A. Образуем последовательность матриц , при помощи рекуррентной формулы (1.39):
(1.39)
где параметры матрицы вращения задаются так, что сумма квадратов внедиагональных элементов матрицы уменьшается, то есть используя формулы вида (1.27-1.29). Продолжаем итерации до тех пор, пока все внедиагональные элементы матрицы не станут достаточно малыми.
При достижении достаточно мылах значений в качестве приближений к собственным числам матрицы принимаются диагональные элементы матрицы , а столбцы матрицы , считаются приближениями к собственным векторам матрицы .
1.3.4 Метод QR
Далее приведем описание метода QR. Данный алгоритм использует плоские преобразования вращений Гивенса. Матрица, определяющая эти преобразования имеет вид . Данная матрица имеет структуру, похожею на матрицу плоских вращений Якоби , только здесь двумерная подматрица из элементов, стоящих на пересечении и ?? строк и столбцов, используется в виде (1.40):
(1.40)
Данные числа ?? и ?? связываются соотношением . Первый полный шаг преобразования Гивенса, применяемого к матрице Хессенберга ??-го порядка в рамках ???? метода, состоит из элементарных подшагов, основной целью которых является последовательное обнуление всех поддиагональных элементов в столбцах от первого до -го [7]. В результате получим разложение матрицы Хессенберга в виде произведения ортогональной и треугольной матриц.
Для определения и на первом промежуточном шаге, используется произведение матриц (1.41):
(1.41)
При этом и берутся такими, что . Это означает, что в результате первого промежуточного шага матрица , полученная при использовании таких и , не будет содержать ненулевых элементов под диагональю в первом столбце.
Аналогично совершается второй промежуточный шаг: матрица получается из предыдущей с помощью матрицы Гивенса , отличающейся от тем, что подматрица смещается на одну позицию вдоль диагонали и угол поворота подбирается так, чтобы в матрице обнулить элементы .
Продолжив далее преобразования Гивенса, в итоге получим правую треугольную матрицу:
(1.42)
Последнее равенство можно переписать в виде:
(1.43)
который позволяет считать выполненным требуемое при разложение
(1.44)
где - ортогональная, а - правая треугольная матрицы. При этом матрица
(1.45)
являющаяся результатом первого полного шага QR метода, сохраняет не только спектр данной матрицы, но и форму Хессенберга, благодаря чему приведение исходной матрицы к почти треугольному виду достаточно сделать только один раз [7].
Очевидно, скалярные параметры и матриц Гивенса , благодаря которым происходит переход от матрицы Хессенберга через матрицы Хессенберга к матрице Хессенберга , можно вычислять на ??-ом промежуточном шаге () по формулам:
(1.46)
(1.47)
Также можно модифицировать данный метод, уменьшив количество итераций, получая метод называемый QR со сдвигом.
Алгоритм со сдвигами состоит в следующем. Положить и выполнить:
1. выбрать сдвиг s вблизи некоторого собственного значения ;
2. вычислить по схеме QR разложение для матрицы
3. определить матрицу
4. если сходимость к собственным значениям не достигнута, положить и перейти к пункту 1.
Матрицы ) и в QR алгоритме со сдвигами будут ортогонально подобными:
(1.48)
Если - точное собственное значение матрицы ), то QR-итерация сойдется за один шаг.
Если не является точным собственным значением, то элемент считается сошедшим, если левый нижний блок ) достаточно мал. Поскольку матрица получается из матрицы ортогональным подобием, то включает в себя ошибки округлений порядка . Таким образом, любой поддиагональный элемент в , по абсолютной величине меньший, чем , мог быть нулем, поэтому его можно заменить на ноль. То есть при условии можно положить
Также следует упомянуть что не обязательно брать верхнюю треугольною матрицу R, точно также можно выбрать нижнюю треугольную. Это приведёт к QL алгоритму, поскольку данное разложение будет представлено в виде A=QL
Если сначала привести диагональную матрицу к трехдиагональной форме и использовать QL-алгоритм, это даст меньшую ошибку округления, поэтому именно он и применяется на практике в большинстве случаев и выбран для дальнейшей реализации.
1.4 Выводы по главе 1
Таким образом, в первой главе данной работы была описана поставленная задача для исследования, приведены основные термины, используемые в ходе рассмотрена задачи.
Для рассматриваемой проблемы представлены основные методы её решения, а именно два точных метода - метод Данилевского и Леверрье-Фадеева, а также четыре итерационных - степенной метод, метод скалярных произведений, метод вращений Якоби и метод QR. Помимо этого для метода QR рассмотрена его модификация - метод QL со сдвигом.
Учитывая специфику исследования, можно утверждать, что точные методы для решения данной задачи не подходят, так как решения характеристического уравнения для матриц большой размерности имеют очень большую вычислительную сложность.
Поэтому для программной разработки в ходе данной работы выбраны степенной метод, метод вращений Якоби и QL метод со сдвигом. В дальнейшем разработанные методы будут использованы для проведения экспериментов и сбора данных, необходимых для сравнительного анализа методов нахождения собственных значений симметричных матриц большой размерности.
2. Программная реализация методов нахождения собственных значений
2.1 Структура программы
Как уже было сказано, для сравнения и реализации были выбраны три основных метода нахождения собственных значений: степенной метод, метод вращений Якоби, и метод QL со сдвигом. Данные методы реализованы на языке C++ в среде Qt Creator.
Qt Сreator - это кроссплатформенная свободная среда программирования для разработки на языках С и С++. Разработана для работы с фреймворком Qt и включает в себя графический интерфейс отладчика и визуальные средства разработки интерфейса.
Основа структуры программы включает в себя:
1) подключенные библиотеки (iostream, math, ctime);
2) функция генерации матрицы заданной размерности, использующая генератор псевдослучайный чисел rand;
3) функция вывода матрицы;
4) функция, реализующая степенной метод;
5) функция, реализующая метод вращения Якоби;
6) функция, реализующая преобразование Хаусхолдера, используемая для реализации метода QL;
7) вспомогательная и основная функции, реализующие метод QL со сдвигом;
8) основная часть;
9) интерфейс программы;
Также в коде предусмотрена возможность замера времени выполнения алгоритма и количества итераций, необходимых для дальнейшего сравнительного анализа методов.
Работа программы начинается с ввода пользователем размерности матрицы. Учитывая, что рассматриваемые матрица должна быть квадратной и симметричной, вызывается функция Gen_matrix, которая, используя генератор псевдослучайных чисел создает симметричную квадратную матрицу заданной размерности заполняя её числами в диапазоне от 0 до 1. Пример сгенерированной матрицы приведен на рисунке 2.1.
Рисунок 2.1 - Пример сгенерированной матрицы
После этого пользователь выбирает один из трех методов нахождения поиска собственных значений представленной в программе функция Pow_Meth, Jacobi_Meth, QL_Meth, каждая из который принимает на вхож матрицу и её размерность, после чего выводится матрица с помощью функции Print_Matrix, выводятся собственные значения матрицы и данные для анализа, а именно время выполнения алгоритма, количество итераций и максимальное собственное значение матрицы. Блок-схема главной функции программы приведена на рисунке 2.2.
Таким образом можно сказать, что программа принимает на вход только размерность матрицы и выбранный метод нахождения собственных значений, а на выходе выдает собственные значения и данные для анализа. Рассмотрим реализация каждого из методов подробнее