Метод квадратного корня

Размещено на http://www.allbest.ru/

Лабораторная работа

по курсу

«Численные методы»

Метод квадратного корня

Задание

Для матрицы А составить процедуры:

– разложения матрицы на множители;

– решения СЛАУ ах=b;

– нахождения определителя матрицы А;

– нахождения обратной матрицы.

методом квадратного корня.

Вывести в файл:

– все входные данные задачи;

– разложенную на множители матрицу А;

– обратную к А матрицу;

– определитель матрицы А;

– решение системы Х;

– следующие проверки:

правильность разложения (невязку),

невязку нахождения решения,

невязку обращенной матрицы.

Входные данные задачи вводятся из файла, содержащего размерность задачи, матрицу А, вектор В.

1. Описание используемого метода

Суть метода квадратного корня состоит в следующем. Пусть дана линейная система, где А=(a_ij) - симметричная квадратная матрица порядка nxn, т.е. А^T=(a_ji)=A; b=(b₁,…, b_n) - вектор правых частей системы, x=(x₁,…, x_n) - вектор-столбец неизвестных.

Расчетные формулы.

a) разложения матрицы на множители.

Рассмотрим матрицу A в виде

A = S^ТDS, (1)

где

	s11	s12	s13	…	s1n
	0	s22	s23	…	s2n
S =	0	0	s33	.	s3n
	…
	0	0	0	…	snn

Sii > 0, i = 2…n

S - верхнетреугольная матрица с положительными элементами на главной диагонали.

D = diag (d11, …, dnn); dii = ±1, i = 1..n - диагональная матрица с элементами на главной диагонали ±1.

S^Т - транспонированная к S матрица (нижнетреугольная).

Тогда, производя перемножение матриц S^Т, D и S, получим уравнения для определения элементов a_ij матрицы A и, преобразовав их, получим следующие уравнения:

 (2)

(3)

(4)

квадратный корень матрица программа

Если s_ii?0, то det(A) = det(S^Т) det(D) det(S) ? 0, и система имеет единственное решение.

Если матрица A является положительно определенной, то матрица D будет единичной.

б) решения СЛАУ Ax=b.

Если разложение получено, то, применив формулу (1), получим:

S^ТDSx=b

Положив DSx=y, сводим решение системы к последовательному решению двух систем с треугольными матрицами: S^Ty=b и DSx=y.

Из первой найдем y_i:

(5)

Затем, зная y_i, из второй системы находим x_i:

(6)

в) нахождения определителя матрицы А.

Как уже отмечалось выше, det(A) = det(S^Т) det(D) det(S) = det(D) (det(S))²

det(S) = (s_11*s_{22*… *}s_nn),

det(D) = (d_11*d_{22*… *}d_nn)

Отсюда находим:

(7)

г) нахождения обратной к А матрицы.

Пусть H - обратная к А матрица. Тогда

AH=E (8)

E - единичная матрица nЧn.

Т.к. уравнение (8) матричное, необходимо записать его в векторном виде, представив матрицы Н и Е в виде систем векторов:

E=[e1, e2,…, en], ei - вектора 1Чn (единичные орты)

H=[h1, h2,…, hn], hi - вектора 1Чn (9)

Используя систему (9), преобразуем матричное уравнение (8) в систему:

, m = (1,2,…, n)

Отсюда



, m = (1,2,…, n)

Таким же образом, как и при решении СЛАУ, положим и , из чего найдем:

, m = (1,2,…, n)

, m = (1,2,…, n)

Таким образом, будет найдена матрица Н, обратная к А.

Перемножение матриц

2. Листинг программы

 // -

#include <vcl.h>

#pragma hdrstop

#include <iostream.h>

#include <fstream.h>

#include <sstream.h>

#include <math.h>

#include <tchar.h>

#include <stdlib.h>

#include <stdio.h>

#include <winbase.h>

#include <windows.h>

 // -

#pragma argsused

const int n=4;

void FindD (float A[n] [n], float S[n] [n], float D[n] [n], int i)

{float c=0, b=0; int l;

for (l=0; l<=i-1; l++)

c=c+S[l] [i]*S[l] [i]*D[l] [l];

b= A[i] [i] - c;

if (b==0) D[i] [i]=0;

else {

if (b<0) D[i] [i]=-1;

else D[i] [i]=1;};

c=0;

}

void FindSij (float A[n] [n], float S[n] [n], float D[n] [n], int i, int j)

{float c=0; int l;

for (l=0; l<=i-1; l++)

c=c+S[l] [i]*D[l] [l]*S[l] [j];

if ((S[i] [i]==0)||(D[i] [i]==0)) {S[i] [j]==0;}

else S[i] [j]=(A[i] [j] - c)/S[i] [i]*D[i] [i];

c=0;

}

void FindSii (float A[n] [n], float S[n] [n], float D[n] [n], int i)

{float c=0; int l;

for (l=0; l<=i-1; l++) {

c=c+S[l] [i]*S[l] [i]*D[l] [l];};

S[i] [i]=sqrt (fabs(A[i] [i] - c));

c=0;

}

void Findy (float A[n] [n], float St[n] [n], float b[n], float y[n], int i)

{float c=0; int k;

for (k=0; k<=i-1; k++) {

c=c+St[i] [k]*y[k];};

y[i]=(b[i] - c)/St[i] [i];

c=0;

}

void Findx (float A[n] [n], float S[n] [n], float D[n] [n], float y[n], float x[n], int i)

{float c=0; int k;

for (k=i+1; k<n; k++) {

c=c+D[i] [i]*S[i] [k]*x[k];};

x[i]=(y[i] - c)/(D[i] [i]*S[i] [i]);

c=0;

}

void FindDet (float S[n] [n], float D[n] [n], float &DetA)

