MFF UK / Ústav teoretické fyziky / Tomáš Ledvinka
Přednášky
. . . . . . . . . . . . . . . . . . . . . . . . . .
Programování pro fyziky (1.r)
  Úvod
  Programy a algoritmy
  Píšeme program
  Píšeme program II
  Procedury a funkce
  Malujeme funkce
  Chyby. Typy I.
  Typy II. Pole a Záznamy
  Pole II.Řetězce.Soubory.
  Gnuplot.Interpolace...
  Matice. Velké O...
    Matice
    Soustava lineárních rovnic
    Časová náročnost algoritmu
    Seznamy
  Fronta,Zásobník. Postscript
  Bin. soubory, ...
  Ukazatele,Objekty, ...
Počítačová algebra
Klasická elektrodynamika (2.r)
Vybrané partie OTR

Cvičení
. . . . . . . . . . . . . . . . . . . . . . . . . .
Programování pro fyziky (1.r)
Teoretická mechanika (2.r)
Klasická elektrodynamika (2.r)


Věda
. . . . . . . . . . . . . . . . . . . . . . . . . .
Diskové zdroje v OTR
Hyperbolické systémy v OTR


Kontakt
. . . . . . . . . . . . . . . . . . . . . . . . . .
Email
Konzultační hodiny


Ostatní
. . . . . . . . . . . . . . . . . . . . . . . . .
Mallorca
Ze společnosti

Řešení soustavy lineárních rovnic (Gaussova - Jordanova eliminace)

Mějme soustavu 4 rovnic pro 4 neznámé:

     9y + 6z + 4w = 3
2x + 3y + 4z      = 3
 x      - 3z + 2w = 0
 x + 4y + 4z -  w = 1

Tato soustava se běžně representuje maticí

Protože matici budme chtít upravit na diagonální tvar, kde nenulové prvky budou jen na diagonále, musíme nejdříve prohodit první dva řádky (tedy první dvě rovnice, to jistě smíme):

Nyní první řádek (tedy rovnici) vydělíme tak aby na diagonále zbyla jednička - normalizace

Nyní od druhého až čtvrtého řádku odečteme vhodný násobek prvního řádku tak aby v prvním sloupci mimo diagonálu zbyly nuly (lineární kombinace rovnic je OK):

Teď již stručně - normalizace 2. řádku (nic nemusíme prohazovat)

odečtení násobků druhého řádku

Normalizace 3.řádku

Odečtení násobků 3. řádku od ostatních

Normalizace 4. řádku

Konečně, odečtení násobků čtvrtého řádku


Připoměňme, že tato matice representuje lineární systém
x = 3/4
y = 23/2
z = -33/4
w = -51/4

Celá metoda tak postupně opakuje tři kroky: podmíněné prohození řádků, normalizaci řádku a eliminaci nediagonálních elementů daného sloupce.
Technický detail: je zvykem prohazovat řádky tak, aby se na diagolále objevil v absolutní hodnotě nejjvyšší použitelný (řádek>=sloupec) člen daného sloupce. To že nestačí jen kontola na 0 tak aby šel řádek normalizovat je zřejmé, i číslo 10^-15 by nám vadilo (vyzkoušejte). To že se vyplatí najít právě v absolutní hodnodě největší prvek je už numerická matematika.


Následuje program, ve kterém je také procedura pro GJ eliminaci. Současným výpočetm pro N pravých stran se spočte inverzní matice (taková aby M M^(-1) = JednotkovaMatice).

Program Maticni;
uses Sysutils;
{Pouzivame assert a ten v Delphi-IDE neni videt bez SysUtils}

type tVektor = array of real;
     tMatice = array of tVektor;

function _v(const a : array of real) : tVektor;
{Udela z pole realnych cisel vektor}
var b : tVektor;
    i : integer;
begin
    SetLength(b,High(a)+1);
    for i :=0 to High(a) do b[i] := a[i];
    _v:=b;
end;

function _m(const a : array of tVektor) : tMatice;
{udela z pole vektoru matici}
var b : tmatice;
    i : integer;
