Точные методы решения систем линейных алгебраических уравнений. Схема Халецкого

Задача и результат выполнения

Код программы:
#include
#include
#include
#include
#include

using namespace std;

int const n=4;
ofstream fout;

int readar(int n, double** x) //Чтение массива с файла

{
ifstream matrix;
matrix.open("Test.txt");
for(int i=0;ifor(int j=0;jmatrix>>x[i][j];
fout<<"Read matrix from the file";
fout<for(int i=0;i{for(int j=0;jfout<fout<fout<}

double arrnorm (int n, double** x) // Поиск нормы матрицы

{
double b[n];
for(int i=0;ib[i]=1;
for(int j=0;jfor(int i=0;ib[j]+=fabs(x[i][j]); // Сумма столбцов матрицы
fout<<"Sum of columns"<for(int i=0;ifout<fout<double max=b[0];
for(int i=0;i{
if(maxmax=b[i];
}
fout<fout<return max;
}

void test_obr(int n, double** A, double** inv, double** q) // Проверка правильности нахождения обратной матрицы

{
double m;
for (int i=0;ifor (int j=0;j{
m=0;
for (int k=0;km=m+A[i][k]*inv[k][j];
q[i][j]=m;
}
}

int SLAU(double** matrica_a, int n, double* massiv_b, double* x) // Вспомогательная функция для нахождения обратной матрицы

{
int i,j,k,r;
double c, M, max, s, **a, *b;
a=new double *[n];
for (i=0;ia[i]=new double[n];
b=new double [n];
for (i=0;ifor (j=0;ja[i][j]=matrica_a[i][j];

for (i=0;i
b[i]=massiv_b[i];

for (k=0;k
{
max=fabs(a[k][k]);
r=k;
for (i=k+1;iif (fabs(a[i][k])>max)
{
max=fabs(a[i][k]);
r=i;
}

for (j=0;j
{
c=a[k][j];
a[k][j]=a[r][j];
a[r][j]=c;
}

c=b[k];
b[k]=b[r];

b[r]=c;
for (i=k+1;i{
for (M=a[i][k]/a[k][k],j=k;ja[i][j]-=M*a[k][j]; b[i]-=M*b[k];
}
}

if (a[n-1][n-1]==0) if(b[n-1]==0)

return -1;
else return -2;

else
{

for(i=n-1;i>=0;i--)
{
for(s=0,j=i+1;js+=a[i][j]*x[j];
x[i]=(b[i]-s)/a[i][i];
}
return 0;
}
for(i=0;idelete [] a[i];
delete [] a;
delete [] b;
}

int INVERSE(double** a, int n, double** y) // Алгоритм нахождения обратной матрицы

{
int i,j,res;
double *b, *x;
b=new double [n];
x=new double [n];
for (i=0;i{
for(j=0;jif (j==i)
b[j]=1;
else
b[j]=0;

res=SLAU(a,n,b,x);

if (res!=0) break;

else for (j=0;j
y[j][i]=x[j];
}

delete [] x;

delete [] b;

if (res!=0)

return -1;
else
return 0;
}

int main()

{

fout.open("Solution.txt");

int result;

double **x=new double* [n];

for(int i=0;ix[i]=new double [n+1];

double **b=new double* [n];

for(int i=0;ib[i]=new double [n];

double **c=new double* [n];

for(int i=0;ic[i]=new double [n+1];

double **unknown=new double* [n];

for(int i=0;iunknown[i]=new double [n];

double **Inv=new double* [n];

for(int i=0;iInv[i]=new double [n];

double **Unite=new double* [n];

for(int i=0;iUnite[i]=new double [n];

readar(n, x); // Чтение матрицы с файла

for(int i=0;i
for(int j=0;junknown[i][j]=x[i][j];

fout<<"Coefficients of unknown variables "<
for(int i=0;i{
for(int j=0;jfout<fout<}

fout<

for(int i=0;i
fout<fout<

double a[n][n+2], y[n], sol[n], det=1, mult[n], T1, T2;

for(int i=0;i
for(int j=0;ja[i][j]=0;

for(int i=0;i
for(int j=0;jif(j==n+1)
{
for(int l=0;la[i][n+1]+=x[i][l];
}
else
{
a[i][j]=x[i][j];
}

for(int i=0;i
for(int j=0;jb[i][j]=0;
for(int i=0;ifor(int j=0;jc[i][j]=0;

for(int i=0; i
{
for(int j=0; j{
b[i][0]=a[i][0];
c[0][i]=a[0][i]/b[0][0];
double sum=0;
for (int k=0; ksum+= b[i][k]*c[k][j];
b[i][j]=a[i][j]-sum;
if (ib[i][j] = 0;
if (i>j)
c[i][j]=0;
sum=0;
for (int k=0; ksum+=b[i][k]*c[k][j];
c[i][j]=(a[i][j]-sum)/b[i][i];
if (i>j)
c[i][j]=0;
}
}

fout<<"Lower triangle matrix"<
for(int i=0;i{
for(int j=0;jfout<fout<}

fout<
for(int i=0;i{
for(int j=0;jfout<fout<}
fout<

y[0]=a[0][n+1]/b[0][0];

for(int i=0;i{
double sum=0;
for (int k=0; ksum+= b[i][k]*y[k];
y[i]=(a[i][n]-sum)/b[i][i];
}
fout<<"Matrix Y"<for(int i=0;ifout<fout<

sol[n-1]=y[n-1]; // Подсчет матрицы решений

for(int i=n-2;i>=0;i--)
{
double sum=0;
for(int k=i+1; ksum+=c[i][k]*sol[k];
sol[i]=y[i]-sum;
}
fout<for(int i=0;ifout<fout<

for(int i=0; i
det*=b[i][i];
fout<

for(int j=0;j
{
mult[j]=0;
for(int i=0;imult[j]+=a[j][i]*sol[i];
}
fout<<"Checking the solution"<for(int i=0;ifout<fout<

fout<

result=INVERSE(unknown, n, Inv); // Подсчет обратной матрицы

fout<<"The inverse matrix of unknown variables "<

for(int i=0;i
{for(int j=0;jfout<fout<

test_obr(n, unknown, Inv, Unite); // Проверка обратной матрицы

fout<
for(int i=0;i{for(int j=0;jfout<fout<fout<double normA=arrnorm(n, x); // Подсчет нормы матрицы
double normInv=arrnorm(n, Inv);

fout<<"Condition number is equal to "<
printf("\n");
system("pause");
return 0;
}


Результаты работы тестовой программы:

Read matrix from the file

4 1 -2 -1 3
2 -3 1 3 14
-2 -1 3 0 2
1 2 -4 1 9

Coefficients of unknown variables

4 1 -2 -1
2 -3 1 3
-2 -1 3 0
1 2 -4 1

Intercepts

3 14 2 9

Lower triangle matrix

4 0 0 0
2 -3.5 0 0
-2 -0.5 1.7143 0
1 1.75 -2.5 1.5417

Upper triangle matrix

1 0.25 -0.5 -0.25
0 1 -0.57143 -1
0 0 1 -0.58333
0 0 0 1

Matrix Y
0.75 -3.5714 1 11.027

Solution matrix
4.2973 11.703 7.4324 11.027

Determinant is equal to -37

Checking the solution
3.6e-015 -7.1e-015 -1.8e-015 0

The inverse matrix of unknown variables

0.38 0.081 0.41 0.14
0.62 -0.081 1.6 0.86
0.46 0.027 1.1 0.38
0.22 0.19 0.95 0.65

Checking the inverse matrix

1 3e-017 -3e-016 0
-2e-016 1 1e-016 0
-2e-016 -5e-017 1 1e-016
2e-016 6e-017 -1e-016 1

Sum of columns

10 8 11 6

matrix norm is equal to 11

Sum of columns
2.6757 1.3784 5.0811 3.027

matrix norm is equal to 5.0811

Condition number is equal to 55.892

Результат работы программы с исходными данными:

Read matrix from the file

5.35 0.82 -0.19 -0.65 0.51 -0.36 0.9 0.07
-0.36 0.6 -0.29 0.15 0.4 4.99 0.66 -0.54
-0.65 0.87 0.76 6.12 -0.49 0.15 0.72 -0.52
-0.19 0.86 5.47 0.76 0.84 -0.29 0.45 -0.61
0.51 -0.58 0.84 -0.49 4.62 0.4 0.67 0.06
0.82 6 0.86 0.87 -0.58 0.6 -0.18 0.19
0.9 -0.18 0.45 0.72 0.67 0.66 5.39 -0.55

Coefficients of unknown variables

5.35 0.82 -0.19 -0.65 0.51 -0.36 0.9
-0.36 0.6 -0.29 0.15 0.4 4.99 0.66
-0.65 0.87 0.76 6.12 -0.49 0.15 0.72
-0.19 0.86 5.47 0.76 0.84 -0.29 0.45
0.51 -0.58 0.84 -0.49 4.62 0.4 0.67
0.82 6 0.86 0.87 -0.58 0.6 -0.18
0.9 -0.18 0.45 0.72 0.67 0.66 5.39

Intercepts

0.07-0.54-0.52-0.61 0.06 0.19-0.55

Lower triangle matrix

5.35 0 0 0 0 0 0
-0.36 0.65518 0 0 0 0 0
-0.65 0.96963 1.185 0 0 0 0
-0.19 0.88912 5.8742 -28.573 0 0 0
0.51 -0.65817 0.55394 -3.0717 4.9087 0 0
0.82 5.8743 3.6039 -17.877 -4.7848 -34.336 0
0.9 -0.31794 0.33503 -0.78254 0.94498 3.2892 4.8823

Upper triangle matrix

1 0.15327 -0.035514 -0.1215 0.095327 -0.06729 0.16822
0 1 -0.46214 0.16219 0.6629 7.5793 1.0998
0 0 1 4.9651 -0.90362 -6.112 -0.20003
0 0 0 1 -0.19517 -1.0101 -0.023768
0 0 0 0 1 1.1624 0.27418
0 0 0 0 0 1 0.15059
0 0 0 0 0 0 1

Matrix Y
0.013084 -0.81701 0.23688 0.044536 -0.097545 -0.12973 -0.071107

Solution matrix
-0.012814 0.076385 -0.12458 -0.065606 0.0603 -0.11902 -0.071107

Determinant is equal to 97666

Checking the solution
1.4e-017 0 -1.1e-016 1.1e-016 -9.7e-017 -5.6e-017 0

The inverse matrix of unknown variables

0.21 0.028 0.028 0.019 -0.025 -0.042 -0.042
-0.042 -0.032 -0.026 -0.036 0.032 0.19 0.02
0.019 0.024 -0.02 0.2 -0.045 -0.036 -0.016
0.028 0.00095 0.18 -0.02 0.021 -0.026 -0.03
-0.025 -0.025 0.021 -0.045 0.24 0.032 -0.02
0.028 0.21 0.00095 0.024 -0.025 -0.032 -0.031
-0.042 -0.031 -0.03 -0.016 -0.02 0.02 0.2

Checking the inverse matrix

1 2e-017 -2e-017 4e-018 2e-017 -5e-017 5e-018
9e-019 1 -4e-018 -8e-018 2e-019 3e-018 1e-017
-1e-017 -8e-018 1 2e-017 -2e-017 3e-017 3e-017
6e-019 -1e-017 1e-018 1 -5e-020 -1e-017 2e-017
3e-017 5e-018 9e-018 -2e-017 1 2e-018 5e-018
-4e-017 9e-018 -3e-017 -6e-018 1e-017 1 1e-017
2e-017 1e-017 6e-018 2e-017 2e-018 -1e-017 1

Sum of columns

9.78 10.91 9.86 10.76 9.11 8.45 9.97

matrix norm is equal to 10.91

Sum of columns
1.3918 1.3544 1.3035 1.3604 1.4069 1.3748 1.3634

matrix norm is equal to 1.4069

Condition number is equal to 15.349

Вывод

В ходе работы удалось познакомиться с методом Халецкого. Метод Халецкого – достаточно быстрый метод решения СЛАУ, основанных на получении матричных разложений. Однако он предъявляет к матрице системы требования симметрии и положительной определённости