Метод вращений для пары симметричных матриц

@kkkk5300 · Регистрация: 04.07.2014

Студворк — интернет-сервис помощи студентам

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

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

Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
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.

Приложение

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

Задание

Дана пара симметричных матриц K и M размерности n. Требуется:
1. Составить алгоритм и программу для определения всех собственных значений и собственных векторов пары матриц K и M с использованием метода вращений (обобщённого метода Якоби).
2. Использовать программу для решения варианта задания.
3. Выполнить проверку и сделать выводы.

Варианты заданий

$https://www.cyberforum.ru/cgi-bin/latex.cgi?K=\begin{bmatrix}50+3i &10-i &3 &4 &-30-i \\\\\\10-i &20+2i &10-i &3 &4 \\\\\\3 &10-i &90-i &10-i &3 \\\\\\4 &3 &10-i &20+2i &10-i \\\\\\-30-i &4 &3 &10-i &50+3i\end{bmatrix}$

$https://www.cyberforum.ru/cgi-bin/latex.cgi?M=\begin{bmatrix}60-i &5+i &3 &1 &10 \\\\5+i &10+2i &5-i &3 &1 \\\\3 &5-i &30+i &5-i &3 \\\\1 &3 &5-i &10+2i &5+i \\\\10 &1 &3 &5+i &60-i\end{bmatrix}$

где i - номер варианта.

Методические указания

Собственные векторы x_i и собственные значения λ_i пара матриц K и M удовлетворяют условию

$https://www.cyberforum.ru/cgi-bin/latex.cgi?(K-\lambda M)x=0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (6.1)$

Задача определения λ_i и x_i, i=1,...,n, называется обобщённой проблемой собственных значений. В методе вращений решения этой задачи матрицы K и M одновременно приводят к диагональному виду с помощью матриц вращения по формулам

$https://www.cyberforum.ru/cgi-bin/latex.cgi?\begin{matrix}K_{s+1}=P_s^TK_sP_s\ ;\\ M_{s+1}=P_s^TM_sP_s\ ,\end{matrix}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (6.2)$

где K₁=K; M₁=M, а несимметричная матрица вращения на s-м шаге P_s имеет вид

$https://www.cyberforum.ru/cgi-bin/latex.cgi?P_s=\begin{bmatrix}1 & & & & \\\\ &1 & &\alpha & \\\\ & &1 & & \\\\ &\gamma & &1 & \\\\ & & & &1 \end{bmatrix}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (6.3)$

У матрицы P_s все диагональные элементы равны единице, а из внедиагональных - отличны от нуля только два элемента: P_ij= γ и P_ji=α. Значения γ и α выбирают так, чтобы после преобразования (6.2) элементы k_ij^(s+1) и m_ij^(s+1) матриц K_s+1 и M_s+1 равнялись нулю. Вычисление элементов γ и α производится с помощью следующего алгоритма:

$https://www.cyberforum.ru/cgi-bin/latex.cgi?\begin{matrix}a_i=k_{ii}^{(s)}m_{ij}^{(s)}-m_{ii}^{(s)}k_{ij}^{(s)}\ ;\\\\a_j=k_{jj}^{(s)}m_{ij}^{(s)}-m_{jj}^{(s)}k_{ij}^{(s)}\ ;\\\\a=k_{ii}^{(s)}m_{jj}^{(s)}-m_{ii}^{(s)}k_{jj}^{(s)}\ ;\\\\\zeta =\frac{a}{2}+\textrm{sign}(a)\sqrt{\left( \frac{a}{2}\right)^2+a_ia_j}\ ;\\\\\gamma =-\,\frac{a_i}{\zeta };\ \ \ \alpha =\frac{a_j}{\zeta }\ .\end{matrix}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (6.4)$

Номера i, j выбирают так, чтобы аннулировать наибольшие по модулю внедиагональные элементы k_ij^(s) и m_ij^(s). В простейшем случае на итерации с номером t допустимо организовать двойной цикл по i=2,...,n и j=1,...,i-1 с тем, чтобы аннулировать поочерёдно внедиагональные элементы k_ij^(s) и m_ij^(s). На новой, (t+1)-й итерации цикл повторяют.

