Итерационный метод вращений Якоби

Размещено на http://www.allbest.ru/

Реферат

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

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

При написании пояснительной записки было использовано 7 литературных источников.

Перечень ключевых слов: собственные значения, собственные вектора, вращение Якоби, симметричная матрица.

Содержание

• Введение
• 1. Теоретическая часть
• 1.1 Собственные значения и собственные вектора матрицы
• 1.2 Итерационный метод вращений Якоби решения симметричной полной проблемы собственных значений
• 1.3 Метод вращений Якоби
• 1.4 Решение обобщенной проблемы собственных значений
• 2. Практическая часть
• 2.1 Пример решения
• 2.2 Реализация на ЭВМ
• Заключение
• Список литературы
• Приложение 1. Листинг программы
• Приложение 2. Пример результата работы программы

Введение

Целый ряд инженерных задач сводится к рассмотрению систем уравнений, имеющих единственное решение лишь в том случае, если известно значение некоторого входящего в них параметра. Этот особый параметр называется характеристическим, или собственным значением системы. С задачами на собственные значения инженер сталкивается в различных ситуациях. Так, для тензоров напряжений собственные значения определяют главные нормальные напряжения, а собственными векторами задаются направления, связанные с этими значениями. При динамическом анализе механических систем собственные значения соответствуют собственным частотам колебаний, а собственные вектора характеризуют моды этих колебаний. При расчете конструкций собственные значения позволяют определять критические нагрузки, превышение которых приводит к потере устойчивости. Выбор наиболее эффективного метода определения собственных значений или собственных векторов для данной инженерной задачи зависит от ряда факторов, таких, как тип уравнений, число искомых собственных значений и их характер. Алгоритмы решения задач на собственные значения делятся на две группы. Итерационные методы очень удобны и хорошо приспособлены для определения наименьшего и наибольшего собственных значений. Методы преобразований подобия несколько сложней, зато позволяют определить все собственные значения и собственные векторы. В данной работе будет рассмотрен метод вращений Якоби. Однако сначала приведем некоторые основные сведения из теории матричного и векторного исчислений, на которых базируются методы определения собственных значений.

1. Теоретическая часть

1.1 Собственные значения и собственные вектора матрицы

Рассмотрим квадратную матрицу n-ого порядка:

Собственные значения _i квадратной матрицы A есть действительные или комплексные числа, удовлетворяющие условию:

,

E - единичная матрица, - собственный вектор матрицы A, соответствующий некоторому собственному значению .

Матрица называется характеристической матрицей матрицы A. Т.к. в матрице по главной диагонали стоят , а все остальные элементы равны нулю, то характеристическая матрица имеет вид:

Определитель этой матрицы называется характеристическим или вековым определителем и равен:



В развернутом виде он является многочленом n-ой степени относительно , т.к. при вычислении этого определителя произведение элементов главной диагонали дает многочлен со старшим членом , т.е.

и называется характеристическим многочленом. Корни этого многочлена - собственные значения или характеристические числа матрицы A. Числа называются коэффициентами характеристического многочлена.

Ненулевой вектор называется собственным вектором матрицы A, если эта матрица переводит вектор X в вектор

,

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

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

Определитель этой системы равен нулю, т.к. из этого условия были определены собственные значения матрицы A. Следовательно, система имеет бесконечное множество решений. Ее можно решить с точностью до постоянного множителя (как систему однородных уравнений). Решив эту систему, мы найдем все координаты собственного вектора X. Подставляя в систему однородных уравнений поочередно , получаем n собственных векторов.

1.2 Итерационный метод вращений Якоби решения симметричной полной проблемы собственных значений

матрица вектор эрмитовый итерационный

Метод вращений применим к эрмитовым матрицам с комплексными элементами. Для простоты рассмотрим частный случай эрмитовых матриц - действительных симметричных.

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

, (1)

где - некоторая ортогональная матрица; - диагональная матрица, элементами которой являются собственные числа матрицы .

Так как для ортогональной матрицы обратная матрица совпадает с транспонированной матрицей: , то равенство (1) равносильно следующему

. (2)

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

, (3)

где - некоторая ортогональная матрица, которая несложно вычисляется; - квадратная матрица с близкими к нулю недиагональными элементами, так как найти в общем случае трудоемко:

. (4)

Естественно предполагать, что коэффициенты , , будут близки к собственным числам матрицы .

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

. (5)

За меру близости матрицы к диагональному виду примем число:

. (6)

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

Таких процессов может быть построено большое число. Рассмотрим метод вращений. Он достаточно прост по вычислительной схеме и обладает быстрой сходимостью.

1.3 Метод вращений Якоби

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

(7)

Матрица вращений получается из единичной матрицы заменой ее элементов на пересечениях - х и - х строк и столбцов записанными элементами.

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

. (8)

Введем обозначения:

.

Ввиду особой структуры матрицы (7) все столбцы матрицы , кроме - го и - го, будут такими же, как и в матрице . Элементы столбцов номеров и будут вычисляться по формулам:

;

. (9)

Аналогично строки матрицы , кроме - й и - й, будут такими же, как в матрице , а элементы ее строк - й и - й вычисляется по формулам:

;

. (10)

Равенства (9) и (10) позволяют несложно вычислить коэффициенты

,

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

. (11)

Сейчас назначим в выражении (11) такое значение угла , чтобы рассматриваемый элемент обратился в нуль. Из этого условия получаем выражение:

,

из которого находим значение угла

(12)

Убедимся, что полученная таким образом матрица ближе к диагональному виду, чем предшествующая . Для этого вычислим - меру близости матрицы к диагональному виду. Пусть матрица имеет меру ее близости к диагональному виду . Матрица симметрична и на основании выражений (9) - (11) можно показать, что

, (13)

так как есть наибольший по модулю недиагональный элемент, отличный от нуля, который исключен при получении матрицы .

Тогда на основании выражения (13) верно неравенство:

.

Следовательно, мера уменьшиться при преобразовании (8).

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

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

и, следовательно,

. (14)

С помощью этого неравенства из выражения (13) получается

, (15)

где

.

Пусть значение . Тогда величина и из неравенства (15) вытекает цепочка неравенств

.

Так как , то для меры близости матрицы к диагональному виду будем иметь оценку

Отсюда следует, что со скоростью не меньшей, чем скорость сходимости геометрической прогрессии со знаменателем .

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

Существуют и другие методы выбора наибольшего по модулю элемента, предназначенного для исключения. Например, можно выбрать наибольшую из сумм (5). Если это есть , то в качестве элемента можно выбрать наибольший по модулю элемент из найденной строки . При таком правиле выбора элемента необходимо выполнить примерно операций. Отметим, что оценка (16) для такого метода выбора наибольшего по модулю элемента, предназначенного для исключения, сохраняется.

1.4 Решение обобщенной проблемы собственных значений

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

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

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

Решение ряда задач сводится к обобщенной проблеме собственных значений с симметричными действительными положительно определенными матрицами , :

. (1)

Тогда все собственные числа обобщенной проблемы положительны (неотрицательны).

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

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

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

Рассмотрим первый метод. С помощью преобразования подобия представим в уравнении (1) матрицу в виде разложения на множители:

, (3)

где - ортогональная матрица; - диагональная матрица, ненулевые элементы которой равны: , ; , - собственные числа матрицы , при этом , ; ; - диагональная матрица, ненулевые элементы которой равны: , .

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

Матрицы и могут быть найдены, например, с помощью процедуры метода вращений, изложенной ранее.

Согласно формуле (3) подставим матрицу в уравнение (1), получим

.

Умножая полученное равенство слева на матрицу , получим:

. (4)

Представим

(5)

и, подставляя это выражение в уравнение (4), приведем его сначала к виду

,

а затем запишем так

. (6)

Уравнение (6) можно представить в стандартном виде

, (7)

где - симметричная матрица.

Так как в уравнении (7) матрица симметричная, то для нахождения ее собственных чисел и собственных векторов можно второй раз применить метод вращений. В результате применения метода вращений получаются множители разложения для матрицы в виде:

(8)

где - ортогональная матрица; - диагональная матрица.

На основании выражений (7), (8) запишем равенство

,

из которого следует

, (9)

где - ортогональная матрица, равная произведению обратимых матриц.

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

Собственные векторы проблемы (7) связаны с векторами исходной проблемы (1) соотношением (5).

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

, (10)

где - правая треугольная матрица разложения по схеме метода квадратных корней (разложение Холецкого, см. формулы (11) подраздела).

Подставляя выражение (10) для матрицы в уравнение (1) и умножая его слева на матрицу , получим

.

Затем, подставляя в это уравнение , получим

, (11)

где - симметричная матрица.

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

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

2. Разработаны методы уточнения приближенного решения полной проблемы собственных значений.

3. Разработаны методики ускорения сходимости при решении проблем собственных значений как частичных, так и полной.

2. Практическая часть

2.1 Пример решения

Найти собственные значения и собственные вектора для матрицы

1. Положим .

2. Выделим максимальный по модулю элемент в над диагональной части: . Так как , то процесс продолжается.

3. Находим угол поворота:

4. Сформируем матрицу вращения:

5. Выполним первую итерацию:

Положим и перейдем к пункту 2.

2(1). Максимальный по модулю над диагональный элемент

Так как , процесс продолжается.

3(1). Найдем угол поворота:

4(1). Сформируем матрицу вращения:

5(1). Выполним вторую итерацию:

Положим и перейдем к пункту 2.

2(2) Максимальный по модулю над диагональный элемент

