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

Всякая симметричная матрица Aможет быть за­­писана в виде

A = V^*DV, (1.43)

где V^*- ортогональная матрица;D- ди­а­го­наль­ная с эле­мен­тамиl₁, l₂,..., l_m- собственными зна­че­ни­я­­ми мат­ри­цыA. После умножения (1.43) слева на V^*получимAV^* = V^*D, и ес­­ли рас­писать этоматричное ра­венство по стол­б­цам, то ока­жет­ся, что каждыйi-й столбец мат­ри­цыV^*является соб­ст­венным век­тором, от­ве­ча­ю­щим соб­ст­вен­ному зна­че­ниюl_i[Ба­х­валов, 1973а]. Пре­­об­­разуя сис­те­му (1.43) к ви­­ду

VAV^*= D,

по­стро­им по­сле­до­ва­­тель­ность ор­то­го­наль­ных мат­риц V_1, ...., V_n так, чтобы приW_n=V_n...V₁,иметь W_n AW_n^* D.Если теперь обо­значитьче­­­резV_kl(f)мат­ри­цу с эле­мен­­та­ми v_kk = v_ll = cosf, v_kl = -v_lk = sin f, v_ii = 1 приi№ k, i№1, v_ij = 0 при ос­таль­ных (i, j),то она будет ортогональной.

Построим последовательность матриц V_nиa⁽ⁿ⁾ =по следующему правилу:

А⁽¹⁾ = V₁ A V₁^* , ..., A⁽ⁿ⁾ = V _nA^(n-1) V_n^* . (1.44)

Матрица

строится в зависимости от матрицы A⁽ⁿ⁾так, что­бы вы­пол­ня­лось неравенство [Бахвалов, 1973а]

. (1.45)

Выражая каждую матрицу A⁽ⁿ⁾через пре­ды­ду­щую, бу­демстроить итерационную по­сле­до­ва­тель­ностьA⁽ⁿ⁾ =WAW_n^*. Согласно оцен­ке (1.45), имеем

.

Пусть матрица D⁽ⁿ⁾получена из A⁽ⁿ⁾простой за­меной всех недиагональных элементов на нули. Тогда не­труд­но доказать [Бахвалов, 1973а], что

, (1.46)

откуда вытекает |l_n⁽ⁿ⁾ - l_n| £ r_n при всехn, где l_i - соб­ст­вен­ные значения матрицыA⁽ⁿ⁾(значит, и мат­рицыA), за­­ну­ме­ро­ванные в порядке убы­ва­ния, аl_n⁽ⁿ⁾ - соб­ствен­ные зна­че­ния матрицы D⁽ⁿ⁾,т.е. ее диагональные эле­мен­ты, также за­ну­­ме­ро­ван­ные в порядке убывания;r_n- верх­няя оцен­ка по­греш­ности.

Столбцы матрицы W⁽ⁿ⁾оказываются при­бли­же­ни­я­ми к собственным векторам матрицыA⁽ⁿ⁾. Все столбцы про­изведения матриц V(f_kl) иA_kl, кро­меk-го и первого, сов­па­дают с соответствующими столб­цами матрицыA⁽ⁿ⁾,k-й и первый столбцы но­вой матрицы яв­ля­ют­ся линейными ком­­би­на­ци­я­ми этих жестолбцов мат­рицы A. Та­когоже ко­л­и­чес­т­­ва дейс­твий по­тре­бу­ет дальнейшее умно­жение на мат­­рицуV^*(f_kl) в со­от­вет­ствии с фор­му­ла­ми (1.44), ког­да про­ис­хо­дит из­ме­не­ниеk-й и пер­вой строк.

Общее количество арифметических операций для по­­лу­­че­ния каждой последующей матрицы A⁽ⁿ⁾со­с­та­вит при этом ~ 6n.

Рассмотренный алгоритм (с некоторой мо­ди­фи­ка­ци­ей, предложенной в работе [Greenstadt, 1961]) ре­­а­ли­зо­ван в про­­це­дуре eigjac,используя которую можно вы­чис­лить все собственные зна­че­ния и соб­ствен­ные век­то­ры за­дан­­ной сим­мет­рич­ной матрицыa[1:n, 1:n]. На вы­хо­­де из про­це­­ду­ры l-й столбец мат­рицыS(на­чаль­ное зна­че­ние мас­си­ваS не­су­щест­венно) со­дер­жит первый изn соб­­­ст­­вен­ных век­­­то­ров дан­ной мат­рицы, ди­а­го­наль­ный эле­ментa[l, l] мас­­сиваAсо­­от­вет­ст­вует первому соб­­ст­вен­ному зна­­­че­нию. Ите­­ра­ци­онный процесс за­­ка­н­чи­ва­ет­ся, ког­да для каж­дого не­­­ди­а­го­наль­но­­го элементаa[i,j] вы­пол­ня­­ется не­ра­вен­ст­во

