Форум на Исходниках.RU

unit Gauss;

interface

type

Matrix_6x6=

array[1..6] of array[1..6] of real; //Матрица 6x6

VectorStolbec6=

array[1..6] of real; //Вектор-столбец свободных членов

function GaussMethod(Matrix:Matrix_6x6; Vector:VectorStolbec6):VectorStolbec6;

implementation

function GaussMethod(Matrix:Matrix_6x6; Vector:VectorStolbec6):VectorStolbec6;

var LowMatrix:Matrix_6x6;

UpMatrix:Matrix_6x6;

i,j,k,l:integer;

TempSum,TempSum2:real;

Y, X:VectorStolbec6; //Векторы Y, X для решения 2 систем уравнений

begin

{!-------------Часть 1. Разложение A=LU-------------!}

//Обнуление матриц и векторов

for i:=1 to 6 do

for j:=1 to 6 do

begin

UpMatrix[i,j]:=0;

LowMatrix[i, j]:=0;

Y[i]:=0;

X[i]:=0;

end;

//Самые примитивные вычисления

for i:=1 to 6 do

LowMatrix[i, 1]:=Matrix[i, 1];

for i:=2 to 6 do

UpMatrix[1, i]:=Matrix[1, i]/LowMatrix[1,1];

for i:=1 to 6 do

UpMatrix[i,i]:=1;

//Начало разложения:

TempSum:=0;

TempSum2:=0;

for i:=2 to 6 do

for j:=2 to 6 do

begin

//Верхняя треугольная матрица U

if (j>i) and (j>1) then

begin

TempSum2:=0;

for l:=1 to i-1 do

TempSum2:=TempSum2+LowMatrix[i,l]*UpMatrix[l,j];

Form1.Caption:=inttostr(i);

UpMatrix[i,j]:=(Matrix[i,j]-TempSum2)/LowMatrix[i,i];

end;

//Нижняя треугольная матрица L

if (j>1) and (j<=i) then

begin

TempSum:=0;

for k:=1 to j-1 do

TempSum:=TempSum+LowMatrix[i,k]*UpMatrix[k, j];

LowMatrix[i,j]:=Matrix[i,j]-TempSum;

end;

{!-------------Конец разложения-------------!}

{!-------------Часть 2. Нахождение корней--------!}

//Сначала решаем систему Ly=b и находим вектор Y

TempSum:=0;

for i:=1 to 6 do

begin

if i>=1 then

begin

TempSum:=0;

for k:=1 to i do

TempSum:=TempSum+LowMatrix[i, k]*Y[k];

Y[i]:=(Vector[i]-TempSum)/LowMatrix[i,i];

end;

{Нужно было для отладки

Form1.Memo1.Lines.Add('Y1: '+floattostr(Y[1]));

Form1.Memo1.Lines.Add('Y2: '+floattostr(Y[2]));

Form1.Memo1.Lines.Add('Y3: '+floattostr(Y[3]));

//----конец

}

//После нахождения вектора Y переходим к нахождению вектора X

//Его находим из системы Ux=y

TempSum:=0;

for i:=6 downto 1 do//Так как ход обратный

begin

TempSum:=0;

for k:=i+1 to 6 do begin

TempSum:=TempSum+UpMatrix[i,k]*X[k];

end;

X[i]:=Y[i]-TempSum;

end;

{!-------------Конец. Корни найдены, возвращаем результат в виде вектора--------!}

Result:=X;

end;

end.