{float c=1, o=1; int k;

for (k=0; k<n; k++)

c=c*S[k] [k];

for (k=0; k<n; k++)

o=o*D[k] [k];

DetA=o*(c*c);

c=1; o=1;

}

void Umn (float A[n] [n], float x[n], float b[n], float r[n])

{float c=0;

for (int i=0; i<n; i++) {

for (int j=0; j<n; j++) {

c = c + (A[i] [j]*x[j]);

}

r[i]= c-b[i]; c=0;};

}

void Umnm (float A[n] [n], float S[n] [n], float M[n] [n])

{

for (int i=0; i<n; i++) {

for (int j=0; j<n; j++) {

M[i] [j]=0;

for (int k=0; k<n; k++)

M[i] [j]= M[i] [j]+ (A[i] [k]*S[k] [j]);

};

};

}

int main (int argc, char* argv[])

{const int n=4;

float A[n] [n];

float S[n] [n], St[n] [n], D[n] [n];

float b[n];

int i, j;

for (i=0; i<n; i++) {

for (j=0; j<n; j++) {

S[i] [j]=0; St[i] [j]=0; D[i] [j]=0;

};

};

FILE* f1, *f2;

f1=fopen («input.txt», «r»);

if((f1=fopen («input.txt», «r»)) == 0) printf («file pust»);

else {

for (i=0; i<n; i++) {

for (j=0; j<n; j++) {

fscanf (f1, «%16f», &A[i] [j]);

fscanf (f1, «\n»);};};

for (i=0; i<n; i++) {

fscanf (f1, «%16f», &b[i]);

fscanf (f1, «\n»);};

fclose(f1);

f2=fopen («output.txt», «w»);

fprintf (f2, «A= \n»);

for (i=0; i<n; i++) {

for (j=0; j<n; j++) {

fprintf (f2, «%16f», A[i] [j]);};

fprintf (f2, «\n»);};

fprintf (f2, «b= \n»);

for (j=0; j<n; j++) {

fprintf (f2, «%16f», b[j]);};

fprintf (f2, «\n»);

 // поиск элементов S, D

bool B=true;

for (i=0; i<n; i++) {

for (j=i; j<n; j++) {

if (i==j) {

FindSii (A, S, D, i);

FindD (A, S, D, i);

if ((S[i] [i]==0)||(D[i] [i]==0)) B=false;

}

else

{

FindSij (A, S, D, i, j);};};

};

if (B==false) fprintf (f2, «Nedopustimaya matrica! DetA=0»);

else {

for (i=0; i<n; i++) { // поиск элементов St

for (j=0; j<n; j++) {

St[i] [j]=S[j] [i];};

};

fprintf (f2, «S= \n»);

for (i=0; i<n; i++) {

for (j=0; j<n; j++) {

fprintf (f2, «%16f», S[i] [j]);};

fprintf (f2, «\n»);};

fprintf (f2, «D= \n»);

for (i=0; i<n; i++) {

for (j=0; j<n; j++) {

fprintf (f2, «%16f», D[i] [j]);};

fprintf (f2, «\n»);};

fprintf (f2, «St= \n»);

for (i=0; i<n; i++) {

for (j=0; j<n; j++) {

fprintf (f2, «%16f», St[i] [j]);};

fprintf (f2, «\n»);};

 // Решение СЛАУ

float y[n], x[n];

for (i=0; i<n; i++)

Findy (A, St, b, y, i);

for (i=n-1; i>=0; i-)

Findx (A, S, D, y, x, i);

fprintf (f2, «y= \n»);

for (j=0; j<n; j++) {

fprintf (f2, «%16f», y[j]);};

fprintf (f2, «\n»);

fprintf (f2, «x= \n»);

for (j=0; j<n; j++) {

fprintf (f2, «%16f», x[j]);};

fprintf (f2, «\n»);

 // Нахождение определителя

float DetA;

FindDet (S, D, DetA);

fprintf (f2, «Det A=»);

fprintf (f2, «%16f», DetA);

fprintf (f2, «\n»);

 // Нахождение обратной матрицы

float H[n] [n];

for (i=0; i<n; i++) {

for (j=0; j<n; j++)

H[i] [j]=0;};

int m;

float z[n], h[n], e[n];

for (m=0; m<n; m++) {

for (i=0; i<n; i++) e[i]=0;

e[m]=1;

for (i=0; i<n; i++)

Findy (A, St, e, z, i);

for (i=n-1; i>=0; i-)

Findx (A, S, D, z, h, i);

for (i=0; i<n; i++) {

H[i] [m]=h[i];

};};

fprintf (f2, «H= \n»);

for (i=0; i<n; i++) {

for (j=0; j<n; j++) {

fprintf (f2, «%16f», H[i] [j]);};

fprintf (f2, «\n»);};

fprintf (f2, «NEVYAZKI: \n»);

float r[n];

Umn (A, x, b, r);

fprintf (f2, «Ax-b= \n»);

for (j=0; j<n; j++) {

fprintf (f2, «%16f», r[j]);};

fprintf (f2, «\n»);

float M[n] [n], F[n] [n];

Umnm (D, S, M);

Umnm (St, M, F);

fprintf (f2, «St*D*S= \n»);

for (i=0; i<n; i++) {

for (j=0; j<n; j++) {

fprintf (f2, «%16f», F[i] [j]);};

fprintf (f2, «\n»);};

float N1 [n] [n];

Umnm (A, H, N1);

fprintf (f2, «A*H= \n»);

for (i=0; i<n; i++) {

for (j=0; j<n; j++) {

fprintf (f2, «%16f», N1 [i] [j]);};

fprintf (f2, «\n»);};

};

};

system («pause»);

return 0;

}

Размещено на Allbest.ru