Нахождение собственных чисел и векторов методом Якоби и степенным методом - файл lab6_nm.doc
Нахождение собственных чисел и векторов методом Якоби и степенным методомДоступные файлы (2):
lab6_nm.doc
Завдання Знайти власні числа методом Якобі та обчислити відповідні власні вектори. Обчислити степеневим методом найбільше та найменше за модулем власні числа матриці на відповідні їм власні вектори. Матриця:
6.7 1.26 0.81 1.18 1.26 3.72 1.3 0.16 0.81 1.3 5.88 2.1 1.18 0.16 2.1 5.66 Текст програми #include
#include
#include
#include
#include
#include
#define max 10
#define eps 0.00001
using namespace std;
void inp_mf (double a[max][max], int *n, char file[max])
{
//reading matrix from file
int flag = 0, i = 0, j = 0;
{
ifstream f(file);
if (!f)
{
cout << "\n\nFile not found.";
getchar();
getchar ();
exit (1);
}
f >> *n;
if (*n > max || *n < 1)
{
cout << "Matrix dimensions overflow.";
getchar();
exit (1);
};
for (i = 0; i < *n; i++)
for (j = 0; j < *n; j++)
{
f >> a[i][j];
if (fabsl(a[i][j]) > eps) flag = 1;
}
if (flag == 0)
{
cout << "Attention! Zero matrix.";
getchar();
exit (1);
}
for (i = 0; i < *n; i++) f >> a[i][*n];
f.close();
}
}
void inp_mk ()
{
int i = 0, j = 0, n = 0;
double t = 0;
//reading matrix from keyboard and writing to file
{
ofstream f("inp_matr.txt");
if (!f)
{
cout << "File not found and cannot be created.";
getchar();
exit (1);
}
cout << "Enter the dimensions of matrix -> ";
cin >> n;
if (n > max || n < 1)
{
cout << "Matrix dimensions overflow";
getchar();
exit (1);
};
f << n <<"\n";
for (i = 0; i < n; i++)
{
for (j = 0; j < n; j++)
{
cout << "a[" << i + 1 << "][" << j + 1 << "] -> ";
cin >> t;
f << t << ' ';
}
f << "\n";
}
f.close();
}
}
void printm (double a[max][max], int n, ofstream& f)
{
//printing matrix
int i = 0, j = 0;
for (i = 0; i < n; i++)
{
for (j = 0; j < n; j++) f << a[i][j] << " ";
f << "\n";
}
}
void mult_matr (double a[max][max], double b[max][max], int n, double ab[max][max])
{
int i=0, j=0, k;
for ( i=0; i
for ( j=0; j
{ ab[i][j] = 0;
for ( k=0; k
ab[i][j] = ab[i][j]+a[i][k]*b[k][j];
}
}
double normaF (double a[max][max], int n)
{
int i, j = 0;
double s = 0;
for (i = 0; i
for (j = 0; j
s = s + pow(a[i][j],2);
return(pow(s,0.5));
}
double delta (double a[max][max], int n)
{
int i = 0;
double s = 0;
for (i = 0; i
s = s + pow(a[i][i],2);
return(s);
}
double w (double a[max][max], int n)
{
int i, j = 0;
double s = 0;
for (i = 0; i
for (j = i+1; j
s = s + pow(a[i][j],2);
return(s);
}
void JacobiMethod(double a[max][max], int n)
{
ofstream f ("out_matr.txt");
f.precision (7);
f.width (14);
double alpha, beta, gama, dzeta, t, c, s, Q[max][max], Q1[max][max], A[max][max], a1[max][max];
int i, j, k, i_max, j_max = 0;
for (i=0; i
for (j=0; j
{
A[i][j] = a[i][j];
if ( i == j ) Q[i][j] = 1;
else Q[i][j] = 0;
}
f<<" NormaF = "<
int p = 0;
do
{
p++;
f<<"\n Iteration number "<
gama = 0;
for (i = 0; i
for (j = i+1; j
if (fabs(gama) < fabs(a[i][j]))
{ gama = a[i][j];
i_max = i;
j_max = j;
alpha = a[i][i];
beta = a[j][j];
}
dzeta = (beta - alpha)/(2*gama);
if (dzeta <= 0) t = 1/(dzeta - sqrt(dzeta*dzeta+1));
else t = 1/(dzeta + sqrt(dzeta*dzeta+1));
c = 1/sqrt(t*t+1);
s = t*c;
f<<"c = "<
f<<"s = "<
f<<"alpha (a["<
f<<"beta (a["<
f<<"gama (a["<
f<<"dzeta = "<
f<<"t = "<
for ( i=0; i
for ( j=0; j
{
a1[i][j] = a[i][j];
Q1[i][j] = Q[i][j];
}
for ( k=0; k
{
a1[k][i_max] = a[k][i_max]*c - a[k][j_max]*s;
}
for ( k=i_max+1; k<=j_max; k++)
{
a1[i_max][k] = a[i_max][k]*c - a[k][j_max]*s;
}
for ( k=j_max+1; k
{
a1[i_max][k] = a[i_max][k]*c - a[j_max][k]*s;
}
for ( k=0; k
{
a1[k][j_max] = a[k][i_max]*s + a[k][j_max]*c;
}
for ( k=i_max+1; k<=j_max; k++)
{
a1[k][j_max] = a[i_max][k]*s + a[k][j_max]*c;
}
for ( k=j_max+1; k
{
a1[j_max][k] = a[i_max][k]*s + a[j_max][k]*c;
}
a1[i_max][i_max] = a[i_max][i_max]*c*c + a[j_max][j_max]*s*s - 2*c*s*a[i_max][j_max];
a1[j_max][j_max] = a[i_max][i_max]*s*s + a[j_max][j_max]*c*c + 2*c*s*a[i_max][j_max];
a1[i_max][j_max] = 0;
for ( k=0; k
for ( i=k+1; i
a1[i][k] = a1[k][i];
for ( i=0; i
for ( j=0; j
a[i][j] = a1[i][j];
f<<" A"<
printm(a, n, f);
printm(Q, n, f);
for ( j=0; j
{
Q1[j][i_max] = c*Q[j][i_max] - s*Q[j][j_max];
Q1[j][j_max] = s*Q[j][i_max] + c*Q[j][j_max];
}
for ( i=0; i
for ( j=0; j
Q[i][j] = Q1[i][j];
f<<" Q"<
printm(Q, n, f);
f<<" Delta = "<
f<<" NormaF = "<
f<<"W = "<
}
while ( w ( a, n ) >= pow(eps,2) );
mult_matr(A,Q,n,a1);
for ( i=0; i
{
f<<"\n"<
f<<" eigenvector :";
for ( j=0; j
f<
f<<"\n Av-lv : ";
for ( j=0; j
f<
}
}
double norma (double x[max], int n)
{
double S = 0;
for (int i = 0; i < n; i++) S = S + pow (x[i], 2);
return pow (S, 0.5);
}
void stepen (double matr [max][max], int n, double *l, double x[max], int s)
{
double y[max], p[max];
int i, j, k = 0;
for (i = 0; i < n; i++) x[i] = 1;
int flag = 0;
while (flag == 0)
{
flag = 1;
for (i = 0; i < n; i++)
{
y[i] = 0;
for (j = 0; j < n; j++) y[i] = y[i] + x[j]*matr[i][j];
}
*l = norma (y, n)/norma(x, n);
if (x[0]*y[0] < 0) *l = -*l;
for (i = 0; i < n; i++) if (fabsl(x[i]-y[i]/(*l)) > eps) flag = 0;
for (i = 0; i < n; i++) x[i] = y[i]/norma(y, n);
}
}
int main()
{
int n = 0, flag = 0, i = 0, j = 0, k = 0;
double a[max][max], A[max][max];
cout << "Please, enter what do you want to do:\n 1. Read matrix slar from file.\n 2. Enter slar from keyboard.\n - > ";
cin >> flag;
char inp[] = "inp_matr.txt";
if (flag == 1) inp_mf (a, &n, inp);
else
{
inp_mk ();
inp_mf (a, &n, inp);
}
for ( i=0; i
for ( j=0; j
A[i][j] = a[i][j];
JacobiMethod(a, n);
cout << "\n\n Power method \n";
double v1[max], v2[max];
for ( i=0; i
v1[i] = 1;
double lmax = 0;
stepen (A, n, &lmax, v1, 1);
for (i = 0; i < n; i++)
{
v2[i] = 0;
for (j = 0; j < n; j++) v2[i] = v2[i] + A[i][j]*v1[j];
}
cout << "\n Proper number_max = " << lmax;
cout << "\n Eigenvector:";
for (i = 0; i < n; i++) cout << v1[i] << " ";
cout << "\n Residual vector:";
for (i = 0; i < n; i++) cout << v2[i] - v1[i]*lmax << " ";
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
{
if (i != j) a[i][j] = A[i][j];
else a[i][j] = A[i][j] - lmax;
}
double lmin = 0;
stepen (a, n, &lmin, v1, 2);
for (i = 0; i < n; i++)
{
v2[i] = 0;
for (j = 0; j < n; j++) v2[i] = v2[i] + A[i][j]*v1[j];
}
lmin = norma (v2, n)/norma(v1, n);
if (v1[0]*v2[0] < 0) lmin = -lmin;
cout << "\n\n Proper number_min = " << lmin;
cout << "\n Eigenvector:";
for (i = 0; i < n; i++) cout << v1[i] << " ";
cout << "\n Residual vector:";
for (i = 0; i < n; i++) cout << v2[i] - v1[i]*lmin << " ";
getchar();
}
Результати роботи
Power method
Proper number_max = 9.25747
Eigenvector:0.562094 0.277295 0.571114 0.530081
Residual vector:-4.9032e-005 -1.03899e-005 2.92896e-005 2.58737e-005
Proper number_min = 2.57038
Eigenvector:0.244928 -0.754764 0.476611 -0.378395
Residual vector:8.01378e-006 -3.41725e-005 -2.73742e-005 3.88716e-005
Jacobi method
1. proper number 5.849002
eigenvector :0.7802453 0.1369973 -0.4812273 -0.3753256
Av-lv : 0.003925614 0.005793835 0.006994363 0.002968846
2. proper number 2.570739
eigenvector :-0.2474624 0.7533496 -0.4797225 0.3756235
Av-lv : -0.01795573 -0.009546406 -0.01980769 -0.01848757
3. proper number 9.261516
eigenvector :0.5552177 0.2949443 0.5812657 0.5165958
Av-lv : 0.02983925 -0.09656342 -0.04755405 0.06247804
4. proper number 4.278744
eigenvector :0.1473532 -0.5715802 -0.4476793 0.6716892
Av-lv : 0.06656181 0.03051914 0.07 0.07007236