4. Схема программы
Схема основной программы.
А b – матрица коэффициентов
и вектор правых частей размером n
Вычисление:
- определителя;
- матриц L и U.
Формирование:
- массива перестановок;
- ifsolve: true – если LU-разложение выполнено успешно,
false – в противном случае.
Ly=Prb
Ux=y
x= Pcz
Схема процедуры Decomp.
true |
|
,
,
|
|
|
false |
false |
|
true |
Схема процедуры Solve
|
Формирование вектора x=Pcz
Решение уравнения Uz=y
Решение уравнения L=yb'
Перестановка элементов вектора b
(формирование Prb)
Y=b
|
5.Текст программы
program lab1;
type rmatr=array[1..4,1..4] of real;
vector=array[1..4] of integer;
vector1=array[1..4] of real;
var
F:text;
eps,det:real;
a:rmatr;
ip,jp:vector;
b,y:vector1;
i,k,n,j,p:integer;
ifsolve:boolean;
filename,s:string[20];
procedure DECOMP (n:integer; var a:rmatr; var ip,jp:vector;var eps,det:real; var ifsolve: boolean);
label 10;
var i,j,k,km,ipm,jpm:integer;
s,norm,am:real;
begin
ifsolve:=true;
det:=1.0;
for i:=1 to n do
begin
ip[i]:=i;
jp[i]:=i
end;
norm:=0.0;
for i:=1 to n do
begin
s:=0.0;
for j:=1 to n do
s:=s+abs(a[i,j]);
if norm<s then norm:=s
end;
eps:=eps*norm;
for k:=1 to n-1 do
begin
am:=abs(a[k,k]);
ipm:=k;
jpm:=k;
for i:=k to n do
for j:=k to n do
if abs(a[i,j]) > am then
begin
am:=abs(a[i,j]);
ipm:=i;
jpm:=j
end;
if ipm<>k then
begin
det:=-det;
for j:=1 to n do
begin
am:=a[k,j];
a[k,j]:=a[ipm,j];
a[ipm,j]:=am
end;
km:=ip[k];
ip[k]:=ip[ipm];
ip[ipm]:=km
end;
if jpm<>k
then
begin
det:=-det;
for i:=1 to n do
begin
am:=a[i,k];
a[i,k]:=a[i,jpm];
a[i,jpm]:=am
end;
km:=jp[k];
jp[k]:=jp[jpm];
jp[jpm]:=km
end;
if abs(a[k,k]) >= eps then
for i:=k+1 to n do
begin
am:=a[i,k]/a[k,k];
a[i,k]:=am;
for j:=k+1 to n do
a[i,j]:=a[i,j]-am*a[k,j]
end
else
begin
ifsolve:=false;
goto 10
end;
end;
10: if (ifsolve and (abs(a[n,n]) >= eps)) then
begin
for k:=1 to n do
det:=det*a[k,k];
end
else
begin
det:=0.0;
ifsolve:=false
end
end;
Procedure SOLVE(n:integer;var a:rmatr;var ip,jp:vector;var b:vector1);
var
x,y:vector1;
i,j:integer;
s:real;
begin
for i:=1 to n do
y[i]:=b[ip[i]];
for i:=2 to n do
begin
s:=y[i];
for j:=1 to i-1 do
s:=s-a[i,j]*y[j];
y[i]:=s;
end;
x[n]:=y[n]/a[n,n];
for i:=n-1 downto 1 do
begin
s:=y[i];
for j:=i+1 to n do
s:=s-a[i,j]*x[j];
x[i]:=s/a[i,i];
end;
for j:=1 to n do
b[jp[j]]:=x[j];
x:=b;
end;
begin
writeln('1 zadat matricu vruchnuyu');
readln(p);
if p=1 then
begin
write('Vvedite razmernoct matrici:');
readln(n);
writeln('Vvedite elementi matrici:');
for i:=1 to n do
for j:=1 to n do
begin
write('a[',i,j,']=');
read(a[i,j])
end;
writeln('Vvedite vektor b');
for i:=1 to n do
readln(b[i]);
end;
for i:=1 to n do
begin
for j:=1 to n do
write(' a[',i,j,']=',a[i,j]:7:4);
writeln;
end;
write('Vvedite eps: eps= ');readln(eps);
decomp(n,a,ip,jp,eps,det,ifsolve);
if ifsolve=true then
begin
solve(n,a,ip,jp,b);
for i:=1 to n do
begin
for j:=1 to n do
write(' d[',i,j,']=',a[i,j]:7:4);
writeln;
end;
writeln('Vector rewenii ');
for i:=1 to 4 do
writeln('x= ',b[i]);
writeln('det=',det);
writeln('eps=',eps);
end else
begin
writeln;
writeln(' Error: nedopystimoe yslovie (norm>eps) ');
end;
end.