http://www.math.hut.fi/teaching/v/3/L/L1.html
Luennot 1-2

Matriisit ja lineaariset yhtälösysteemit,Maple ja Matlab-perusteita

Laode
Kreyszig : 7.4. s. 344 ->
Strang
Simo Kivelä Algebra ja Geometria luku 4
1. Matriisialgebraa

Kertaa KRE 7.1 - 7.3 tai GRE 10.1 - 10.3

  
    Matriisitulo:  A (m x n), B (n x p), C=AB (m x p)   

                       c[i,j]=sum(a[i,k]*b[k,j],k=1..n) 
   
 Matriisitulo ei yleisesti kommutoi. (A B # B A)               
   
 Transpoosi A^T
   
 Yksikkömatriisi Id:=diag(1,1,...1);
   
 Käänteismatriisi A^-1	
   
 Determinantit ??


 Matriisien käsittely Maple:lla 

Kts. [HAM] s. 147 ->
	  
> with(linalg):  # Ladataan lineaarialgebrapakkaus
> A:=matrix(3,4,[a[1,1],a[1,2],...,a[3,4]]); 
          3 x 4 - matriisi täytetään vaakariveittäin.
      
> B:=matrix(4,5,b); # 
Toinen tapa on antaa vaakalistojen listana (lähellä Matlab-syntaksia).
matrix([[1,2,3],[4,5,6],[7,8,9]]);

Matriisialgebra toimii + ja - laskuissa normaaliin tapaan. Skalaarilla
kertomisoperaatio on *

Matriisikertolaskun symboli on &*
Normaali tapa laskea on
> C:=evalm(A &* B);

Maplessa tarvitaan aina loppumerkki (; tai :), kaksoispiste (:) estää
tulostuksen.

Maple käsittelee sekä symbolisia että numeerisia matriiseja.
Rationaaliluvuilla se operoi tarkasti, ellei erityisesti pyydetä
liukulukulaskentaa.


 Matriisien käsittely Matlab:lla 


>> A=[1,2,3,4;5,6,7,8;9,10,11,12]  % Pilkun sijasta käy erottimena väli
>> B=A'
>> C=A*B


Matlab:ssa ei tarvita loppumerkkiä. Jos tulostus halutaan estää, lopetetaan
puolipisteeseen (;)
Matlab käsittelee vain numeerisia matriiseja, mutta niitä sitten tehokkaasti.

Symbolisessa muodossa operointi ei ole mahdollista (symbolic toolbox
laajentaa tähän suuntaan, mutta sitä emme tässä käsittele).
Ratinaaliaritmetiikkaa voidaan käyttää komennolla 

 format rational  



 2. Lineeriset yhtälösysteemit, Gaussin eliminaatio 

   A x =b, missä
     A on kerroinmatriisi (m x n) (m yhtälöä, n tuntematonta)
     x on tuntemattomien muodostama sarakevektori     (n x 1)
     b on annettujen lukujen muodostama sarakevektori (m x 1)



Ratkaisujen lukumäärä voi olla 0, 1 tai ääretön, kuten jo 2 x 2 systeemeistä
tiedetään (yhdensuunt. suorat, leikkaavat suorat, yhtyvät suorat).

 Riviajattelu: Suorien leikkauspiste  

 Sarakeajattelu: Etsitään annetun vektorin (b) koordinaatit
matriisin sarakevektorien avulla esitettynä.

 Gaussin eliminaatio, s. 347

 Alkeisrivioperaatiot 

   1. rivi[i] <-->  rivi[k]  (rivien vaihto)
   2. rivi[i] <--  c*rivi[i] , c ei ole 0
   3. rivi[i] <--  rivi[i]+c*rivi[k] , c ei ole 0


Rivioperaatioissa yhtälösysteemi säilyy ekvivalenttina (samat ratkaisut).
Tähän perustuu ratkaisu Gaussin menettelyllä.
Muunnetaan systeemi yläkolmiomuotoon. Rivioperaatiot voidaan helpoimmin
kohdistaa suoraan matriisiin, jossa A:n perään liitetään b-sarake. Tätä
sanotaan liitännäismatriisiksi, "augmented matrix" 
Matlab 
Rivioperaatiot:
A([i k],:)=A([k i],:)
A(i,:)=c*A(i,:)
A(i,:)=A(i,:)+c*A(i,:)
% Testimatriisi (laode):
e2_3_8



 Maple 

> with(linalg):
> A:=matrix(m,n,[a[1,1],...,a[m,n]]);
> b:=vector([b[1],...,b[m]);
> Ab:=augment(A,b);

 Matlab 

>> A=[1 2 3;4 5 6; 7 8 9 ];b=[1;1;1];Ab=[A b]



Luennolla katsotaan esimerkkejä, kts. myös
Esimerkki rivioperaatioista Maplella tehtynä .

 Porrasmuoto, Echelon form 

Gaussin algoritmi johtaa aina muotoon:

   a[1,1]  a[1,2]       ...           a[1,n]   b[1]
     0     c[2,2]       ...           c[2,n]   bmato[2]
     0       0
     0       0
     
     

     0       0   0.. 0  k[r,r] ...      k[r,n]  bmato[r]
     0       0    ....    0    ...        0     bmato[r+1]

     0       0    ....    0    ...        0     bmato[m]



Tässä "päälävistäjän" alkiot  a[1,1], c[2,2], ..., k[r,r] 
ovat nollasta poikkeavia.

Tämä muoto saadaan aikaan sarakevaihtojen avulla siten, että aina kun
käsiteltävä tukialkio (pivot) = 0 ja sen alapuolella oleva sarake on
niinikään pelkkää 0:aa, niin ko. sarake siirretään loppuun (siis A-
osan loppuun).

Edellä saadusta muodosta on helppo lukea kaikki mahdolliset tapaukset.
Itse ratkaisun kannalta sarakevaihdot aiheuttavat ylimääräisen komplikaation,
sillä tuntemattomien järjestyksestä täytyy silloin pitää kirjaa.
Niinpä käytännössä niitä ei yleensä tehdä, vaan annetaan päälävistäjän
porrastua tarpeen mukaan pitempiin askelmiin.

Olennaista on, että siirryttäessä seuraavalle riville alaspäin,
1. nollasta poikkeva alkio siirtyy oikealle ainakin 1. askeleen.

Ratkaisujen olemassaolo ja lukumääräkysymykset ovat kaikki luettavissa
tästä kaaviosta. Samme peruslauseen (esitämme "kauniimman" version 
L3:ssa:


Lause (linsys) 

 (a) r < m ja jokin bmato[i] # 0, r+1 <= i <= m
     Ei ratk.  
Esimerkki  

(b) r=n ja bmato[i]=0, i=r+1 .. m (tai nämä (i > r)bmadot puuttuvat (tap.  m=n))
      yksikäs. ratk. 
Esimerkki 

 (c) r < n ja bmato[i]=0, i=r+1 .. m (tai "kriittiset" bmadot puuttuvat
                                      (tap. m < n ja r = m) )
      äärettömän monta ratkaisua 
Esimerkki  (Myös (a)-kohdassa oli jo yksi.)



 Havainnollistus 

 Tapaus a) r < m  ja B2 # 0  

  a)   r < m  ja B2 # 0
                       n                          n
      ------------------               ----------
     |                 |  |            |        |  |
     |                 |  | B1         |        |  | B1                  |  |            |        |
     |                 |  |            |        |  |
  r   ------------------             r ----------
     |  0 ...       0  |  | B2         | 0..0   |  | B2
     |  0 ...       0  |  |            | 0..0   |  | 
  m   ------------------             m ----------


Saadaan ristiriitainen yhtälö 0=B2(i) jollain i.
 Ei ratkaisuja 

 Tapaus b) r = n  ja (B2 = 0 tai puuttuu)  

