Определитель матрицы (+метод Зейделя)
04.12.2011, 08:58. Показов 1139. Ответов 0
Здравствуйте, сам метод Зейделя было реализовать легко, но возникла проблема с проверкой на определитель матрицы, я попытался воспользоваться шаблоном, но он не подключился, переписать в обычные процедуры тоже не особо удалось.
С шаблонами работаю в первый раз
ошибка идет на не соответствие типов в шаблоне и в коде основной программы
"det = GaussDeterminant<double>( a[i][j],n);"
пока что необходимо найти определитель, алгоритм используется следующий:
Квадратная матрица преобразуется методом Гаусса в треугольную верхнюю матрицу, после чего нахождение детерминанта сводится к произведению элементов главной диагонали.
| 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
130
131
132
133
134
135
136
137
138
| #include <stdafx.h>
#include <stdio.h>
#include <math.h>
#include <conio.h>
#define eps 0.0001
using namespace std;
//условие выхода из счетчика
bool converge(double *xk, double *xkp, int n)
{
for (int i = 0; i < n; i++)
{
if (fabs(xk[i] - xkp[i]) >= eps)
return false;
}
return true;
}
//-------------
template<class T>
T GaussDeterminant(T M, int N) //в исходнике **const M, но он также не рабочий
{
T determinant(1);
int exchanges(0);
for(int l1(0); l1<N-1; ++l1)
{ //Перебираю все строки матрицы, кроме последней
int maxN( l1 );
T maxValue( fabs( M[l1][l1] ) );
for(int l2(l1+1); l2<N; ++l2)
{ //Нахожу строку с максимальным по модулю элементом
T const value( fabs( M[l2][l1] ) );
if( value > maxValue ) { maxN=l2; maxValue=value; }
}
if( maxN > l1 )
{ //Нужен обмен
T *const temp(M[l1]); M[l1]=M[maxN]; M[maxN]=temp;
++exchanges;
} else { //Обмен не нужен, но нужна проверка на ноль
if( maxValue == T(0) ) return maxValue;
}
T const value( M[l1][l1] );
determinant *= value;
for(int l2(l1+1); l2<N; ++l2)
{ //Вычитаю строку из всех последующих
T const k( M[l2][l1]/value );
M[l2][l1] = T(0);
for(int c(l1+1); c<N; ++c) M[l2][c] -= M[l1][c]*k;
}
}
determinant *= M[N-1][N-1];
if(exchanges%2) return -determinant; else return determinant;
}
int main()
{
int i, j, n;
double det;
double a[sizeof(n)][sizeof(n)], b[sizeof(n)], p[sizeof(n)], x[sizeof(n)],c[sizeof(n)];
GaussDeterminant<double>;
printf("\nInput n value(dim of Ax=c):");
scanf_s("%d", &n);
if (n <= 0)
{
printf("Your number is <=0, in this program n will be only >0\n");
getch();
return 1;
}
//ввод матрицы и вектора
printf("Now input the matrix a(i,j),i,j=0,...,%d:\n", n - 1);
for (i = 0; i < n; i++)
{
printf("Input string #%d:\n",i);
for (j = 0; j < n; j++)
{
scanf_s("%lf", &a[i][j]);
}
}
det = GaussDeterminant<double>( a[i][j],n);
//проверка определителя
/* if (n == 2)
{
del = a[0][0]*a[1][1]-a[0][1]*a[1][0];
if (del <= 0)
{
printf("delta<=0\n");
getch();
return 1;
}
}
if (n == 3)
{
del = a[0][0]*a[1][1]*a[2][2]+a[0][0]*a[1][2]*a[2][1]+a[0][1]*a[1][0]*a[2][2]-a[0][2]*a[1][1]*a[2][0]-a[0][1]*a[1][2]*a[2][0]-a[1][0]*a[2][1]*a[0][2];
if (del <= 0)
{
printf("delta<=0\n");
getch();
return 1;
}
}*/
//---------------
printf("Now input the matrix b(i),i=0,...,%d:\n", n - 1);
for (i = 0; i < n; i++)
scanf_s("%lf", &b[i]);
//--------
x[j] = b[i];
//выносим главную диагональ матрицы в отдельный вектор
for (i=0; i<n; i++)
{
c[i]=a[i][i];
}
//-------
do
{
//нахождение x[i]
for (int i = 0; i < n; i++)
{
double var = 0;
for (int j = 0; j < n; j++)
if (j != i) var += (a[i][j] * x[j]);
p[i] = x[i];
x[i] = (b[i] - var) / c[i];
}
}
//сравнение x[i] с предыдущим значением p[i]
while (!converge(x, p, n));
//вывод
printf("Solve... x_i=\n");
for (i = 0; i < n; i++)
printf("%f\n", x[i]);
getch();
return 0;
//-----------
} |
|
| 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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
| #include <stdafx.h>
#include <conio.h>
#include <math.h>
#include <iomanip>
#include <iostream>
#include <vector>
using namespace std;
#define eps 0.000000001
bool converge(std::vector<double> xk, std::vector<double> xkp, int n) {
for (int i = 0; i < n; i++) {
if (fabs(xk[i] - xkp[i]) >= eps)
return false;
}
return true;
}
template<typename T>
T GaussDeterminant(std::vector<std::vector<T> > &array)
{
T determinant(1);
int exchanges(0);
for (size_t i(0); i < array.size() - 1; ++i) { //Перебираю все строки матрицы, кроме последней
size_t maxN(i);
T maxValue(fabs(array[i][i]));
for (size_t j(i + 1); j < array.size(); ++j) { //Нахожу строку с максимальным по модулю элементом
T const value(fabs(array[j][i]));
if (value > maxValue) {
maxN = j;
maxValue = value;
}
}
if (maxN > i) { //Нужен обмен
array[i].swap(array[maxN]);
++exchanges;
} else { //Обмен не нужен, но нужна проверка на ноль
if (maxValue == T(0))
return maxValue;
}
T const value(array[i][i]);
determinant *= value;
for (size_t l2(i + 1); l2 < array.size(); ++l2) { //Вычитаю строку из всех последующих
T const k(array[l2][i] / value);
array[l2][i] = T(0);
for (size_t c(i + 1); c < array.size(); ++c)
array[l2][c] -= array[i][c] * k;
}
}
determinant *= array[array.size() - 1][array.size() - 1];
if (exchanges % 2) {
determinant = (- determinant);
}
return determinant;
}
int main() {
int arraysize;
double number;
std::vector<double> current, b, c, x, p;
std::vector<std::vector<double> > a;
cout << "Input n value(dim of Ax=c): ";
cin >> arraysize;
if (arraysize <= 0) {
cout << "Your number is <= 0, in this program n will be only > 0" << endl;
getch();
return 1;
}
cout << "Input the matrix A(i,j), i,j=0,...,: " << (arraysize - 1) << ": " << endl;
for (int i = 0; i < arraysize; i++) {
current.erase (current.begin(),current.end());
cout << "Input line #" << i << ": ";
for (int j = 0; j < arraysize; j++) {
cin >> number;
current.push_back(number);
}
a.push_back(current);
}
cout << "Matrix A: " << endl;
for (size_t i = 0; i < a.size(); i++) {
for (size_t j = 0; j < a[i].size(); j++) {
cout << setw(9) << a[i][j];
}
cout << endl;
}
std::vector<std::vector<double> > d(a.size());
std::copy(a.begin(), a.end(), d.begin());
double det = GaussDeterminant<double>(d);
if (det<=0)
{
cout << "Determinant matrix is <=0, metod Zeidel cann't solve this Matrix" <<endl;
getch();
return 1;
}
cout << "Input the matrix B(i),i=0,...," << (arraysize - 1) << ": " << endl;
cout << "Input line: ";
for (size_t i = 0; i < a.size(); i++) {
cin >> number;
b.push_back(number);
}
cout << "Matrix B: " << endl;
for (size_t i = 0; i < a.size(); i++) {
cout << setw(9) << b[i];
}
cout << endl;
// выносим главную диагональ в отдельный вектор
for (size_t i = 0; i < a.size(); i++) {
c.push_back(a[i][i]);
}
// задаём необходимые вектора
for (size_t i = 0; i < a.size(); i++) {
x.push_back(0);
p.push_back(0);
}
for (size_t i = 0; i < a.size(); i++) {
x.push_back(b[i]);
}
do {
for (size_t i = 0; i < a.size(); i++) {
double var = 0;
for (size_t j = 0; j < a[i].size(); j++)
if (j != i) {
var += (a[i][j] * x[j]);
}
p[i] = x[i];
x[i] = (b[i] - var) / c[i];
}
} while (!converge(x, p, arraysize));
cout << "Solved Xi: " << endl;
for (size_t i = 0; i < a.size(); i++) {
cout << setw(9) << x[i];
}
getch();
cout << endl;
return 0;
} |
|
Мало ли вдруг кому пригодится...
0
|
(
