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

Найти обратную матрицу методом Гауса - C++

Восстановить пароль Регистрация
 
stawerfar
 Аватар для stawerfar
141 / 55 / 4
Регистрация: 14.12.2010
Сообщений: 347
Записей в блоге: 1
19.06.2012, 20:42     Найти обратную матрицу методом Гауса #1
Всем привет. Задание следующее,есть матрица
C++
1
2
3
4
5
typedef GLdouble GLTDoubleMatrix[4][4];
GLTDoubleMatrix tempm =  {-2,  1,  3,  2,
                       1,  2, -1,  1,
                       1, -2, -1,  2,
                       2,  2,  1,  1};
Нужно найти обратную матрицу, я это сделал , частичный код приведен ниже.
Проблема заключается в том что при умножении матрицы на обратную матрицу не получается единичная матрица. Компилятор при сложении чисел в некоторых местах теряет разряды и в итоге получается число типа 5.51342343e + 017. Например при сложении чисел 0.39130434782608697 и
-0.39130434782608692 за место 0 получается 5.5511151231257827e-017.
Да вы можете заметить что числа отличаются в последнем разряде 7 и 2 , но я тут не причастен компилятор так считает.
Или причастен?
Кто сталкивался подскажите пожалуйста.
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
117
118
119
120
121
122
123
124
125
126
127
128
129
double gltDeterminantMatrix(const GLTDoubleMatrix const m)
{
 
    //разложим матрицу по по строкам и столбцам
    const int ciMAX = 4;
    double d1 = 0.0,d2 = 0.0,d3 = 0.0,d4 = 0.0;
 
                
    //подсчитаем определитель первой матрицы
    d1 = m[0][0]*((m[1][1] * m[2][2] *m[3][3] + m[1][2] * m[2][3] * m[3][1] + m[2][1] * m[3][2] * m[1][3])   - (m[3][1] * m[2][2] * m[1][3] + m[1][1] * m[3][2] * m[2][3] + m[2][1] * m[1][2] * m[3][3]));
    //подсчитаем определитель второй матрицы
    d2 = -m[0][1]*((m[1][0] * m[2][2] * m[3][3] + m[2][0] * m[1][3] * m[3][2] + m[1][2] * m[2][3] * m[3][0]) - (m[3][0] * m[2][2] * m[1][3] + m[3][2] * m[2][3] * m[1][0] + m[2][0] * m[1][2] * m[3][3]));
    //подсчитаем определитель третьей матрицы
    d3 = m[0][2]*((m[1][0] * m[2][1] * m[3][3] + m[2][0] * m[1][3] * m[3][1] + m[1][1] * m[2][3] * m[3][0])  - (m[3][0] * m[2][1] * m[1][3] + m[3][1] * m[2][3] * m[1][0] + m[2][0] * m[1][1] * m[3][3]));
    //подсчитаем определитель четверной матрицы
    d4 = -m[0][3]*((m[1][0] * m[2][1] * m[3][2] + m[2][0] * m[1][2] * m[3][1] + m[3][0] * m[1][1] * m[2][2]) - (m[3][0] * m[2][1] * m[1][2] + m[3][1] * m[2][2] * m[1][0] + m[2][0] * m[1][1] * m[3][2]));
    
    //подсчитать детерминант матрицы 4х4
    return d1+d2+d3+d4;
}
 