(B2 puuttuu <==> r=m, joten tässä puutostilanteessa r=n=m)
                       n       
     -------------------
     |                 |  |
     |                 |  |            
     |                 |  | B1         
     |                 |  |            
     |                 |  |
     |                 |  |
  r   ------------------             
     |  0 ...       0  |  0 B2    Tämä siis voi puuttua,    
     |  0 ...       0  |  0       jolloin m=r=n     
  m   ------------------             

 Yksikäsitteinen ratkaisu.  Tällöin porrasmuoto
on joka kohdassa "yksiaskelinen". Takaisinsijoitusvaiheessa
jokaisessa yhtälössä yksi uusi tuntematon. Siis todellakin
yksikäs. ratkaisu.

 Tapaus c) r < n  ja (B2 = 0 tai puuttuu)  

(B2 puuttuu <==> r=m)

                     r    n       
     ----------------------
     |               |    |  |            
     |               |    |  | B1         
     |    E_{r x r}      | P  |  |           E_{r x r}   
     |               |    |  |      on yläkolmiomatriisi, jonka 
     |               |    |  |      päälävistäjän alkiot # 0
  r   ---------------------             
     |  0 ...       0| 0  |  0 B2    Tämä voi puuttua.    
     |  0 ...       0| 0  |  0      
  m   --------------------             

