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

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

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

Содержание

Введение

1. Основные определения матричного исчисления

2. Основные свойства собственных значений

3. Преобразование подобия матриц

4. Матрица вращения

5. Метод Гивенса

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

7. Характеристический многочлен матрицы

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

Заключение

Список используемых источников

Приложение



Введение

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

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



1. Основные определения матричного исчисления

Матрица A называется симметричной, если а_ij = а_ij, где i, j = 1, 2, . . ., n.

Отсюда следует симметрия относительно диагонали а_kk, где k = 1, 2, . n.

Матрица

1	4	5
4	3	7
5	7	2

является примером симметричной.

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

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

Матрица A называется ортогональной, если

А^ТА = Е

где А^Т - транспонированная матрица A, а Е - единичная матрица. Очевидно, матрица, обратная ортогональной, эквивалентна транспонированной.

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

В = Р-1АР

2. Основные свойства собственных значений

Если А - квадратная матрица размерности n и матричное равенство

А • Х_i = л_i • X_i

Выполнено при некотором Х_i ? 0, тогда л_i называют собственным значением (собственным числом) матрицы А, а ненулевой вектор Х_i , соответствующий л_i называют собственным вектором матрицы А:

Х_i =(x_i1 , x_i2 ,…, x_in)^T

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

3. Преобразование подобия матриц

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

Преобразованием подобия называют преобразование вида:

С = Т^-1АТ, | Т | ? 0

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

4. Матрица вращения

Элементарной матрицей вращения Т_ij называется единичная матрица Е, в которой изменены четыре элемента: вместо единиц в позициях и поставлено число с, в позиции поставлено число s, а в позиции поставлено -s:

При этом числа с и s удовлетворяют соотношению: .

Свойства матрицы вращения:

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

1. Матрица ортогональная

2. При умножении матрицы А слева на матрицу вращения •А изменяются только i и j строки матрицы А. Если , то

(1)

При умножении справа (при согласованных размерностях) изменяются только i и j столбцы матрицы А. Если , то

(2)

Преобразование вида (1) над расширенной матрицей системы линейных алгебраических уравнений приводит к расширенной матрице эквивалентной системы.

3. Если:

То в матрице •А в позиции будет стоять ноль.

4. Суммы квадратов элементов соответствующих столбцов матриц •А и А равны.

Действительно, из (2) получаем

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

5. Метод Гивенса

Метод Гивенса относится к прямым методам. Метод Гивенса позволяет за N шагов привести матрицу А к трехдиагональному виду, если А - симметричная, либо к левому почти треугольному виду, если А - произвольная квадратная матрица. Этот метод основан на подобных преобразованиях исходной матрицы с помощью матриц вращения:



где i - 2,3,…n - 1; j = i+1, i+2,…n; k = 1,2,…N,

A0 = A; Tij - матрица вращения.

В итоге получаем:

если А была симметричной, либо

если А - была произвольной матрицей.

В результате преобразования вида обнуляется элемент матрицы Аk-1, т.е. . Для получения матрицы АN потребуется

таких преобразований.

Пусть T = T23 • T24 • … • T2n • … • Tn-1,n , тогда и можно записать: AN = TT • A • T , где Т - ортогональная матрица. Это означает, что матрица AN подобна А и имеет те же собственные значения. Следовательно, для исходной матрицы будет также решена проблема собственных значений.

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

Рассмотрим квадратную матрицу A = (ai,j) порядка n. Пусть Tij матрица вращения порядка n с элементами c и s такими, что c2 + s2 = 1.

Построим последовательность из подобных матриц, с начальной матрицей A, такую, что у последней матрицы в позициях (1,3), (1,4), ...(1, n), (2,4), (2,5), ..., (2, n), ..., (n ? 2, n) стояли бы нули. Другими словами, все элементы последней матрицы в позициях (i, j) таких, что i + 1 < j были бы нули, а последняя матрица была трехдиагональной (матрица A симметричная). Число матриц последовательности равно числу обнуляемых элементов. Последнюю матрицу такой последовательности условно можно записать так:

(1)

В этой записи присутствуют все матрицы набора.

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

(2)

при i = 2, j = 3. Пусть B = A•Ti,j , тогда . Все столбцы матрицы A совпадают со столбцами матрицы B = (bij) за исключением i, j столбцов.