.

3(2). Найдем угол поворота:

4(2). Сформируем матрицу вращения

5(2). Выполним третью итерацию и положим и перейдем к пункту 2:

2(3). Максимальный по модулю над диагональный элемент

3(3). Найдем угол поворота:

4(3). Сформируем матрицу вращения:

5(3). Выполним четвертую итерацию и положим и перейдем к пункту 2:

2(4). Так как , процесс повторяется

3(4). Найдем угол поворота

4(4). Сформируем матрицу вращения:

5(4). Выполним пятую итерацию и положим и перейдем к пункту 2:

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

Собственные значения:

Для нахождения собственных векторов вычислим

отсюда

или после нормировки

2.2 Реализация на ЭВМ

Алгоритм метода вращений

1) Положить и задать .

2) Выделить в верхней треугольной над диагональной части матрицы максимальный по модулю элемент .

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

.

Собственные векторы находятся как i-e столбцы матрицы, получающейся в результате перемножения:

Если , процесс продолжается.

3) Найти угол поворота по формуле.

4) Составить матрицу вращения .

5) Вычислить очередное приближение

Положить и перейти к пункту 2.

Замечания:

Используя обозначение

,

можно в пункте 3 алгоритма вычислять элементы матрицы вращения по формулам

,

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

Заключение

При выполнении данной курсовой работы был изучен метод вращения Якоби решения симметрически полной проблемы собственных значений.

Этот метод при небольших размерах матриц дает неплохие результаты, и позволяет с большой точностью вычислять собственные пары. При больших размерах матриц реализация метода наталкивается на существенные потери машинных ресурсов (из-за необходимости поиска максимального элемента матрицы, и перемножения матриц большой размерности). Так же резко сужает возможность применить данный метод на практике то, что он может быть применен только к симметрическим, вещественным матрицам. Ввиду вышесказанного, этот метод на практике применяется редко, чаще применяется циклический метод Якоби с барьерами. Скорость сходимости которого так же асимптотически квадратичная.

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

Бахвалов Н.С. Численные методы. М.: Наука, 1987. 600 с.

Бусленко H. П., Голенко Д.И., Соболь И.М., Срагович В.Г., Шрейдер Ю.А.. Метод статистических испытаний (метод Монте-Карло). М.: Физматгиз, 1962. 332 с.

Бусленко Н.П., Шрейдер Ю.А. Метод статистических испытаний (Монте-Карло) и его реализация на цифровых вычислительных машинах. М.: Государственное издательство физико-математической литературы, 1961. 226 с.

Демидович Б.П., Марон И.А. Основы вычислительной математики. М.: Наука, 1966. 664 с.

Ермаков С.М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1975. 294 с.

Калиткин Н.Н. Численные методы. М.: Наука, 1978. 512 с.

Копчёнова Н.В., Марон И.А. Вычислительная математика в примерах и задачах. М.: Наука, 1972. 368 с.

Приложение 1. Листинг программы

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

namespace Числаки

