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

Метод вращений (Якоби) - C++

Восстановить пароль Регистрация
 
Рейтинг: Рейтинг темы: голосов - 19, средняя оценка - 4.74
Fantom.AS
 Аватар для Fantom.AS
2 / 1 / 0
Регистрация: 17.11.2010
Сообщений: 121
20.12.2011, 17:51     Метод вращений (Якоби) #1
Почему то данный алгоритм работает не правильно. На часть тестов дает левый ответ.
Подскажите, пожалуйста, в чем ошибка?
например,
на матрице
3 0
2 -1

выдает
d = 3 1 //правильно
v =
1 0
0 1 // не правильно

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
//***************************************************************************//
// 
//  Метод Якоби нахождения собственных значений и собственных векторов
//                       для симметричных матриц
//    Входные данные:
//        n - размерность матрицы
//        a - исходная матрица. В процессе работы наддиагональные элементы
//            будут изменены, но их легко восстановить по поддиагональным
//    Выходные данные:
//        d - массив собственных значений
//        v - массив собственных векторов
//
//***************************************************************************//
 
void jacobi ( const int n, double * const * a, double * d, double * const * v )
{
    if ( n == 0 ) return;
    double * b = new double[n+n];
    double * z = b + n;
    unsigned int i, j;
    for ( i = 0; i < n; ++i )
    {
        z[i] = 0.;
        b[i] = d[i] = a[i][i];
        for ( j = 0; j < n; ++j ) v[i][j] = i == j ? 1. : 0.;
    }
    for ( i = 0; i < 50; ++i )
    {
        double sm = 0.;
        unsigned int p, q;
        for ( p = 0; p < n - 1; ++p )
        {
            for ( q = p + 1; q < n; ++q ) sm += fabs ( a[p][q] );
        }
        if ( sm == 0 ) break;
        const double tresh = i < 3 ? 0.2 * sm / ( n*n ) : 0.;
        for ( p = 0; p < n - 1; ++p )
        {
            for ( q = p + 1; q < n; ++q )
            {
                const double g = 1e12 * fabs ( a[p][q] );
                if ( i >= 3 && fabs ( d[p] ) > g && fabs ( d[q] ) > g ) a[p][q] = 0.;
                else
                if ( fabs ( a[p][q] ) > tresh )
                {
                    const double theta = 0.5 * ( d[q] - d[p] ) / a[p][q];
                    double t = 1. / ( fabs(theta) + sqrt(1.+theta*theta) );
                    if ( theta < 0 ) t = - t;
                    const double c = 1. / sqrt ( 1. + t*t );
                    const double s = t * c;
                    const double tau = s / ( 1. + c );
                    const double h = t * a[p][q];
                    z[p] -= h;
                    z[q] += h;
                    d[p] -= h;
                    d[q] += h;
                    a[p][q] = 0.;
                    for ( j = 0; j < p; ++j )
                    {
                        const double g = a[j][p];
                        const double h = a[j][q];
                        a[j][p] = g - s * ( h + g * tau );
                        a[j][q] = h + s * ( g - h * tau );
                    }
                    for ( j = p+1; j < q; ++j )
                    {
                        const double g = a[p][j];
                        const double h = a[j][q];
                        a[p][j] = g - s * ( h + g * tau );
                        a[j][q] = h + s * ( g - h * tau );
                    }
                    for ( j = q+1; j < n; ++j )
                    {
                        const double g = a[p][j];
                        const double h = a[q][j];
                        a[p][j] = g - s * ( h + g * tau );
                        a[q][j] = h + s * ( g - h * tau );
                    }
                    for ( j = 0; j < n; ++j )
                    {
                        const double g = v[j][p];
                        const double h = v[j][q];
                        v[j][p] = g - s * ( h + g * tau );
                        v[j][q] = h + s * ( g - h * tau );
                    }
                }
            }
        }
        for ( p = 0; p < n; ++p )
        {
            d[p] = ( b[p] += z[p] );
            z[p] = 0.;
        }
    }
}
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
20.12.2011, 17:51     Метод вращений (Якоби)
Посмотрите здесь:

C++ метод вращений
C++ метод якоби
C++ Метод итерации( Якоби)
C++ Метод вращений Якоби с++
Метод вращений с построением КЮЭР-разложения C++
C++ КУЭР-разложение методом вращений
Метод Якоби. Выводит результат -1.INF и -1.IND C++
C++ Итерационные методы. Метод Якоби

Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
После регистрации реклама в сообщениях будет скрыта и будут доступны все возможности форума.
Ответ Создать тему
Опции темы

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