Аналогично, все строки G = (gij) совпадают со строками B, за исключением i, j строк. Элементы матриц G и В в позиции (i - 1, j) совпадают gi?1,j = bi?1,j . Потребуем, чтобы подобное преобразование (2) обнулило элемент gi?1,j матрицы G в позиции (i ? 1, j). Так как элемент gi?1,j лежит только на одном изменяющимся j-ом столбце, то, в силу совпадения элементов gi?1,j = bi?1,j , получим Требуем, чтобы

Отсюда, с учетом с2 + s2 = 1, находим

(3)

Подставим найденные элементы с и s в матрицу вращения в (2) и получим вторую матрицу последовательности G.

Аналогично проводим следующее умножение в (1) и получаем третью матрицу

(4)

при i = 2, j = 4. Отсюда и из (2), подставляя соответствующие значения индексов, получим часть формулы (1)

Важно отметить, что элементы с и s для матрицы вращения в (4) вычисляются по формулам (3), но только с помощью элементов второй матрицы G.

Продолжая подобным образом, получим последнюю матрицу S, обладающую нужными свойствами.

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

В случае матрицы размерности n*n метод Гивенса требует n-2 основных шагов, на каждом из которых выполняется ряд преобразований, число которых зависит от числа нулей, которое хотят получить в данном столбце или строке. На k-м шаге обращают в нули элементы, стоящие вне трех диагоналей k-й строки и k-го столбца, сохраняя в то же время нулевые элементы, полученные на предыдущих шагах. Таким образом, перед началом k-го шага преобразованная матрица является трехдиагональной, если ограничиться рассмотрением ее первых k-1 строк и столбцов. По мере преобразований симметричная матрица размерности 5х5 приобретает следующие формы:

исходная матрица

после первого основного шага, состоящего из трех преобразований:

матричный исчисление гивенс бисекция



после второго основного шага, состоящего из двух преобразований:

после третьего основного шага, состоящего и одного преобразования. Теперь матрица имеет трехдиагональный вид:

На каждом основном шаге изменяются лишь те элементы матрицы, которые расположены в ее правой нижней части. Таким образом на k-м шаге преобразуется только матрица порядка (n - k + 1), занимающая правый нижний гол исходной матрицы. Ясно, что на каждой следующей стадии выполняется меньшее число преобразований, чем на предыдущей. Всего для приведения матрицы к трехдиагональному виду требуется выполнить (n2 - 3n + 2)/2 преобразований.



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

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

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

7. Характеристический многочлен матрицы

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

(1)

Минор k-го порядка , составленный из элементов матрицы, стоящих на пересечении одноименных строк и столбцов, называется главным минором. В формуле (1) коэффициент при равен сумме главных миноров k-го порядка, в частности, след матрицы - это сумма главных миноров 1-го порядка, определитель матрицы - это главный минор n-го порядка.

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

Найдем собственные числа симметричной трехдиагональной матрицы

Пусть Sk матрица k-го главного минора матрицы Sn, Pk(л) - характеристический многочлен матрицы Sk :

Разложим определитель Pk(л) по последней строке, получим следующее рекуррентное соотношение:

Добавим к этому списку еще два многочлена P0(л) = 1 и P1(л) = б1 - л.

Если вi ? 0, i = 2, 3, … , n, то корни полинома Pk-1(л) строго разделяют корни полинома Pk(л) для любого k = :

Все собственные числа матриц Sk, k ? n, находятся в пределах от a = ?||Sn|| и до b = ||Sn|| так как из соотношений Sk•u = лu, ||u|| = 1 и очевидного неравенства ||Sk|| ? ||Sn|| вытекает

|л| = |л|||u|| = ||л u|| = ||Sk•u|| ? ||Sk||||u|| ? ||Sk|| < ||Sn||

Поэтому ?||Sn|| ? л ? ||Sn||.

Свойство разделения корней многочленов Pk(л) позволяет определить интервалы изоляции корней каждого из многочленов. Так, многочлен P1(л) = б1?л обладает одним корнем л = б1 на отрезке [a,b]. Так как этот корень разделяет корни л1 и л2 многочлена P2(л): л1 < б1 < л2, то интервалы изоляции корней многочлена P2(л) будут

отрезки [a, б1] и [б1, b]. Найдя корни многочлена P2(л) находим интервалы изоляции корней многочлена P3(л) и так далее. Зная интервалы изоляции корней, можно найти корни всех многочленов, например, методом деления отрезка пополам (метод бисекций).

Если вm = 0 при некотором m ? n, то определитель Pn(л) можно представить как произведение двух определителей и для каждого из них применить метод бисекции.



Заключение

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

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

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

Сам код приведен в Приложении



Список используемых источников

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