//нахождение обратной матрицы
bool gltInverseMatrix(const GLTDoubleMatrix const m,GLTDoubleMatrix mResult)
{
    //нахождение детерминанта
    double det= gltDeterminantMatrix(m);
    //если детерминант матрицы не равен 0 то значит обратная матрицы существует
    
    const int ciMAX = 4;
 
    if( det != 0)
    {
 
        //тогда найдем алгиброические дополнения
 
            //A11 в матрице  вычеркиваем строку 1 и столбец 1.
            mResult[0][0] =  ((m[1][1] * m[2][2] *m[3][3] + m[1][2] * m[2][3] * m[3][1] + m[2][1] * m[3][2] * m[1][3])  - (m[3][1] * m[2][2] * m[1][3] + m[1][1] * m[3][2] * m[2][3] + m[2][1] * m[1][2] * m[3][3]));
            //A12 в матрице  вычеркиваем строку 1 и столбец 2.
            mResult[1][0] = -((m[1][0] * m[2][2] * m[3][3] + m[2][0] * m[1][3] * m[3][2] + m[1][2] * m[2][3] * m[3][0]) - (m[3][0] * m[2][2] * m[1][3] + m[3][2] * m[2][3] * m[1][0] + m[2][0] * m[1][2] * m[3][3]));
            //A13 в матрице  вычеркиваем строку 1 и столбец 3.
            mResult[2][0] =  ((m[1][0] * m[2][1] * m[3][3] + m[2][0] * m[1][3] * m[3][1] + m[1][1] * m[2][3] * m[3][0]) - (m[3][0] * m[2][1] * m[1][3] + m[3][1] * m[2][3] * m[1][0] + m[2][0] * m[1][1] * m[3][3]));
            //A14 в матрице  вычеркиваем строку 1 и столбец 4.
            mResult[3][0] = -((m[1][0] * m[2][1] * m[3][2] + m[2][0] * m[1][2] * m[3][1] + m[3][0] * m[1][1] * m[2][2]) - (m[3][0] * m[2][1] * m[1][2] + m[3][1] * m[2][2] * m[1][0] + m[2][0] * m[1][1] * m[3][2]));
            
            //A21 в матрице  вычеркиваем строку 2 и столбец 1.
            mResult[0][1] =  -((m[0][1] * m[2][2] * m[3][3] + m[2][1] * m[0][3] * m[3][2] + m[0][2] * m[2][3] * m[3][1]) - (m[3][1] * m[2][2] * m[0][3] + m[3][2] * m[2][3] * m[0][1] + m[2][1] * m[0][2] * m[3][3]));
            //A22 в матрице  вычеркиваем строку 2 и столбец 2.
            mResult[1][1] =   ((m[0][0] * m[2][2] * m[3][3] + m[2][0] * m[3][2] * m[0][3] + m[3][0] * m[0][2] * m[2][3]) - (m[3][0] * m[2][2] * m[0][3] + m[0][0] * m[3][2] * m[2][3] + m[2][0] * m[0][2] * m[3][3]));
            //A23 в матрице  вычеркиваем строку 2 и столбец 3.
            mResult[2][1] =  -((m[0][0] * m[2][1] * m[3][3] + m[2][0] * m[3][1] * m[0][3] + m[2][3] * m[0][1] * m[3][0]) - (m[3][0] * m[2][1] * m[0][3] + m[0][0] * m[3][1] * m[2][3] + m[3][3] * m[2][0] * m[0][1]));
            //A24 в матрице  вычеркиваем строку 2 и столбец 4.
            mResult[3][1] =   ((m[0][0] * m[2][1] * m[3][2] + m[2][0] * m[3][1] * m[0][2] + m[2][2] * m[0][1] * m[3][0]) - (m[0][2] * m[2][1] * m[3][0] + m[2][0] * m[0][1] * m[3][2] + m[3][1] * m[2][2] * m[0][0])); 
 
            //A31 в матрице  вычеркиваем строку 3 и столбец 1.
            mResult[0][2] =   ((m[0][1] * m[1][2] * m[3][3] + m[1][1] * m[3][2] * m[0][3] + m[1][3] * m[0][2] * m[3][1]) - (m[3][1] * m[1][2] * m[0][3] + m[3][2] * m[1][3] * m[0][1] + m[1][1] * m[0][2] * m[3][3]));
            //A32 в матрице  вычеркиваем строку 3 и столбец 2.
            mResult[1][2] =  -((m[0][0] * m[1][2] * m[3][3] + m[1][0] * m[3][2] * m[0][3] + m[1][3] * m[0][2] * m[3][0]) - (m[3][0] * m[1][2] * m[0][3] + m[0][0] * m[3][2] * m[1][3] + m[1][0] * m[0][2] * m[3][3]));
            //A33 в матрице  вычеркиваем строку 3 и столбец 3.
            mResult[2][2] =   ((m[0][0] * m[1][1] * m[3][3] + m[1][0] * m[3][1] * m[0][3] + m[1][3] * m[0][1] * m[3][0]) - (m[3][0] * m[1][1] * m[0][3] + m[3][1] * m[1][3] * m[0][0] + m[1][0] * m[0][1] * m[3][3]));
            //A34 в матрице  вычеркиваем строку 3 и столбец 4.
            mResult[3][2] =  -((m[0][0] * m[1][1] * m[3][2] + m[1][0] * m[3][1] * m[0][2] + m[1][2] * m[0][1] * m[3][0]) - (m[3][0] * m[1][1] * m[0][2] + m[3][1] * m[1][2] * m[0][0] + m[1][0] * m[0][1] * m[3][2]));
            
            //A41 в матрице  вычеркиваем строку 4 и столбец 1.
            mResult[0][3] =  -((m[0][1] * m[1][2] * m[2][3] + m[1][1] * m[2][2] * m[0][3] + m[1][3] * m[0][2] * m[2][1]) - (m[2][1] * m[1][2] * m[0][3] + m[2][2] * m[1][3] * m[0][1] + m[1][1] * m[0][2] * m[2][3]));
            //A42 в матрице  вычеркиваем строку 4 и столбец 2.
            mResult[1][3] =   ((m[0][0] * m[1][2] * m[2][3] + m[1][0] * m[2][2] * m[0][3] + m[1][3] * m[0][2] * m[2][0]) - (m[2][0] * m[1][2] * m[0][3] + m[2][2] * m[1][3] * m[0][0] + m[1][0] * m[0][2] * m[2][3]));
            //A43 в матрице  вычеркиваем строку 4 и столбец 3.
            mResult[2][3] =  -((m[0][0] * m[1][1] * m[2][3] + m[1][0] * m[2][1] * m[0][3] + m[1][3] * m[0][1] * m[2][0]) - (m[2][0] * m[1][1] * m[0][3] + m[2][1] * m[1][3] * m[0][0] + m[1][0] * m[0][1] * m[2][3]));
            //A44 в матрице  вычеркиваем строку 4 и столбец 4.
            mResult[3][3] =   ((m[0][0] * m[1][1] * m[2][2] + m[1][0] * m[2][1] * m[0][2] + m[1][2] * m[0][1] * m[2][0]) - (m[2][0] * m[1][1] * m[0][2] + m[2][1] * m[1][2] * m[0][0] + m[1][0] * m[0][1] * m[2][2])); 
 
        //умножим получившуюся матрицу на 1/ det
            
            det = 1/det;
 
            for(int i=0; i< ciMAX ;i++)
            {
                for(int j = 0; j< ciMAX ; j++)
                {
                    mResult[i][j]*=det;
                }
            }
    }
    else
    {
        MessageBox(NULL,"Для текущей матрицы не существует обртатной матрицы","Ошибка",MB_OK);
        return true;
    }
}
 
