Метод вращений для нахождения собственных чисел и векторов

94K

29 августа 2014 года

1 / / 29.08.2014

Здравствуйте!
У меня возникла проблема с реализацией метода вращения для определения собственных значений и векторов пары матриц.
Для матриц размерность N=2 считает вроде бы верно, для N=3 иногда выдает верные ответы, а для N>3 выдает совсем неверно.
В приложении представлен алгоритм, нужно реализовать именно его.

Прошу помочь, найти ошибку. Может я что-то не понял, пояснить.
Вот словесное пошаговое описание алгоритма (так, как я его понял):

Замечания:
1) K и M - симметричные матрицы
2) Для транспонирования матрицы P можно просто поменять местами alpha и gamma
3) массив L - это массив собственных чисел (лямбды)
4) массив PL - это массив собственных чисел, вычисленных на предыдущем шаге
5) матрица X - массив собственных векторов (в алгоритме тоже обозначена как X)
6) цикл с постусловием в алгоритме обозначен итерацией t
7) циклы i и j вместе в алгоритме обозначены итерацией s
8) циклы i, j и цикл с постусловием все вместе обозначены итерацией r

Инициализируем элементы:
a) P, X - единичные матрицы
b) L = K / M, i=1..n (по формуле 6.5)

1. Начинаем цикл с постусловием (проверка по формуле 6.6)
1.1. Начинаем цикл i=2,..,n
1.1.01. Начинаем цикл j=1,..,i-1
1.1.02. Находим ai, aj, a, z, gamma, alpha по формулам 6.4
1.1.03. Присваиваем P[j] = alpha, P[j] = gamma (строим сразу транспонированную матрицу)
1.1.04. K = P * K (матричное умножение; P мы построили так, что это транспонированная матрица; часть формулы 6.2)
1.1.05. M = P * M (часть формулы 6.2)
1.1.06. Присваиваем P = gamma, P[j] = alpha (строим нормальную матрицу, не транспонированную)
1.1.07. K = K * P (оставшаяся часть формулы 6.2)
1.1.08. M = M * P (оставшаяся часть формулы 6.2)
1.1.09. X = X * P (по формуле 6.7)
1.1.10. Зануляем P[j] и P[j] (теперь у нас матрица P снова единичная)
1.1.11. Присваиваем PL = L
1.1.12. Обновляем значения массива L по формуле 6.5
1.1.13. Выходим из цикла j
1.2. Выходим из цикла i
2. Проверяем постусловие цикла (начало цикла в п.1): если для всех i=1..n выполняется abs(L - PL) < eps, то завершаем цикл.

Теперь ответ будет лежать в L и X.

Вот, собственно, код, который я написал:

Свернуть исходник

Код:

program task4;

{$APPTYPE CONSOLE}

uses
SysUtils,
Windows;

const
SZ = 10;
eps = 1e-5;

type
arrM = array[0..SZ] of real;
matrM = array[0..SZ] of array[0..SZ] of real;

procedure my_writeln(const s: string);
var NewStr: string;
begin
SetLength(NewStr, Length(s));
CharToOem(PChar(s), PChar(NewStr));
Writeln(NewStr);
end;

procedure multiplication(n: byte; a, b: matrM; var c: matrM); overload;
var i, j, k: byte;
begin
for i := 1 to n do begin
for j := 1 to n do begin
c[i][j] := 0;
for k := 1 to n do
c[i][j] := c[i][j] + a[i][k] * b[k][j];
end;
end;
end;

function Sign(a: real): integer;
begin
if (a >= 0) then
Result := 1
else
Result := -1;
end;

function method(n: byte; K, M: matrM; var L: arrM; var X: matrM): boolean;
var
i, j, e: byte;
ai, aj, a, z, gamma, alpha: real;
P: matrM;
PL: arrM;
check: boolean;

begin
for i := 1 to n do begin
L[i] := K[i][i] / M[i][i];
P[i][i] := 1;
end;

for i := 2 to n do
for j := 1 to i - 1 do begin
P[i][j] := 0;
P[j][i] := 0;
end;

X := P;
repeat
for i := 2 to n do begin
for j := 1 to i - 1 do begin
ai := K[i][i] * M[i][j] - M[i][i] * K[i][j];
aj := K[j][j] * M[i][j] - M[j][j] * K[i][j];
a := K[i][i] * M[j][j] - M[i][i] * K[j][j];
z := a / 2 + sign(a) * Sqrt((a / 2) * (a / 2) + ai * aj);
gamma := - ai / z;
alpha := aj / z;

P[i][j] := alpha;
P[j][i] := gamma;
multiplication(n, P, K, K);
multiplication(n, P, M, M);

P[i][j] := gamma;
P[j][i] := alpha;
multiplication(n, K, P, K);
multiplication(n, M, P, M);
multiplication(n, X, P, X);

P[i][j] := 0;
P[j][i] := 0;

PL := L;
for e := 1 to n do
L[e] := K[e][e] / M[e][e];
end;
end;

check := true;
for i := 1 to n do
if (abs(L[i] - PL[i]) >= eps) then
check := false;
until (check);
end;

var
n, i, j: byte;
K, M, X: matrM;
L: arrM;

begin
my_writeln('Введите размерность матрицы N: ');
readln(n);

my_writeln('Введите матрицу K[N][N]:');
for i := 1 to n do
for j := 1 to n do
read(K[i][j]);

my_writeln('Введите матрицу M[N][N]:');
for i := 1 to n do
for j := 1 to n do
read(M[i][j]);

method(n, K, M, L, X);

writeln;
my_writeln('Собственные числа: ');
for i := 1 to n do
write(Format('%0.2f ', [L[i]]));
writeln;

writeln;
for i := 1 to n do begin
for j := 1 to n do
write(Format('%0.2f ', [X[i][j]]));
writeln;
end;

sleep(20000);
end.