Algoritmické riešenie matíc

V lineárnej algebre sa sústavy vyššieho stupňa zapisujú v skrátenej forme - vo forme matíc. Vznikajú tak dvojrozmerné číselné útvary - matice, ktoré majú špecifické operácie a vlastnosti. V nasledujúcom článku sa dozviete postupy, ako algoritmicky upravovať matice do "koreňového" tvaru.


Matematický pohľad:

Na úvod predpokladajme sústavu N lineárnych rovníc o N neznámych, ktorú možno vo všeobecnosti upraviť a symbolicky zapísať do tvaru:

a11x1 + a12x2 + a13x3 + a14x4 + .... a1nxn = y1
a21x1 + a22x2 + a23x3 + a24x4 + .... a2nxn = y2
a31x1 + a32x2 + a33x3 + a34x4 + .... a3nxn = y3
   ...
an1x1 + an2x2 + an3x3 + an4x4 + .... annxn = yn
Riešením takejto rovnice sú korene x1 až xn. Ak uvažujeme N-rozmerný priestor, jednotlivé rovnice budú v tomto priestore reprezentované priamkami, ktoré sa pretnú práve v jednom bode danom súradnicami x1, x2 až xn.

Maticový zápis takejto sústavy má tvar "tabuľky" s N riadkami a N+1 stĺpcami, pričom v matici sa nachádzajú čisto koeficienty a11 až ann špecifické pre každú sústavu.

a11 a12 a13 a14 ... a1n |y1
a21 a22 a23 a24 ... a2n |y2
a31 a32 a33 a34 ... a3n |y3
   ...
an1 an2 an3 an4 ... ann |yn
Pre naše potreby uvažujme konkrétnu maticu 3x3. ktorá má zápis:
Matica:           Maska:

  5 3 1 111       1 1 1 1
  3 1 2 205       1 1 1 1
100 5 3 410       1 1 1 1
Maska matice
Každá matica má svoju "masku". Maska je vlastne určitá šablóna logických hodnôt, ktorá nesie informáciu o nulových na nenulových hodnotách "matice". Nulové hodnoty matice sú v maske reprezentované hodnotou 0 (false) a nenulové hodnoty 1 (true). Našim cieľom je dovolenými úpravami vytvoriť maticu, kde budú všetky hodnoty a11 až ann okrem hodnôť na diagonále nulové:
Matica:           Maska:

  1 0 0   1       1 0 0 1
  0 1 0   2       0 1 0 1
  0 0 1 100       0 0 1 1
Z upravenej matice je možné vyčítať korene 1*x1 = 1, 1*x2 = 2, 1*x3 = 100.

Hodnosť matice
Samozrejme, existujú matice ktoré nemajú riešenie. Ináč povedané, neexistuje N-tica koreňov, ktorá by vyhovovala všetkým N-rovniciam. V tejto súvislosti treba spomenúť tazvanú hodnosť matice. Hodnosťou matice sa rozumie počet lineárne nezávislých riadkov. Riadok je lineárne závislí vtedy, ak je lineárnou kombináciou ľubovolných dvoch ďalších riadkov. Ak v priebehu "úprav" dostaneme riadok so samými nulami, okamžite ho môžeme z matice vyhodiť von, pričom hodnota matice sa zníži o 1.

Frobeniova veta
Frobeniova veta hovorí, že matica je riešiteľná iba vtedy, ak hodnota matice (a11 až ann) je rovnaká ako hodnota rozšírenej matice (a11 až yn).

Úpravy a ich algoritmická implementácia:

V našom programe bude maticu preprezentovať dvojrozmerné pole reálnych hodnôt. Pričom predpokladáme že maximálny rozmer matice nepresiahne hodnotu 100.
     const short int max = 100;
     float[max][max] matrix;
Násobenie riadka konštantou
Násobenie riadka konštantou považujeme za elementárnu operáciu, pri ktorej sa každý riadok matice vynásobí nenulovou reálnou konštantou. V našej matici reprezentovanej poľom matrix bude riadok N určený hodnotami matrix[1][n]matrix[n][n+1].

void nasobenie (int riadok, float konstanta) {
     for (int i = 1; i < = n + 1; i++)
         matrix[i][riadok] *= konstanta;
}
Výmena riadkov
Výmena riadkov je triviálna operácia, ktorá má význam iba pri ručnom počítaní matice (z dôvodu prehľadnosti zápisu). Pri algoritmických výpočtoch stráca význam.