//умножение матриц
void gltMatrixMultiplication(const GLTDoubleMatrix const m1, const GLTDoubleMatrix const m2, GLTDoubleMatrix mResult)
{
    const int ciMAX = 4;
    //обнуление выходящего  массива
    for(int i=0;i< ciMAX ;i++)
    {
        for(int k =0;k<ciMAX;k++)
        {
            mResult[i][k] = 0.0;    
        }
    }
 
    for(int i=0;i< ciMAX ;i++)
    {
        for(int k =0;k<ciMAX;k++)
        {
            for(int j =0;j<ciMAX;j++)
            {
                mResult[i][k]+= m1[i][j] * m2[j][k];
            }
 
        }
    }
}
 
int main(void)
{              typedef GLdouble GLTDoubleMatrix[4][4];
               
        
GLTDoubleMatrix tempm =  {-2,  1,  3,  2,
                       1,  2, -1,  1,
                       1, -2, -1,  2,
                       2,  2,  1,  1};
                GLTDoubleMatrix res,one;
        gltInverseMatrix(tempm,res);
                gltMatrixMultiplication(tempm,res,one);
return 0;
}
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
19.06.2012, 20:42     Найти обратную матрицу методом Гауса
Посмотрите здесь:

Для матрицы а(n, n) найти обратную матрицу C++
C++ С помощью метода отражения найти обратную матрицу
Найти обратную матрицу и умножить ее на вектор C++
решение системы уравнений методом Гауса C++
C++ Как привести матрицу к треугольному виду по методу гауса?
Матрицы: определить обратную матрицу C++
C++ Найти обратную матрицу
C++ Как найти систему методом гауса

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

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

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