Ääretön määrä ratkaisuja, n-r vapaata parametria.


  
 Tärkeä seurauslause homogeenisille systeemeille Ax=0 

 Seuraus (Hsys)  
Jos homogeeniyhtälössä (b=0) on  m < n ,
niin aina on äärettömän monta ratkaisua. Sillä tällöinhän bmato=0
ja r <=m < n .



 Gauss-resepti 

 Askel 1: Oletetaan, että a[1,1] # 0 (uudelleenjärjestämällä).
          Valitaan kertoimet c[2]=-a[2,1]/a[1,1],...,c[m]=-a[m,1]/a[1,1]
          ja suoritetaan rivioperaatiot:
           rivi[2] <--  c[2]*r[1] + r[2]
           ...

           rivi[m] <--  c[m]*r[1] + r[m]

          Näin saadaan alkion a[1,1] alle 0-sarake.

 Askel 2: Jos amato[2,2] ja kaikki sen alapuoliset alkiot ovat
          sattumoisin = 0,  niin kaikki on hyvin. Tällöin syntyy
          "pitkä porras" ja siirrytään katsomaan samalla silmällä
          amato[2,3]:a ja sen alapuolista saraketta.
         
          Ellei, niin vaihdetaan tarvittaessa rivejä niin, että
          uusi amato[2,2] # 0. (Numeerisessa käytännössä vaihdetaan
          itseisarvoltaan mahdollisimman suuri.)
          Sitten jatketaan kohdasta [2,2] alkavalle alimatriisille
          samoin kuin Askel 1:ssä.

 Näin jatketaan niin kauan, kunnes kaikki on nollaa alapuolella tai
 törmätään alareunaan (r=m) .
 (Jos pysähdytään pystyseinään (r=n), niin silloinkin kaikki on
 nollaa alapuolella, eli tätä ei tarvitse sanoa erikseen lopetusehtona.)

 Huom  Yllä esiintyvää r-lukua ei sattumalta merkitä r:llä,
vaan se on (alkuperäisen) matriisin A rangi, eli säännöllisyysaste 
 "rank" . (Matlab:ssa ja Maplessa on molemmissa 
rank-funktio)

 Määritelmä, rangi, rank  Matriisin rangi r(A) on maksimimäärä
LRT sarakevektoreita.

 Lause  r(A) on 

 
  max määrä LRT rivivektoreita
 
Gaussin algoritmin luku r (alimman nollasta poikkeavan rivivektorin
  rivi-indeksi).



 Takaisinsijoitus 

Kun systeemi on saatu porrasmuotoon, saadaan ratkaisut heti lasketuksi,
tai nähdään välittömästi, että ratkaisuja ei ole.


Gaussin algoritmi on sikäli miellyttävä, että sen avulla saadaan sekä
ratkaisut lasketuksi (jos niitä on) että olemassaolo- ja lukumäärätieto
(ellei numeerisia
vaikeuksia oteta huomioon).