С увеличением количества вращений по формулам (6.2)-(6.4) матрицы K_s+1 и M_s+1 приобретают диагональный вид. Тогда приближение к i-у собственному значению после s-го вращения определяют по формуле

$https://www.cyberforum.ru/cgi-bin/latex.cgi?\lambda _i^{(s+1)}=\frac{k_{ii}^{(s+1)}}{m_{ii}^{(s+1)}},\ \ i=1,...,n\ .\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (6.5)$

Оценка сходимости процесса производится на основе вычисленных приближений λ_i^(s+1) и λ_i^(s), а также проверки малости всех внедиагональных элементов матриц K_s+1 и M_s+1.

В простейшем случае вращения продолжают, пока по окончании итерации с номером (t+1) не выполнится условие

$https://www.cyberforum.ru/cgi-bin/latex.cgi?\left| \lambda _i^{(t+1)}-\lambda _i^{(t)}\right|<\varepsilon , \ \ \ i=1,...,n\ ,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (6.6)$

где ε - заданная малая погрешность вычислений (принять ε=10^-5).

Окончательно принимается λ_i=λ_i^(r+1), где r - номер последнего вращения, а матрица X=[x₁,...,x_n], состоящая из собственных векторов, представляется в виде произведения

X=P₁P₂...P_r.

@kkkk5300 · 29.08.2014, 14:59 **[ТС]**

Что-то никто не отвечает, у меня все правильно что-ли?

Выкладываю словесное описание алгоритма (так, как я его понял):

Кликните здесь для просмотра всего текста

Замечания:
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.

@Cyborg Drone · 29.08.2014, 18:41

Пока заметил непорядок с функцией Sign. Так, что ли, надо:

Delphi
1
2
3
4
5
6
7
8
function Sign(a: real): integer;
begin
  if (a > 0) 
    then Result := 1
    else if a < 0
      then Result := -1
      else Result := 0
end;

Не знаю...

@kkkk5300 · 29.08.2014, 21:50 **[ТС]**

Да, но функцию Sign я специально так написал.

Предположим, что $https://www.cyberforum.ru/cgi-bin/latex.cgi?K[i][i] * M[j][j] = M[i][i] * K[j][j]$ , значит $https://www.cyberforum.ru/cgi-bin/latex.cgi?a = 0$
Тогда, если функция Sign описана так, как описали вы, то $https://www.cyberforum.ru/cgi-bin/latex.cgi?z = 0$ и далее деление на нуль.

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

@Cyborg Drone · 29.08.2014, 22:04

Мда... с ζ нехорошо получается... Будет сегодня время, гляну, что там можно сделать в программе... Посижу, поотлаживаю...

@kkkk5300 · 30.08.2014, 09:31 **[ТС]**

Спасибо, что помогаете.

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

@Cyborg Drone · 30.08.2014, 11:47

Сообщение от kkkk5300

я не нашел ничего схожего в интернете

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

Сообщение от kkkk5300

Я так понял, что метод основан на диагонализации пар матриц

Нет. И не обязательно пар матриц. Просто в Вашей задаче используются две матрицы.

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

Метод обычно применяется для приведения матриц к треугольному, трёхдиагональному, диагональному виду.

@Cyborg Drone · 02.09.2014, 14:16

Не по теме:

Что-то никак не получается найти ошибку... Попробую сегодня ночью разобраться. Наверное, просто свой вариант напишу.