{

class RotationMethod

{

public static double Max(double[,] matrix, out int index_i, out int index_j)

{

index_i = 0;

index_j = 1;

double max = Math.Abs(matrix[1,0]);

for (int i = 0; i < matrix.GetLength(0); i++)

for (int j = i + 1; j < matrix.GetLength(1); j++)

if(Math.Abs(matrix[i, j]) > max)

{

max = Math.Abs(matrix[i, j]);

index_i = i;

index_j = j;

}

return max;

}

public static double[,] CreateRotationMatrix(int rows, int columns, int i, int j, double cos_phi, double sin_phi)

{

double[,] matrix_h = new double[rows, columns];

matrix_h[i,i] = cos_phi;

matrix_h[j,j] = cos_phi;

matrix_h[i,j] = - sin_phi;

matrix_h[j,i] = sin_phi;

for (int k = 0; k < rows; k++)

for (int l = 0; l < columns; l++)

{

if (k == l && k != i && k != j)

{

matrix_h[k, l] = 1;

}

}

return matrix_h;

}

public static double[,] Transpose(double[,] matrix)

{

double[,] matrix_transpose = new double[matrix.GetLength(1), matrix.GetLength(0)];

for (int i = 0; i < matrix.GetLength(0); i++)

for (int j = 0; j < matrix.GetLength(1); j++)

matrix_transpose[j, i] = matrix[i, j];

return matrix_transpose;

}

public static double[,] Multiply(double[,] matrix_1, double[,] matrix_2)

{

if(matrix_1.GetLength(0) != matrix_2.GetLength(1)) throw new Exception("Некорректная размерность матриц");

double[,] matrix_multiply = new double[matrix_1.GetLength(0), matrix_2.GetLength(1)];

for (int i = 0; i < matrix_multiply.GetLength(0); i++)

for (int j = 0; j < matrix_multiply.GetLength(1); j++)

for (int k = 0; k < matrix_1.GetLength(0); k++)

{

matrix_multiply[i, j] += matrix_1[i, k] * matrix_2[k, j];

}

return matrix_multiply;

}

public static double[] Rotation(double[,] matrix, double eps, out double[,] matrix_vector)

{

double[] vector_eigenvals = new double[matrix.GetLength(0)];

double[,] matrix_a = matrix;

int i, j;

int count = 0;

double max = Max(matrix_a, out i, out j);

Console.WriteLine("Max={0:0.###},\ti={1:0.###},\tj={2:0.###}", max, i, j);

double cos_phi = Math.Sqrt(0.5 * (1 + 1 / Math.Sqrt(1 + (2 * matrix_a[i, j] / (matrix[i, i] - matrix[j, j])) * (2 * matrix_a[i, j] / (matrix[i, i] - matrix[j, j])))));

double sin_phi = Math.Sqrt(0.5 * (1 - 1 / Math.Sqrt(1 + (2 * matrix_a[i, j] / (matrix[i, i] - matrix[j, j])) * (2 * matrix_a[i, j] / (matrix[i, i] - matrix[j, j])))));

Console.WriteLine("Cos_phi={0}\tSin_phi={1}", cos_phi, sin_phi);

string title_h = "Матрица H_" + count + ":";

double[,] matrix_h = CreateRotationMatrix(matrix_a.GetLength(0), matrix_a.GetLength(1), i, j, cos_phi, sin_phi);

PrintMatrix(title_h, matrix_h);

string title_a = "Матрица A_" + count + ":";

matrix_a = Multiply(Multiply(Transpose(matrix_h), matrix_a), matrix_h);

PrintMatrix(title_a, matrix_a);

matrix_vector = matrix_h;

while (max > eps)

{

count++;

max = Max(matrix_a, out i, out j);

Console.WriteLine("Max={0:0.###},\ti={1:0.###},\tj={2:0.###}", max, i, j);

cos_phi = Math.Sqrt(0.5 * (1 + 1 / Math.Sqrt(1 + (2 * matrix_a[i, j] / (matrix_a[i, i] - matrix_a[j, j])) * (2 * matrix_a[i, j] / (matrix_a[i, i] - matrix_a[j, j])))));

sin_phi = Math.Sqrt(0.5 * (1 - 1 / Math.Sqrt(1 + (2 * matrix_a[i, j] / (matrix_a[i, i] - matrix_a[j, j])) * (2 * matrix_a[i, j] / (matrix_a[i, i] - matrix_a[j, j])))));

Console.WriteLine("Cos_phi={0}\tSin_phi={1}", cos_phi, sin_phi);

title_h = "Матрица H_" + count + ":";

matrix_h = CreateRotationMatrix(matrix_a.GetLength(0), matrix_a.GetLength(1), i, j, cos_phi, sin_phi);

PrintMatrix(title_h, matrix_h);

title_a = "Матрица A_" + count + ":";

matrix_a = Multiply(Multiply(Transpose(matrix_h), matrix_a), matrix_h);

PrintMatrix(title_a, matrix_a);

matrix_vector = Multiply(matrix_vector, matrix_h);

}

for (int l = 0; l < vector_eigenvals.GetLength(0); l++)

vector_eigenvals[l] = matrix_a[l, l];

return vector_eigenvals;

}

static void PrintMatrix(string header, Array matrix)

{

int i = 0;

Console.WriteLine();

Console.WriteLine(header);

foreach (object x in matrix)

{

if (i == matrix.GetLength(1))

{

Console.WriteLine();

i = 0;

}

i++;

Console.Write("{0,11:0.#######}", x);

}

Console.WriteLine();

Console.WriteLine();

}

static void PrintVector(string header, string index, Array vector)

{

int i = 1;

Console.WriteLine();

Console.WriteLine(header);

foreach (object x in vector)

{

Console.WriteLine("{0,6}[{1}]={2}", index, i, x);

i++;

}

}

static void Main()

{

try

{

double[,] matrix = { { 5, 1, 2 }, { 1, 4, 1 }, { 2, 1, 3 } };

double eps = 1e-3;

double[,] matrix_vector = new double[matrix.GetLength(0), matrix.GetLength(1)];

double[] vector_eigenvals = new double[matrix.GetLength(0)];

PrintMatrix("Исходная матрица:", matrix);

vector_eigenvals = Rotation(matrix, eps, out matrix_vector);

PrintVector("Собственные значения", "lambda", vector_eigenvals);

PrintMatrix("Собственные вектора:", matrix_vector);

Console.ReadLine();

}

catch(Exception e)

{

Console.WriteLine(e.Message);

Console.ReadLine();

}

}

}

}

Приложение 2. Пример результата работы программы



Размещено на Allbest.ru