begin
    SetLength(b,High(a)+1);
    for i :=0 to High(a) do begin
      Assert(High(a[i])=High(A[0]),'_m: Matice musi byt obdelnikova!');
      b[i] := a[i];
    end;
    _m:=b;
end;

Procedure WriteM(const M : tMatice;W:integer=9;D:integer=5);
var r,s : integer;
begin
  for r:=0 to High(M) do begin
    for s:=0 to High(M[r]) do Write(M[r,s]:W:D);
    Writeln;
  end;
end;

Function IdentM(N:integer):tMatice;
var Q  : tMatice;
    r,s: integer;
begin
  SetLength(Q,N,N);
  for r:=0 to N-1 do
    for s:=0 to N-1 do
      if r=s then Q[r,s]:=1 else Q[r,s]:=0;
  IdentM:=Q;
end;

Function TranspM(const A:tMatice):tMatice;
var Q  : tMatice;
    N,M,r,s: integer;
begin
  N := High(A);
  M := High(A[0]);
  SetLength(Q,M+1,N+1);
  for r:=0 to M do
    for s:=0 to N do
      Q[r,s]:=A[s,r];
  TranspM:=Q;
end;


function GJElim(const A,B: tMatice): tMatice;
{Samozrejme resi A.x_i=b_i pro vice pravych stran}
{Cviceni: zkuste funkci pretizit a pridat dalsi pro jedinou pravou stranu}
var M,Z : tMatice;
    t : tVektor;
    r1,r,s,N,L,rmax : integer;
    u : real;

begin
    N := High(A);  // rozmery ctvercove matice A
    L := High(A)+High(B)+1;

    Assert(N = High(A[0]),'GJElim: Matice musi byt ctvercova!');
    Assert(N = High(B[0]),'GJElim: Vektory v B musi byt spravne dlouhe!' );

    {1. A a B spojime do jedine matice}
    SetLength(M, N+1{x},L+1);

    for r := 0 to N do {radek}
      for s := 0 to N do  {sloupec}
        M[r,s] := A[r,s];
    for r := 0 to N do
      for s := 0 to High(B) do
        M[r,s+N+1] := B[s,r];

    {2. Gauss-Jordanova eliminace}
    for r := 0 to N do begin{pro vsechny radky}
      {2.1 Najdi nejvetsi |prvek| v r-tem slopecku na radcich r..N}
      rmax := r;
      s := r;     {v r-tem sloupecku, ze}
      for r1 := r+1 to N do if abs(M[r1,s])>abs(M[rmax,s]) then rmax:=r1;
      {2.2 Prehod r-ty a rmax-ty radek}
      if r<>rmax then begin
         t := M[r];
         M[r] := M[rmax];
         M[rmax] := t;
      end;
      {2.3 Normalizuj radek r}
      u := M[r,r];
      Assert(u<>0,'GJElim: Matice neni regularni');
      for s := 0 to L do begin
         M[r,s] := M[r,s]/u;
      end;

      {2.4 odecti od ostatnich radku r-ty*A[r,s] }
      for r1:=0 to N do if r<>r1 then begin
         u := M[r1,r];
         M[r1,r]:=0;
         for s := r+1 to L do
            M[r1,s]:=M[r1,s]-u*M[r,s];
      end;
    end;

    {3. Vratit vysledek}
    SetLength(Z,High(B)+1,N+1);
    for r := 0 to N do
      for s := 0 to High(B) do
        Z[s,r]:=M[r,s+N+1];

    GJElim := Z;
end;

Function InvM(const M:tMatice):tMatice;
begin
   InvM:=TranspM(GJElim(M,IdentM(High(M)+1)))
end;

var M   : tMatice;


begin

   M := _m([_v([ 0, 4, 0]),
            _v([ 2, 0, 0]),
            _v([ 0, 8, 8])] );

   WriteM(M);
   Writeln;
   WriteM(InvM(M));

   Readln;
end.

{ A takto si můžeme vyzkoušet, zda funguje transpozice 
   WriteM(_m([_v([ 11, 12, 13]),
              _v([ 21, 22, 23])] ) );
   Writeln;
   WriteM(TranspM( _m([_v([ 11, 12, 13]),
                       _v([ 21, 22, 23])] ) ));
   Writeln;
}

 

.