Systemaattinen tapa laskea ratkaisut porrasmuodosta lähtien on katsoa ensin alhaalla
olevat 0-rivit. Jos jonkun 0-rivin kohdalla bmato-sarakkeessa
on nollasta poikkeava luku, ei homma jatku, vaan todetaan, että ratkaisuja
ei ole. (Tapaus (a) yllä )

Jos ristiriitaisia yhtälöitä ei ole, käydään käsiksi 
alimmaiseen ei-0-riviin, rivi[r] .
Siitä lasketaan x[r] (joko 1-käs tai ottamalla vapaiksi parametreiksi
x[r+1], ... , x[n] ). Näin edetään alhaalta ylös.





 Matlab ja Maple 

Kun Gaussin algoritmia on joitakin kertoja harjoiteltu rivi riviltä,
 addrow  ja  swaprow  - tyyppisillä komennoilla
(Maple) tai Ab([i,k],:)=Ab([k i],:); 
Ab(i,:)=Ab(i,:)+c*Ab(k,:)  (Matlab)
(vrt. Perus-Gauss ), voidaan ryhtyä käyttämään
komentoja
     gausselim, backsub, gaussjord  (Maple)
     rref  (Matlab) 




 3.  Gauss-Jordan, redusoitu porrasmuoto, rref 

AG s. 60, KRE s. 366, Exa 1, GRE s. 13

Takaisinsijoitus voidaan myös toteuttaa rivioperaatioina, jolloin
päädytään muotoon, josta ratkaisut ovat luettavissa vielä suoremmin
kuin alakolmiomuodosta. Siitä käytetään nimitystä  redusoitu
porrasmuoto, "reduced row echelon form" (siitä lyhenne rref, jota mekin
käytämme).
Tässä nollataan yläkolmiota myös niin
paljon kuin mahdollista. Ei-singulaarisen neliösysteemin tapauksessa päädytään
diagonaalimatriisiin.

Tämä on erityisen hyödyllinen tapa silloin, kun on laskettava annetun
neliömatriisin käänteismatriisi.

 Redusoitua porrasmuotoa luonnehditaan näin 


   1. Jokaisen ei-nollarivin 1. nollasta poikkeava alkio = 1.
       Tätä sanotaan "johtavaksi ykköseksi".
  
 rivi[i+1]:n johtava 1 on aidosti oikealla rivi[i]:n johtavaa 1:tä.
  
 Sarake, jolla johtava 1 sijaitsee (johtava sarake), on muuten pelkkää
       nollaa.
  
 Kaikki 0-rivit ovat matriisin pohjalla.


Kaikki nämä matriisit ovat rref-muodossa, johtavat sarakkeet on merkitty *:llä.

 *    *    *       *    *            *         *

[1    0    0]     [1    0    2]     [1    5    0    0]
[           ]     [           ]     [                ]
[0    1    0]     [0    1    6]     [0    0    1    2]
[           ]     [           ]
[0    0    1]     [0    0    0]     [0    0    0    0]

                                    [0    0    0    0]
                              
     *         *    *    *

    [1    1    0    0    0     6]
    [                           ]
    [0    0    1    0    0    -1]
    [                           ]
    [0    0    0    1    0     3]
    [                           ]
    [0    0    0    0    1     1]


 Esimerkki rivioperaatioista Matlabilla 

format compact
format rational
A=[1 -1 1;-1 1 -1;0 10 25;20 10 0],b=[0;0;90;80]
Ab=[A b]
% 1. sarake
Ab(2,:)=Ab(2,:)+Ab(1,:)
Ab(4,:)=Ab(4,:)-20*Ab(1,:)
% Siirretään 0-rivi alimmaiseksi
Ab([4 2],:)=Ab([2 4],:)
% 2. sarake
Ab(3,:)=Ab(3,:)-1/3*Ab(2,:)
% Porrasmuoto on valmis ja ratkaisut voidaan kirjoittaa
% alhaalta ylöspäin:
x3=Ab(3,4)/Ab(3,3)
x2=(Ab(2,4)-Ab(2,3)*x3)/Ab(2,2)
x1=(Ab(1,4)-Ab(1,3)*x3-Ab(1,2)*x2)/Ab(1,1)
x=[x1;x2;x3]
% Tarkistus:
[A*x b]

 Komentojen tuottama tulostus: 