2. Мостовской А.П. Численные методы и система Mathematica / А.П.Мостовской - Мурманск: 2009. - 249с.

3. Молчанов И.Н. Машинные методы решения прикладных задач. Алгебра, приближение функций - Киев: Наукова Думка, 1987. - 288с.

4. Вержбицкий В.М. Численные методы. Линейная алгебра и нелинейные уравнения - Москва: Высшая школа, 2000. - 263с.

5. Ивтб.рф/exams/вычмат/24.htm

6. http://www.matcabi.net/matrix_ch.php

7. http://www.mathforyou.net/AnyPowSol.html

8. http://matrix.reshish.ru/multiplication.php



Приложение

uses crt;

const w=3;

var A,G,T,T_t,B:array[1..w,1..w]of real;

acb:array[1..w+2] of real;

xxx:array[1..w+2] of real;

i,j,z,x,y,m,k:integer; a1,b1,c,s,pp,e,sym,ev_norm:real;

procedure ym_m(a,b:array[1..w,1..w]of real; var r:array[1..w,1..w]of real);

var i,j,z:integer; k:real;

begin

k:=0;

for i:=1 to w do

for j:=1 to w do begin

for z:=1 to w do // умножение матриц

k:=a[i,z]*b[z,j]+k;

r[i,j]:=k;

k:=0;

end;

end;

function det(x:real):real;

begin

if x=11 then det:=A[2,2]*A[3,3]-A[2,3]*A[3,2];

if x=22 then det:=A[1,1]*A[3,3]-A[1,3]*A[3,1]; //определители матрицы 2х2

if x=33 then det:=A[1,1]*A[2,2]-A[1,2]*A[2,1];

end;

function f(y:real):real;

begin

if k=1 then f:=A[1,1]-y;

if k=2 then f:=y*y+(A[1,1]+A[2,2])*(-y)+det(33); //характер многочл

if k=3 then f:=(-1)*y*y*y+(A[1,1]+A[2,2]+A[3,3])*(y*y)+

+(det(11)+det(22)+det(33))*(-y)+(A[1,1]*A[2,2]*A[3,3]+A[1,3]*A[2,1]*A[3,2]+

+A[1,2]*A[2,3]*A[3,1]-A[1,3]*A[2,2]*A[3,1]-A[1,1]*A[2,3]*A[3,2]-

-A[1,2]*A[2,1]*A[3,3]);

end;

begin

clrscr;

writeln('Введите элементы матрицы 3х3');

for i:=1 to w do

for j:=1 to w do //заполняю матрицу A

if i<=j then readln(a[i,j]) else a[i,j]:=a[j,i];

for i:=1 to w do

for j:=1 to w do

G[i,j]:=A[i,j]; //передаю массив G=A

for x:=1 to w do begin

for y:=1 to w do

T[x,y]:=0;

T[x,x]:=1; // делаю Т единичной

end;

i:=2; j:=3;

c:=G[i-1,i]/sqrt(sqr(G[i-1,i])+sqr(G[i-1,j]));

s:=-G[i-1,j]/sqrt(sqr(G[i-1,i])+sqr(G[i-1,j]));

Т[i,i]:=c; Т[j,j]:=c; Т[i,j]:=s; Т[j,i]:=-s;

for x:=1 to w do

for y:=1 to w do // транспонирование T_t

T_t[x,y]:=T[y,x];

ym_m(T_t,G,B);

ym_m(B,T,G);

writeln('Трехдиагональная симметричная:') ;

for i:=1 to w do begin //вывод массива

for j:=1 to w do

write(g[i,j]:8:3);

writeln;

end;

// нахождение корней

sym:=0;

for i:=1 to w do

for j:=1 to w do sym:=sqr(A[i,j])+sym; // эвклидова норма матрицы

ev_norm:=sqrt(sym);

acb[1]:=-ev_norm;

acb[2]:=ev_norm;

e:=0.0001;

k:=1;

while k<4 do begin

for i:=1 to k do begin

a1:=acb[i]; b1:=acb[i+1];

while abs(b1-a1)>e do begin

c:=(a1+b1)/2;

if f(a1)*f(c)<0 then b1:=c

else

if f(c)=0 then break else a1:=c;

end;

xxx[i+2]:=ev_norm;

xxx[i+1]:=c;

end;

for j:=2 to i+2 do

acb[j]:=xxx[j];

k:=k+1;

end;

writeln;

writeln('Корни уравнения: ',acb[2]:7:3, acb[3]:7:3, acb[4]:7:3);

end.

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