Справочник функций

Ваш аккаунт

Войти через: 
Забыли пароль?
Регистрация
Информацию о новых материалах можно получать и без регистрации:

Почтовая рассылка

Подписчиков: -1
Последний выпуск: 19.06.2015

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

94K
29 августа 2014 года
kkkk5300
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.
Прикрепленные файлы:
50 Кб
Загрузок: 998
21 Кб
Загрузок: 999
39 Кб
Загрузок: 957
41 Кб
Загрузок: 973
57 Кб
Загрузок: 959

Знаете кого-то, кто может ответить? Поделитесь с ним ссылкой.

Ваш ответ

Реклама на сайте | Обмен ссылками | Ссылки | Экспорт (RSS) | Контакты
Добавить статью | Добавить исходник | Добавить хостинг-провайдера | Добавить сайт в каталог