LU dekompozice matice

From TKwiki

Hledání vlastních čísel čtvercové matice iterativní numerickou metodou

Download

program LUdecomposition;
 
{$APPTYPE CONSOLE}
 
uses
  SysUtils;
 
CONST
   np = 4;
   presnost = 10E-16;
   maxkroku = 1000; //max pocet iteraci
 
TYPE
   RealArrayNPbyNP = ARRAY [1..np,1..np] OF real;
   RealArrayNP = ARRAY [1..np] OF real;
   IntegerArrayNP = ARRAY [1..np] OF integer;
 
VAR
  velikost:integer;
 
 
// LU dekompozice matice A na L a U bez pivotingu
PROCEDURE ludec(VAR a: RealArrayNPbyNP; VAR l: RealArrayNPbyNP; VAR u: RealArrayNPbyNP; n: integer);
var i,j,k:integer;
    r:real;
begin
  for i:=1 to n do
    for j:=1 to n do begin
      u[i,j]:=0;
      l[i,j]:=0;
      if i=j then l[i,j]:=1;
    end;
 
  for i:=1 to n do begin
    for j:=1 to n do begin
      //hledam Uij
      if i<=j then begin
        r:=a[i,j];
        if i>1 then
          for k:=1 to i-1 do
            r:=r-(l[i,k]*u[k,j]);
        u[i,j]:=r;
      end;
 
      //hledam Lij
      if i>j then begin
        r:=a[i,j];
        if i>1 then
          for k:=1 to j-1 do
            r:=r-(l[i,k]*u[k,j]);
        l[i,j]:=r/u[j,j];
      end;
 
    end;
  end;
 
end;
 
// udela A:=B*C pro matice vel. n
PROCEDURE NasobMatice(VAR a: RealArrayNPbyNP; VAR b: RealArrayNPbyNP; VAR c: RealArrayNPbyNP; n: integer);
var i,j,k:integer;
    r:real;
begin
  for i:=1 to n do begin
    for j:=1 to n do begin
 
      r:=0;
      for k:=1 to n do
        r:=r+(b[i,k]*c[k,j]);
      a[i,j]:=r;
 
    end;
  end;
 
end;
 
//zkopiruje B do A
procedure PriradMatici(var A:RealArrayNPbyNP;B:RealArrayNPbyNP;n:integer);
var i,j:integer;
begin
  for i:=1 to n do
    for j:=1 to n do
      A[i,j]:=B[i,j];
end;
 
//vypise diagonalu
procedure VypisDiagonalu(var A:RealArrayNPbyNP;n:integer);
var i:integer;
begin
  writeln;
  writeln('Vlastni cisla:');
  for i:=1 to n do
      writeln(A[i,i]:20:16);
  writeln;
  writeln('  s presnosti ',presnost);
end;
 
// vrati velikost nejvyssiho rozdilu A[i,j] a B[i,j] pro vsechny prvky
function MaxRozdilMatic(var A:RealArrayNPbyNP;B:RealArrayNPbyNP;n:integer):real;
var i,j:integer;
    r:real;
begin
  r:=0;
  for i:=1 to n do
    for j:=1 to n do
      if abs(A[i,j]-B[i,j])>r then
        r:=abs(A[i,j]-B[i,j]);
  MaxRozdilMatic:=r;
end;
 
 
//nacte matici z txt souboru:
// 1.radek - velikost matice (integer)
// 2. - n+1. radek - n hodnot oddelenych mezerami/tabulatory
procedure NactiMatici(var A:RealArrayNPbyNP; filename:string; var n:integer);
var f:text;
    i,j:integer;
    r:real;
begin
  AssignFile(f, filename);
  Reset(f);
 
  read(f,n);
  for i:=1 to n do
    for j:=1 to n do begin
      read(f,r);
      A[i,j]:=r;
    end;
 
  CloseFile(f);
end;
 
//vypise matici se zadanym uvodnim popisem s
procedure VypisMatici(A:RealArrayNPbyNP;s:string;n:integer);
var i,j:integer;
begin
  writeln(s,':');
  for i:=1 to n do begin
    for j:=1 to n do begin
      write(A[i,j]:2:2,'     ');
    end;
    writeln('');
  end;
  writeln('');
end;
 
 
var
  rozdil: real;
  A,A2,l,u:RealArrayNPbyNP;
  krok:integer;
begin
  NactiMatici(A,'matice.txt',velikost); //ted uz mame nactenou matici
 
  krok:=0;
  repeat
    inc(krok);
 
    PriradMatici(A2,A,velikost); // zkopiruje A do A2
    ludec(a,l,u,velikost); //dekompozice
    NasobMatice(a,u,l,velikost); // slozime zpet do A
 
    rozdil:=MaxRozdilMatic(A,A2,velikost);
    writeln('krok: ',krok:5,'  max rozdil: ',rozdil:16:16); // jaky je rozdil A2 a A (konverguje?)
 
    //VypisMatici(A,'A vynasobene',velikost);
 
  until ((rozdil<presnost) or (krok>maxkroku)); // koncime kdyz matice konvergovaly nebo kdyz jsme udelali hodne kroku a nekonverguji
 
  if krok>maxkroku then writeln('Divna matice :)')
  else VypisDiagonalu(A,velikost);
 
  readln;
end.