Форум программистов, компьютерный форум CyberForum.ru

метод вращений - C++

Восстановить пароль Регистрация
 
Рейтинг: Рейтинг темы: голосов - 25, средняя оценка - 4.96
papirus
 Аватар для papirus
2 / 2 / 0
Регистрация: 21.07.2009
Сообщений: 49
05.12.2010, 18:31     метод вращений #1
нужно найти собственные значения и векторы: вот прога тока она кажется путает индесы элемента http://www.cyberforum.ru/cgi-bin/latex.cgi?a_{ij} из-за этого не правильно считает, пытался другую матрицу задать, она вообще зацикливается
C++
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
#include <iostream>
#include <math.h>
#include <stdio.h>
#define eps 0.001
const int n = 2;
double max_el(int& i, int& j);
void multi_m(double m[n][n], double p[n][n], double (&t)[n][n] );
void print_m(double m[n][n]);
void tr_m(double (&m)[n][n]);
void sled(double m[n][n]);
 
double a[n][n] = {
/*  {1.8, 1.6, 1.7, 1.8}, 
    {1.6, 2.8, 1.5, 1.3},
    {1.7, 1.5, 3.8, 1.4},
    {1.8, 1.3, 1.4, 4.3}
    {5,1,2},                //если на этом примере, то прога зацикливается
    {1,4,1},
    {2,1,3}
*/
    {2,1},   //а на этом вообще не работает
    {1,3}
    };
int main() {
int im , jm;
int p = 1;
printf("\n Матрица A:");
print_m(a);
while (fabs(max_el(im, jm)) > eps) {
    printf("\nИтерация: %d" ,p);
    sled(a);
 
    double maxij = max_el(im,jm);
    printf("\nmax:a[%d,%d]= %4.8f", im, jm, maxij);
    double fi = atan(2*maxij/(a[im][im]-a[jm][jm]))/2;
    printf("\nfi=  %4.8f",fi);
    double cosfi = cos(fi), sinfi = sin(fi);
    printf("\ncosfi=%4.8f sinfi=%4.8f",cosfi,sinfi);
    double hi[n][n];
 
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++) {
            if (i == j)
                hi[i][j] = 1;
            else hi[i][j] = 0;
        }
    hi[im][im] = cosfi;
    hi[jm][jm] = cosfi;
    hi[im][jm] = -sinfi;
    hi[jm][im] = sinfi;
 
    printf("\nматрица вращения H:");
    print_m(hi);
 
    double hic[n][n];  //вспомогательная матрица для умножения
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            hic[i][j] = hi[i][j];
 
    tr_m(hi);
 
    double b[n][n];
    multi_m(hi, a, b);
    multi_m(b, hic, a);
    printf("\nСледующее приближение:");
    print_m(a);
    ++p;
}
}
 
void sled(double m[n][n]) {
    double s = 0;
    for (int i = 0; i < n; i++)
        s += a[i][i];
    printf("\nслед= %f",s);
};
 
double max_el(int& f, int& g) {
double max = a[0][1];
for (int j = 1; j < n; j++)
    for (int i = 0; i < j; i++) {
        if (a[i][j] > max) {
            max = a[i][j];
            f = i; g = j;
        }
    }
return max;
 }
 
void tr_m(double (&m)[n][n]) {
    for (int i = 1; i < n; i++)
        for (int j = 0; j < i; j++) {
            double buf = m[i][j];
            m[i][j] = m[j][i];
            m[j][i] = buf;
        }
}
 
void multi_m(double m[n][n], double p[n][n], double (&t)[n][n]) {
    for (int i = 0; i < n; i++)
    for (int j = 0; j < n; j++) {
            double s = 0;
            for (int l = 0; l < n; l++)
                s += m[i][l]*p[l][j];
            t[i][j] = s;
    }
 }
 
void print_m(double m[n][n]) {
    printf("\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++)
            printf("%4.8f ", m[i][j]);
    printf("\n");
    };
 }
После регистрации реклама в сообщениях будет скрыта и будут доступны все возможности форума.
papirus
 Аватар для papirus
2 / 2 / 0
Регистрация: 21.07.2009
Сообщений: 49
08.12.2010, 15:39  [ТС]     метод вращений #2
помогите плз
Кирилл7785
Сообщений: n/a
26.02.2011, 19:05     метод вращений #3
Замени код своей функции с таким же именем следующей