Новые блоги и статьи Все статьи Все блоги /
Советы по крайней бережливости. Внимание, это ОЧЕНЬ длинный пост. Programma_Boinc 28.12.2025 Советы по крайней бережливости. Внимание, это ОЧЕНЬ длинный пост. Налог на собак: https:/ / ********/ gallery/ V06K53e Финансовый отчет в Excel: https:/ / ********/ gallery/ bKBkQFf Пост отсюда. . .	Кто-нибудь знает, где можно бесплатно получить настольный компьютер или ноутбук? США. Programma_Boinc 26.12.2025 Нашел на реддите интересную статью под названием Anyone know where to get a free Desktop or Laptop? Ниже её машинный перевод. После долгих разбирательств я наконец-то вернула себе. . .	Thinkpad X220 Tablet — это лучший бюджетный ноутбук для учёбы, точка. Programma_Boinc 23.12.2025 Рецензия / Мнение/ Перевод Нашел на реддите интересную статью под названием The Thinkpad X220 Tablet is the best budget school laptop period . Ниже её машинный перевод. Thinkpad X220 Tablet —. . .	PhpStorm 2025.3: WSL Terminal всегда стартует в ~ and_y87 14.12.2025 PhpStorm 2025. 3: WSL Terminal всегда стартует в ~ (home), игнорируя директорию проекта Симптом: После обновления до PhpStorm 2025. 3 встроенный терминал WSL открывается в домашней директории. . .	Как объединить две одинаковые БД Access с разными данными VikBal 11.12.2025 Помогите пожалуйста !! Как объединить 2 одинаковые БД Access с разными данными.
Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro	Музыка, написанная Искусственным Интеллектом volvo 05.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .	От async/await к виртуальным потокам в Python IndentationError 23.11.2025 Армин Ронахер поставил под сомнение async/ await. Создатель Flask заявляет: цветные функции - провал, виртуальные потоки - решение. Не threading-динозавры, а новое поколение лёгких потоков. Откат?. . .	Поиск "дружественных имён" СОМ портов Argus19 22.11.2025 Поиск "дружественных имён" СОМ портов На странице: https:/ / norseev. ru/ 2018/ 01/ 04/ comportlist_windows/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .	Сколько Государство потратило денег на меня, обеспечивая инсулином. Programma_Boinc 20.11.2025 Сколько Государство потратило денег на меня, обеспечивая инсулином. Вот решила сделать интересный приблизительный подсчет, сколько государство потратило на меня денег на покупку инсулинов. . . .

@kkkk5300 0 / 0 / 0 Регистрация: 04.07.2014 Сообщений: 4
	29.08.2014, 14:59 [ТС]
	Что-то никто не отвечает, у меня все правильно что-ли? Выкладываю словесное описание алгоритма (так, как я его понял): Кликните здесь для просмотра всего текста Замечания: 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. 0

@kkkk5300 0 / 0 / 0 Регистрация: 04.07.2014 Сообщений: 4
	29.08.2014, 21:50 [ТС]
	Да, но функцию Sign я специально так написал. Предположим, что $https://www.cyberforum.ru/cgi-bin/latex.cgi?K[i][i] * M[j][j] = M[i][i] * K[j][j]$ , значит $https://www.cyberforum.ru/cgi-bin/latex.cgi?a = 0$ Тогда, если функция Sign описана так, как описали вы, то $https://www.cyberforum.ru/cgi-bin/latex.cgi?z = 0$ и далее деление на нуль. Я все же попробовал написать функцию Sign так, как надо, и при делении на нуль переходить к следующей итерации. Но таким образом даже для тех матриц, для которых считало правильно, теперь считает неверно. 0

@Cyborg Drone Модератор 10373 / 5661 / 3398 Регистрация: 17.08.2012 Сообщений: 17,298
	29.08.2014, 22:04
	Мда... с ζ нехорошо получается... Будет сегодня время, гляну, что там можно сделать в программе... Посижу, поотлаживаю... 0

@kkkk5300 0 / 0 / 0 Регистрация: 04.07.2014 Сообщений: 4
	30.08.2014, 09:31 [ТС]
	Спасибо, что помогаете. А вообще, метод странный, я не нашел ничего схожего в интернете. Я так понял, что метод основан на диагонализации пар матриц. Так вот и про диагонализацию через именно такую матрицу я ничего не нашел. И возникает вопрос, правильно ли подобрана матрица для диагонализации. 0

@Cyborg Drone
	02.09.2014, 14:16
	Не по теме: Что-то никак не получается найти ошибку... Попробую сегодня ночью разобраться. Наверное, просто свой вариант напишу. 0