A =
      1           -1            1      
     -1            1           -1      
      0           10           25      
     20           10            0      
b =
      0      
      0      
     90      
     80      
Ab =
      1           -1            1            0      
     -1            1           -1            0      
      0           10           25           90      
     20           10            0           80      
Ab =
      1           -1            1            0      
      0            0            0            0      
      0           10           25           90      
     20           10            0           80      
Ab =
      1           -1            1            0      
      0            0            0            0      
      0           10           25           90      
      0           30          -20           80      
Ab =
      1           -1            1            0      
      0           30          -20           80      
      0           10           25           90      
      0            0            0            0      
Ab =
      1           -1            1            0      
      0           30          -20           80      
      0            0          95/3        190/3     
      0            0            0            0      
x3 =
      2      
x2 =
      4      
x1 =
      2      
x =
      2      
      4      
      2      
ans =
      0            0      
      0            0      
     90           90      
     80           80      
» 



» rref(Ab)
ans =
      1            0            0            2      
      0            1            0            4      
      0            0            1            2      
      0            0            0            0      


Tästä muodosta voidaan ratkaisut lukea suoraan.

 Takaisinsijoitus rivioperaatioilla, rref-muoto 

Edellä se jo näkyi, mutta tehdään "käsin".
Aloitetaan normeeraamalla rivit niin, että pivotit saadaan ykkösiksi.
Pudotetaan alin 0-rivi pois (sitä ei tarvita) ja tehdään elegantti pikku
temppu: Kerrotaan vasemmalta diagonaalimatriisilla, jonka diagonaalilla ovat
pivottien käänteisluvut. 

Ab=Ab(1:3,:)   % 0- rivi pois
d=diag(Ab)     % Poimitaan pivotit
D=diag(1./d)   % Muodostetaan kertojamatriisi.
Ab=D*Ab        % Elegantti normeeraus!
% 3. sarake
Ab(2,:)=Ab(2,:)-Ab(2,3)*Ab(3,:)
Ab(1,:)=Ab(1,:)-Ab(1,3)*Ab(3,:)
% 2. sarake
Ab(1,:)=Ab(1,:)-Ab(1,2)*Ab(2,:)

Näin on rref-muoto valmis, kuten tulostuksista näkyy:
>> Ab=Ab(1:3,:)
Ab =
      1           -1            1            0      
      0           30          -20           80      
      0            0          95/3        190/3     
>> d=diag(Ab)
d =
      1      
     30      
    95/3     
>> diag(1./d)
ans =
      1            0            0      
      0           1/30          0      
      0            0           3/95    
>> Ab=Ab(1:3,:)   % 0- rivi pois
Ab =
      1           -1            1            0      
      0           30          -20           80      
      0            0          95/3        190/3     
>> d=diag(Ab)     % Poimitaan pivotit
d =
      1      
     30      
    95/3     
>> D=diag(1./d)   % Muodostetaan kertojamatriisi.
D =
      1            0            0      
      0           1/30          0      
      0            0           3/95    
>> Ab=D*Ab
Ab =
      1           -1            1            0      
      0            1          -2/3          8/3     
      0            0            1            2 
>> Ab(2,:)=Ab(2,:)-Ab(2,3)*Ab(3,:)
Ab =
      1           -1            1            0      
      0            1            0            4      
      0            0            1            2      
>> Ab(1,:)=Ab(1,:)-Ab(1,3)*Ab(3,:)
Ab =
      1           -1            0           -2      
      0            1            0            4      
      0            0            1            2      
>> Ab(1,:)=Ab(1,:)-Ab(1,2)*Ab(2,:)
Ab =
      1            0            0            2      
      0            1            0            4      
      0            0            1            2  


 Maple ja Matlab 

Maplen linalg-pakkauksessa on funktio  gaussjord  ja Matlab:ssa
funktio  rref . 


 Käänteismatriisista ja sen laskemisesta 


Tässähän sitä hiukan


This page created by
<
Heikki.Apiola@hut.fi>

Last update 
10.9.2000