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

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

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

по курсу

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

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

Задание

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Расчетные формулы.) разложения матрицы на множители.

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

= S^ТDS,                                  (1)

где

ii > 0, i = 2…n

S - верхнетреугольная матрица с положительными элементами на главной диагонали.= diag (d11, …, dnn); dii = ±1, i = 1..n - диагональная матрица с элементами на главной диагонали ±1.^{Т -} транспонированная к 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), получим:

^Т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),(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 argsusedint n=4;FindD (float A[n] [n], float S[n] [n], float D[n] [n], int i)

{float c=0, b=0; int l;(l=0; l<=i-1; l++)=c+S[l] [i]*S[l] [i]*D[l] [l];= A[i] [i] - c;(b==0) D[i] [i]=0;{(b<0) D[i] [i]=-1;D[i] [i]=1;};=0;

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

{float c=0; int l;(l=0; l<=i-1; l++)=c+S[l] [i]*D[l] [l]*S[l] [j];((S[i] [i]==0)||(D[i] [i]==0)) {S[i] [j]==0;}S[i] [j]=(A[i] [j] - c)/S[i] [i]*D[i] [i];=0;

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

{float c=0; int l;(l=0; l<=i-1; l++) {=c+S[l] [i]*S[l] [i]*D[l] [l];};[i] [i]=sqrt (fabs(A[i] [i] - c));=0;

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

{float c=0; int k;(k=0; k<=i-1; k++) {=c+St[i] [k]*y[k];};[i]=(b[i] - c)/St[i] [i];=0;

}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;(k=i+1; k<n; k++) {=c+D[i] [i]*S[i] [k]*x[k];};[i]=(y[i] - c)/(D[i] [i]*S[i] [i]);=0;

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

{float c=1, o=1; int k;(k=0; k<n; k++)=c*S[k] [k];(k=0; k<n; k++)=o*D[k] [k];=o*(c*c);=1; o=1;

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

{float c=0;(int i=0; i<n; i++) {(int j=0; j<n; j++) {= c + (A[i] [j]*x[j]);

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

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

{(int i=0; i<n; i++) {(int j=0; j<n; j++) {[i] [j]=0;(int k=0; k<n; k++)[i] [j]= M[i] [j]+ (A[i] [k]*S[k] [j]);

};

};

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

{const int n=4;A[n] [n];S[n] [n], St[n] [n], D[n] [n];b[n];i, j;(i=0; i<n; i++) {(j=0; j<n; j++) {[i] [j]=0; St[i] [j]=0; D[i] [j]=0;

};

};* f1, *f2;=fopen («input.txt», «r»);((f1=fopen («input.txt», «r»)) == 0) printf («file pust»);{(i=0; i<n; i++) {(j=0; j<n; j++) {(f1, «%16f», &A[i] [j]);(f1, «\n»);};};(i=0; i<n; i++) {(f1, «%16f», &b[i]);(f1, «\n»);};(f1);=fopen («output.txt», «w»);(f2, «A= \n»);(i=0; i<n; i++) {(j=0; j<n; j++) {(f2, «%16f», A[i] [j]);};(f2, «\n»);};(f2, «b= \n»);(j=0; j<n; j++) {(f2, «%16f», b[j]);};

fprintf (f2, «\n»);

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

bool B=true;(i=0; i<n; i++) {(j=i; j<n; j++) {(i==j) {(A, S, D, i);(A, S, D, i);((S[i] [i]==0)||(D[i] [i]==0)) B=false;

}

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

};(B==false) fprintf (f2, «Nedopustimaya matrica! DetA=0»);{(i=0; i<n; i++) { // поиск элементов St(j=0; j<n; j++) {[i] [j]=S[j] [i];};

};(f2, «S= \n»);(i=0; i<n; i++) {(j=0; j<n; j++) {(f2, «%16f», S[i] [j]);};(f2, «\n»);};(f2, «D= \n»);(i=0; i<n; i++) {(j=0; j<n; j++) {(f2, «%16f», D[i] [j]);};(f2, «\n»);};

(f2, «St= \n»);(i=0; i<n; i++) {(j=0; j<n; j++) {(f2, «%16f», St[i] [j]);};(f2, «\n»);};

// Решение СЛАУy[n], x[n];(i=0; i<n; i++)(A, St, b, y, i);(i=n-1; i>=0; i-)(A, S, D, y, x, i);(f2, «y= \n»);(j=0; j<n; j++) {(f2, «%16f», y[j]);};(f2, «\n»);(f2, «x= \n»);(j=0; j<n; j++) {(f2, «%16f», x[j]);};

fprintf (f2, «\n»);

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

FindDet (S, D, DetA);(f2, «Det A=»);(f2, «%16f», DetA);(f2, «\n»);

// Нахождение обратной матрицыH[n] [n];

for (i=0; i<n; i++) {(j=0; j<n; j++)[i] [j]=0;};m;z[n], h[n], e[n];(m=0; m<n; m++) {(i=0; i<n; i++) e[i]=0;[m]=1;(i=0; i<n; i++)(A, St, e, z, i);(i=n-1; i>=0; i-)(A, S, D, z, h, i);(i=0; i<n; i++) {[i] [m]=h[i];

};};(f2, «H= \n»);(i=0; i<n; i++) {(j=0; j<n; j++) {(f2, «%16f», H[i] [j]);};(f2, «\n»);};(f2, «NEVYAZKI: \n»);r[n];(A, x, b, r);(f2, «Ax-b= \n»);(j=0; j<n; j++) {(f2, «%16f», r[j]);};(f2, «\n»);M[n] [n], F[n] [n];(D, S, M);(St, M, F);(f2, «St*D*S= \n»);(i=0; i<n; i++) {(j=0; j<n; j++) {(f2, «%16f», F[i] [j]);};(f2, «\n»);};N1 [n] [n];(A, H, N1);(f2, «A*H= \n»);(i=0; i<n; i++) {(j=0; j<n; j++) {(f2, «%16f», N1 [i] [j]);};(f2, «\n»);};

};(«pause»);

return 0;

}