,

где параметр e,связанный с оценкой по­греш­нос­тиr_nв фор­му­ле (1.46) и введенный вместоr_nдля оценки ка­чес­т­ва ре­ше­­ния в работе [Gre­en­stadt, 1961], и есть за­да­ва­емая точ­ность.

Формальные параметры про­цедуры. Входные:n(тип in­teger) - порядок матрицыA;eps(типreal) - за­да­­ваемая точ­­­ность;a[1:n,1:n] (типreal) - ис­ход­ная сим­­метричная мат­­­рица.Вы­ход­ные:a[1:n, 1:n] (типre­al) - матрица, ди­а­го­наль­­ные элементы ко­то­­ройa[i,i] яв­ляются со­от­вет­ст­вен­ноl_i-ми соб­ст­вен­ны­ми зна­че­­ни­­я­ми;s[1:n, 1:n] (типre­al) - мат­­ри­ца,k-й столбец ко­­то­рой содержитk-й изnсоб­ст­вен­ных векторов ис­ход­­ной матрицы A.

procedure eigjac (n:integer;eps:real;

var a, s:mas1);

var norm1,norm2,thr,mu,omega,sint,cost,

int1,v1, v2,v3,t : real;

i,j,p,q,ind : integer;

label main, main1;

begin

{ *** формирование в массиве s

единичной матрицы *** }

for i:=1 to n do

for j:=1 to i do

if i=j then

s[i,j]:=1

else

begin

s[i,j]:=0;

s[j,i]:=0;

end;

{ *** вычисление начальной нормы norm1,

конечной нормы norm2 и порога thr *** };

int1:=0;

for i:=2 to n do

for j:=1 to i-1 do

int1:=int1+2*a[i,j]*a[i,j];

if int1=0 then

exit;

norm1:=sqrt(int1);

thr:=sqrt(int1);

norm2:=(eps/n)*norm1;

ind:=0;

main:

thr:=thr/n;

{**начало обработки недиагональных элементов**};

main1:

for q:=2 to n do

for p:=1 to q-1 do

if abs(a[p,q]) >= thr then

begin

ind:=1;

v1:=a[p,p];

v2:=a[p,q];

v3:=a[q,q];

mu:=0.5*(v1-v3);

if mu=0 then

omega:=-1

else

omega:=-sign(mu)*

v2/sqrt(v2*v2+mu*mu);

sint:=omega/sqrt(2*(1+

sqrt(1-omega*omega));

cost:=sqrt(1-sint*sint);

for i:=1 to n do

begin

if i<>p AND i<>q then

begin

int1:=a[i,p];

mu:=a[i,q];

a[q,i]:=int1*sint+mu*cost;

a[i,q]:=int1*sint+mu*cost;

a[p,i]:=int1*cost-mu*sint;

a[i,p]:=int1*cost-mu*sint;

end;

int1:=s[i,p];

mu:=s[i,q];

s[i,q]:=int1*sint+mu*cost;

s[i,p]:=int1*cost-mu*sint

end {***i***};

mu:=sint*sint;

omega:=cost*cost;

int1:=sint*cost;

t:=2*v2*int1;

a[p,p]:=v1*omega+v3*mu-t;

a[p,p]:=v1*mu+v3*omega+t;

a[p,q]:=(v1-v3)*int1+v2*(omega-mu);

a[q,p]:=(v1-v3)*int1+v2*(omega-mu);

end ;

{***проверка достижения заданной

точности***};

if ind=1 then

begin

ind:=0;

go to main1

end;

if thr > norm2 then go to main;

end.

Процедура eigjac получена из процедуры jacobi [Агеев, 1976] переводом последней с языка AL­G­OL на язык PAS­CAL, при этом было сделано не­сколь­ко тож­дественных пре­образований, оп­ти­ми­зи­ру­ю­щих зат­раты машинного вре­мени. Про­­це­­­ду­раeig­jac бы­ла про­­ве­ре­на для тех же ис­­ход­ных данных, что и про­цедураjacobi [Агеев, 1976]:n = 4,eps=10^-8и

.

Вычисления, выполненные на IBM PC/AT-386 с час­то­той 40 МГц, дали следующие результаты, которые сов­падают с ре­зуль­та­та­ми тес­ти­рования

_;

,

при­ве­денными в ра­бо­­те Агеева [Аге­ев и др., 1976], и име­ют мак­си­маль­ную по­греш­ность по от­но­ше­нию к соб­ст­вен­ным зна­че­ни­ям l_i и собственным век­­­то­рамXис­­ход­ной мат­рицы [в сравнении с их ве­­­ли­чи­на­ми из ра­бо­­ты Фадеева и др. (1963) ] max r_i= 4*10^-8, maxS_i = 3*10^-6.