Лабораторные работы по численным методам / 1 / 1

uses crt, Arrays;
const n=40;
var i,j,k,l: integer;
aa,a,aCopy,edn,check: array [1..n,1..n] of real;
bbb,bb,b,ab,c,d,dd: array[1..n] of real;
prog,z,det,obdet,max,maxo,maxn,sum: real;
f,f1,f2: text;

begin

var start:=System.DateTime.Now;

assign(f,'исходная матрица.txt');
Reset(f);
assign(f1,'обратная матрица.txt');
rewrite(f1);
assign(f2,'корни уравнений.txt');
rewrite(f2);

//считывание исходной матрицы
for i:=1 to n do
begin
for j:=1 to n do
read(f,a[i,j]);
readln(f,b[i]);
end;

//находим норму для исходной матрицы
max:=0;
for i:=1 to n do
begin
sum:=0;
for j:=1 to n do sum:=sum+abs(a[i,j]);
d[i]:=sum;
if d[i]>max then max:=d[i];
end;

bb:=b;
aa:=a;
aCopy:=a;

for i:=1 to n do
edn[i,i]:=1;
for i:=1 to n do
begin
c[i]:= a[i,i];
b[i]:=b[i]/c[i];
for j:=1 to n do
begin
edn[i,j]:=edn[i,j]/c[i];
a[i,j]:=a[i,j]/c[i];
end;
for l:=i+1 to n do
begin
z:=a[l,i];
for k:=1 to n do
begin
a[l,k]:=a[l,k]-a[i,k]*z;
edn[l,k]:=edn[l,k]-edn[i,k]*z;
end;
b[l]:=b[l]-b[i]*z;
end;
end;

for i:=n downto 1 do
begin
for l:=i-1 downto 1 do
begin
z:=a[l,i];
for k:=1 to n do
begin
a[l,k]:=a[l,k]-a[i,k]*z;
edn[l,k]:=edn[l,k]-edn[i,k]*z;
end;
b[l]:=b[l]-b[i]*z;
end;
end;


for i:=1 to n do
for j:=1 to n do
begin
check[i,j]:=0;
for l:= 1 to n do
check[i,j]:=check[i,j]+aCopy[i,l]*edn[l,j];
end;
for i:=1 to n do
begin
for j:=1 to n do
write(f1,edn[i,j]:3:5,' ');
writeln(f1);
end;

//находим норму для обратной матрицы
maxo:=0;
for i:=1 to n do
begin
sum:=0;
for j:=1 to n do sum:=sum+abs(edn[i,j]);
dd[i]:=sum;
if dd[i]>maxo then maxo:=dd[i];
end;

//вывод корней в файл
for i:=1 to n do
writeln(f2,'x[',i,'] = ',b[i]:8:5);

//определитель
obdet:=1;
for i:=1 to n do
obdet:=obdet*c[i];
for k:=1 to n do
for i:=k+1 to n do
for j:=k+1 to n do
begin
edn[i,j]:=edn[i,j]-edn[i,k]*edn[k,j]/edn[k,k];
det:=edn[1,1];
end;

//исходную матрицу умножаем на корни уравнений
for i:=1 to n do
begin
ab[i]:=0;
for j:=1 to n do
ab[i]:=ab[i]+aa[i,j]*b[j];
end;

// вычитаем свободные члены(невязка)
for j:=1 to n do
bbb[j]:=ab[j]-bb[j];

//находим норму невязки
maxn:=abs(bbb[1]);
for j:=1 to n do
if abs(bbb[j])>maxn then maxn:=abs(bbb[j]);

// норму обратной матрицы * на норму невязку
prog:=maxn*maxo;
writeln('Прогрешность:',prog);

for i:=2 to n do
det:=det*edn[i,i];
writeln('Определитель исходной матрицы = ',obdet);
writeln('Определитель обратной матрицы = ',det);
writeln('Проверка.Произведение определителей = ',obdet*det);
writeln('Норма исходной матрицы =',max);
writeln('Норма обратной матрицы =',maxo:7:5);
writeln('Число обусловленности = ',(max*maxo):7:5);

close(f);
close(f1);
close(f2);

var finish:=System.DateTime.Now;
writeln('Время выполнения = ',(finish - start).TotalSeconds,' секунд');

end.