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

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

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

по курсу

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

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

Задание

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Суть метода квадратного корня состоит в следующем. Пусть дана линейная система, где А=(aij) — симметричная квадратная матрица порядка nxn, т.е. АT=(aji)=A; b=(b1,…, bn) — вектор правых частей системы, x=(x1,…, xn) — вектор-столбец неизвестных.

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

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

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

A = SТDS, (1)

где

S =

Sii > 0, i = 2…n

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

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

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

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

(2)

(3)

(4)

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

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

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

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

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

SТDSx=b

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

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

(5)

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

(6)

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

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

det(S) = (s11*s22*… *snn),

det(D) = (d11*d22*… *dnn)

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

(7)

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

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

AH=E (8)

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

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

E=, ei-вектора 1Чn (единичные орты)

H=, hi-вектора 1Чn (9)

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

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

Отсюда

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

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

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

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

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

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

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

// —

#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 , float S , float D , int i)

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

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

c=c+S *S *D ;

b= A — c;

if (b==0) D =0;

else {

if (b<0) D =-1;

else D =1;};

c=0;

}

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

{float c=0; int l;

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

c=c+S *D *S ;

if ((S ==0)||(D ==0)) {S ==0;}

else S =(A — c)/S *D ;

c=0;

}

void FindSii (float A , float S , float D , int i)

{float c=0; int l;

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

c=c+S *S *D ;};

S =sqrt (fabs(A — c));

c=0;

}

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

{float c=0; int k;

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

c=c+St *y;};

y=(b — c)/St ;

c=0;

}

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

{float c=0; int k;

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

c=c+D *S *x;};

x=(y — c)/(D *S );

c=0;

}

void FindDet (float S , float D , float &DetA)

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

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

c=c*S ;

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

o=o*D ;

DetA=o*(c*c);

c=1; o=1;

}

void Umn (float A , float x, float b, float r)

{float c=0;

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

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

c = c + (A *x);

}

r= c-b; c=0;};

}

void Umnm (float A , float S , float M )

{

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

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

M =0;

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

M = M + (A *S );

};

};

}

int main (int argc, char* argv)

{const int n=4;

float A ;

float S , St , D ;

float b;

int i, j;

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

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

S =0; St =0; D =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 );

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

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

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

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 );};

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

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

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

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

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 ==0)||(D ==0)) B=false;

}

{

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 =S ;};

};

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

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

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

fprintf (f2, «%16f», S );};

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

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

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

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

fprintf (f2, «%16f», D );};

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

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

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

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

fprintf (f2, «%16f», St );};

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

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

float y, x;

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);};

fprintf (f2, «\n»);

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

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

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

fprintf (f2, «\n»);

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

float DetA;

FindDet (S, D, DetA);

fprintf (f2, «Det A=»);

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

fprintf (f2, «\n»);

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

float H ;

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

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

H =0;};

int m;

float z, h, e;

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

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

e=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 =h;

};};

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

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

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

fprintf (f2, «%16f», H );};

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

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

float r;

Umn (A, x, b, r);

Добавить комментарий

Закрыть меню