Lineárna kombinácia
Lineárna kombinácia je operácia, pri ktorej riadok nahradíme lineárnou kombináciou tohto riadku s iným riadkom. Lineárna kombinácia riadkov spočíva v sčítaní ich hodnôt. Teda, kombinácia riadkov 1 a 2 bude mať tvar ( a11 + a12 ; a21 + a22 ; .... ; an1 + an2 )
Vstupom nasledujúcej funkcie sú čísla "zdrojového" a "cieľového" riadka v poli, pričom funkcia "cieľový" riadok nahradí lineárnou kombináciou oboch riadkov.
void kombinacia (int zdroj, int ciel) {
     for (int i = 1; i < = n + 1; i++)
         matrix[i][ciel] *= matrix[i][zdroj];
}
Postup riešenia:
Riešenie matice spočíva v kombinácií úprav tak, aby bol výsledný tvar "diagonálny". Vhodnými lineárnymi kombináciami vytvárame postupne nuly na pozíciách mimo diagonály.
 -- 12 11 10  --
  1 --  9  8  --
  2  3 --  7  --
  4  5  6 --  --
Vhodnými úpravami sa najskôr snažíme dostať nulu na pozícií označenej jednotkou matrix[1][2]. Potom postupne vytvárame nuly na pozíciách 2, 3, 4 ... 12.
Vytvorme si teda cyklus, ktorý postupne "prebehne" všetky indexy prvkov, ktoré treba upraviť na nulový tvar:
   for (int y=2; y<=n; y++)
      for (int x=1; x<=y-1; x++) {
        // vytvoriť nulu pre matrix[x][y]

   int currline = 0;
   for (int y=n-1; y>=1; y--) {
      currline++;
      for (int x=n; x>=(n-currline+1); x--) {
        // vytvoriť nulu pre matrix[x][y]
   }
Vytvorenie nulového prvku na pozícií matrix[x][y] spočíva vo vhodnej lineárnej kombinácií tohto riadka (y) s iným riadkom tak. Predpokladajme, že riadok vhodný na kombináciu riadkom y má číslo z. Potom musí platiť

matrix[x][y] == - matrix[x][z]

Napíšme si teda funkciu, ktorá upraví riadky y a z tak, aby platila horeuvedená podmienka:
void uprav(int x, int y, int z) {
     float prve = matrix[x][y];
     float druhe = matrix[x][z];
     for (int i=1; i<=n+1; i++) {
         matrix[i][y] *=  fabs(druhe);
         matrix[i][z] *= -fabs(prve);
     }
}
Príklad: riadky y a z upravujeme pomocou funkcie uprav
riadok y:  3 -1  4  8  9
riadok z: -2 -9  1  2  3
Násobenie riadkov na najmenší spoločný násobok:
       y:  6 -2  8 16 18    / * abs(-2)
       z: -6 27  3  6  9    / * -abs(3)
Keď sú riadky upravené jeden z nich (y) nahradíme lineárnou kombináciou oboch riadkov. Na tento učel si vytvorme funkciu combine, ktorá vykoná samotnú lineárnu kombináciu:
void combine (int y, int z) {
     for (int i=1; i<=n+1; i++)
         matrix[i][z] += matrix[i][y];
}
Záver:
Ak funkcie na úpravu a kombináciu matíc vložíme do vnoreného cyklu, ktorý postupne prebehne vhodné indexy poľa, na konci cyklu dostávame maticu upravenú do diagonálneho tvaru, ktorú možno považovať za riešenie sústavy lineárnych rovníc.

Na záver prikladám hotový program v c++, ktorý načíta vstupnú maticu zo súboru ./matrix.txt a všetky úpravy, ako aj výslednú maticu vypíše na štandardný výstup. Keďže program rieši lineárne sústavy N x N, je nevyhnutné aby vstup obsahoval maticu s N riadkami a N+1 stĺpcami. Hodnoty môžu byť oddelené medzerou alebo \n, prípadne \t.

Príklad vstupu:
  5 3 1 111
  3 1 2 205
100 5 3 410
Príklad výstupu:
1.0 0.0 0.0 1.0
0.0 1.0 0.0 2.0
0.0 0.0 1.0 100.0
Zdrojový súbor v cpp:
#include <iostream.h>
#include <stdlib.h>
#include <fstream.h>
#include <math.h>
#include <stdio.h>
#include <string>
#include <cassert>
#define DEBUG(x) cout<<"Debug:"<<#x<<"="<<x<<endl;
#define TEMPLAT templat_out();
#define MATRIX matrix_out();
const char* source = "./matrix.txt";  // vstup
const int max = 100;   // max rozmery  matice
const char* format_string = "%3.1f "; // formátovanie výstupu
float* create_datastream();
bool check_stream(int);
void templat_out();
int nn1(float);
void create_matrix(float*);
void matrix_out();
int find0 (int, int);
int find0b (int, int);
int linear(int, int, int);
void make_x_positive(int, int);
void negative_row(int);
void combine(int, int, int);
void saveline (int);
void restoreline (int);
int matrix_value(int, int);
bool isnull(float);

using namespace std;
ifstream sourcefile;
float matrix[max][max];
bool templ[max][max];
float templine[max];
float* datastream = create_datastream();
int n;

bool is_null(float num) {
     bool flag;
     fabs(num) <= 0.00000001 ?
               flag = true : flag = false;
     return flag;
}

float* create_datastream () {
     sourcefile.open(source);
     float* stream = new float;
     float x; int index = 0;
     while (sourcefile >> x)
           stream[++index] = x;
     stream[0] = index;
     if (!check_stream(index)) {
        cout << "\nZly vstup";
        assert(false);
     }
     return stream;
}

bool check_stream(int p) {
    bool flag = false;
    (float)((-1+sqrt(1+4*p)))/2 ==
    floor ((-1+sqrt(1+4*p)))/2 ?
                 flag = true : flag = false;
    return flag;
}

int nn1(float p) {
    int out = (int)floor((-1+sqrt(1+4*p))/2);
    int n; out >= 0 ? n = out : n = - out;
    return n;
}

void create_matrix(float* datastream) {
     n = nn1(datastream[0]);
     int i = 0;
     for (int y = 1; y <=n; y++)
        for (int x = 1; x <= n+1; x++) {
            matrix[x][y] = datastream[++i];
            templ[x][y] = true;
        }
}

void matrix_out() {
     n = nn1(datastream[0]);
     for (int y = 1; y <=n; y++) {
        for (int x = 1; x <= n+1; x++)
            printf (format_string, matrix[x][y]);
        cout << endl;
     }
     cout << endl;
}

void templat_out() {
     n = nn1(datastream[0]);
     for (int y = 1; y <=n; y++) {
        for (int x = 1; x <= n+1; x++)
            cout << templ[x][y] << " ";
        cout << endl;
     }
     cout << endl;
}

int find0 (int x,int y) {
    for (;y>=2;)
        if (templ[x][--y]==true) break;
    return y;
}

int find0b (int x,int y) {
    for (;y<=n-1;)
        if (templ[x][++y]==true) break;
    return y;
}
void make_x_positive(int x, int y) {
     if (matrix[x][y]<0)
        for (int i=1; i<=n+1; i++)
            matrix[i][y]*=-1;
}

void negative_row(int y) {
     for (int i=1; i<=n+1; i++) matrix[i][y]*=-1;
}

void combine (int dy, int sy, int x) {
     for (int i=1; i<=n+1; i++)
         matrix[i][dy]+=matrix[i][sy];
     matrix[x][dy] == 0 ? templ[x][dy] = false : 0;
}

int matrix_value(int xm, int ym) {
    return 0;
}

int linear(int dy, int sy, int x) {
    make_x_positive(x,dy);
    make_x_positive(x,sy);
    bool do_it = false;
    if (!is_null(matrix[x][dy])) do_it = true;

    if (do_it) {
       float sy_origin = matrix[x][sy];
       float dy_origin = matrix[x][dy];

       saveline(sy);
       for (int i=1; i<=n+1; i++) {
           matrix[i][sy]*=dy_origin;
           matrix[i][dy]*=sy_origin;
       }
       negative_row(dy);
       combine (dy, sy, x);
       restoreline(sy);
    }
    templ[x][dy] = false;
}

void saveline(int sy) {
    for (int i=1; i<=n+1; i++)
        templine[i] = matrix[i][sy];
}
void restoreline(int sy) {
    for (int i=1; i<=n+1; i++)
        matrix[i][sy] = templine[i];
}

int main() {
   create_matrix(datastream);
   MATRIX;

   for (int iy=2; iy<=n; iy++)
      for (int ix=1; ix<=iy-1; ix++) {
        linear(iy, find0(ix,iy),ix);
        cout << "Matrix content:" << endl;
        MATRIX;
        cout << "Matrix mask:" << endl;
        TEMPLAT;
      }

   int currline = 0;
   for (int iy=n-1; iy>=1; iy--) {
      currline++;
      for (int ix=n; ix>=(n-currline+1); ix--) {
        linear(iy, find0b(ix,iy),ix);
        cout << "Matrix content:" << endl;
        MATRIX;
        cout << "Matrix mask:" << endl;
        TEMPLAT;
      }
   }
   cout << "---- SOLUTION ----\n";
   for (int i=1; i<=n; i++) {
       matrix[n+1][i] /= matrix[i][i];
       matrix[i][i] = 1;
   }
   MATRIX;
   return 0;
}


Dev, noemail@noemail.com
Stiahnuté z Developer.sk