double max_el(int& f, int& g) {
double max = a[0][1];
f=0; g=1; // стартовое значение
for (int j = 1; j < n; j++)
for (int i = 0; i < j; i++) {
if (a[i][j] > max) {
max = a[i][j];
f = i; g = j;
}
}
return max;
}

Ошибка была в том что нужно было присвоить начальные значения. А вообще спасибо большое за код я буду его использовать.
011
9 / 9 / 0
Регистрация: 28.11.2013
Сообщений: 149
06.12.2015, 20:15     метод вращений #4
Кто подскажет, в чем ошибка этого кода? (с учетом замечания, сделанного Кирилл7785)
Дело в том, что точность не соблюдается.
Пусть дана симметричная матрица
Код
{1.8, 1.6, 1.7, 1.8},
{1.6, 2.8, 1.5, 1.3},
{1.7, 1.5, 3.8, 1.4},
{1.8, 1.3, 1.4, 4.3}
На 9ой итерации данный код остановится и выведет такую матрицу:
Код
0,35821852 0,00000000 -0,08988536 -0,06623520
0,00000000 1,71300058 -0,00187514 -0,00027005
-0,08988536 -0,00187514 2,70509201 0,00000000
-0,06623520 -0,00027005 0,00000000 7,92368889
При заданной точности eps = 0.001.
Предполагается, что собственные значения на главной диагонали.

Однако Wolfram Mathematica 10.0 говорит о других числах (с точностью до первых 9 знаков после запятой):
{7.924268858, 2.708531888, 1.712997043, 0.3542022109}.

Как видно, точность не соблюдена.
Я специально привел пример с относительно небольшой точностью (вдруг, дело в double).
Однако, я переписал такой же код, только с использованием длинной арифметики, результат такой же.

Не по теме:

(Зря вечер убил, стоило более обстоятельно проверить код)



Условием продолжения процесса, как я вижу, служит условие
Код
fabs(max_el(im, jm)) > eps
Однако, в приведенном мною выше примере (с матрицей), на 9-ой итерации (последней) выводится такое значение
Код
max:a[0,1]= 0,02782967
Что никак не меньше eps.
Полный вывод 9ой итерации

Код
Итерация: 9
след= 12,700000
max:a[0,1]= 0,02782967
fi=  -0,02054759
cosfi=0,99978891 sinfi=-0,02054614
матрица вращения H:
0,99978891 0,02054614 0,00000000 0,00000000
-0,02054614 0,99978891 0,00000000 0,00000000
0,00000000 0,00000000 1,00000000 0,00000000
0,00000000 0,00000000 0,00000000 1,00000000

Следующее приближение:
0,35821852 0,00000000 -0,08988536 -0,06623520
0,00000000 1,71300058 -0,00187514 -0,00027005
-0,08988536 -0,00187514 2,70509201 0,00000000
-0,06623520 -0,00027005 0,00000000 7,92368889


Я немного поменял порядок вычислений, вроде матрица вращения на последней итерации стала "приятнее", но значения точнее не стали. Кстати, на одну итерацию стало больше.
Код
...
double maxij = max_el(im, jm);
	while (fabs(maxij) > eps) {
...
}
Полный вывод 10ой итерации

Код
Итерация: 10
след= 12,700000
max:a[2,3]= 0,00000000
fi=  -0,00000000
cosfi=1,00000000 sinfi=-0,00000000
матрица вращения H:
1,00000000 0,00000000 0,00000000 0,00000000
0,00000000 1,00000000 0,00000000 0,00000000
0,00000000 0,00000000 1,00000000 0,00000000
0,00000000 0,00000000 -0,00000000 1,00000000

Следующее приближение:
0,35821852 0,00000000 -0,08988536 -0,06623520
0,00000000 1,71300058 -0,00187514 -0,00027005
-0,08988536 -0,00187514 2,70509201 -0,00000000
-0,06623520 -0,00027005 0,00000000 7,92368889


Буду благодарен за любую помощь!
Возможно, Вы знаете правильную реализацию этого алгоритма и можете поделиться ссылкой? (язык реализации не важен, если в коде нет сложных конструкций)
Yandex
Объявления
06.12.2015, 20:15     метод вращений
Ответ Создать тему
Опции темы

Текущее время: 21:11. Часовой пояс GMT +3.
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2016, vBulletin Solutions, Inc.
Рейтинг@Mail.ru