Ohjelmointi ja numeeriset menetelmät, luento 10, 29.3.

Ohjelman kokonaisrakenne

        program p      
        ! paaohjelma
        ...
        call sub
        end

        subroutine sub 
        ! aliohjelma
        ..
        end

Pääohjelma ja aliohjelmat toisistaan riippumattomia. Voivat olla peräkkäin samassa tiedostossa tai erillisissä tiedostoissa.

Pääohjelma ja aliohjelmat käännetään erikseen ja linkitetään yhteen.

Erillisiä aliohjelmia helppo linkittää yhteen.

Aliohjelmakirjasto on vain aliohjelmia sisältävä tiedosto.

Parametrit (viiteparametreja!) ovat vain osoitteita. Käännösaikana ei yleensä mahdollista tarkistaa, että todelliset parametrit ovat oikeaa tyyppiä tai että niitä on oikea määrä. Suoritus voi kaatua hämäräperäiseen ilmoitukseen (segmentation fault tms.)

Globaalit muuttujat on kerrottava erikseen jokaiselle aliohjelmalle common-alueen avulla.

Aliohjelmien työtilan varaaminen hankalaa. Tavallinen ratkaisu on varata työtilat kutsuvassa ohjelmassa ja välittää parametreina alihohjelmalle. Aliohjelmakutsut usein hyvin mutkikkaita.

Fortran 90

Tiedostossa voi olla:

1) Pääohjelma ja proseduureja

        program p      
        ! paaohjelma
        call sub
        ..
        end

        contains
        subroutine sub 
        ! aliohjelma

        end
        end program

Sisäiset proseduurit näkevät pääohjelman globaalit muuttujat automaattisesti.

Parametrien oikeellisuus tarkistettavissa.

2) Ulkoinen proseduuri

Erillisillä proseduureilla voi myös olla omia sisäisiä proseduureja:

        subroutine sub 
        ...
        call inner
        ...
        contains
          subroutine inner
          ...
          end
        end subroutine

Ei useita sisäkkäisiä tasoja, kuten Pascalissa.

3) Moduuli

        module m
        ! tyyppien ja muuttujien maarittelyt                
        ...
        contains
        subroutine sub 
        ! aliohjelma

        end
        ...
        end module

Myös moduulin proseduureilla voi olla omia sisäisiä proseduureja.

        module m2
        ...
        contains
        subroutine sub 
          ..
          contains
          subroutine inner
            ..
          end subroutine inner
        end subroutine sub
        ...
        end module

4) Mitä tahansa ohjelmatekstiä

Esimerkiksi usein toistuva osa voidaan kirjoittaa erilliseen tiedostoon ja ottaa mukaan include-lauseella:

        include 'fyle.f'

Erittäin tarpeellinen f77:ssa esimerkiksi globaalien muuttujien (common-alueiden) määrittelyssä. F90:ssä ei vastaavaa ongelmaa.

Paikallisten muuttujien säilyminen

Jos paikallinen muuttuja alustetaan määrittelyn yhteydessä, sillä on automaattisesti tämä ominaisuus. Sitä EI alusteta uudestaan joka kerta aliohjelmaan tultaessa:

    program savetest
      write(*,*) f(1.0)
      write(*,*) f(1.0)
    contains
    real function f(x)
      real, intent(in) :: x
      real :: y=0.0
      y=y+1
      f=x+y
   end function
   end program
  
   2.000000    
   3.000000    

Jos tällaista sivuvaikutusta ei haluta, muuttujat on alustettava sijoituslauseilla:

    real function f(x)
      real, intent(in) :: x
      real :: y
      y=0.0
      ...
 

Muuttujien attribuutit

Seuraava on täydellinen luettelo kaikista muuttujan attribuuteista.

dimension: määrittelee taulukon koon

parameter: muuttuja on vakio, jota ei saa muuttaa

save: muuttuja varataan staattisesti ja sen arvo säilyy aliohjelmakutsujen välillä

allocatable: dynaamisesti varattava taulukko

pointer: osoitin

target: muuttujaan saatetaan viitata osoittimen avulla

intent: aliohjelman argumentin käyttötapa

optional: valinnainen argumentti; voi puuttua proseduurin kutsusta

external: muuttuja on proseduurin nimi

intrinsic: muuttuja on varusfunktion nimi

private: moduulin yksityinen muuttuja

public: moduulin julkinen muuttuja

Osittaisdifferentiaaliyhtälöt

Kaikki edellä esiintyneet yhtälöt ovat tavallisia differentiaaliyhtälöitä. Kahden tai useamman muuttujan yhtälö, joka sisältää osittaisderivaattoja, on osittaisdifferentiaaliyhtälö.

Hyvin monia ilmiöitä voidaan kuvata toisen kertaluvun lineaarisilla yhtälöillä

a (\pd²u / \pd²x) + b (\pd²u / \pd x \pd y) + c (\pd²u / \pd²y) + d (\pd u \pd x) + e (\pd u \pd y) + fu + g = 0.

Tällainen yhtälö on jotakin kolmesta tyypistä riippuen toisten derivaattojen kertoimista:

1. Elliptinen, jos b²-4ac < 0.
2. Parabolinen, jos b²-4ac = 0.
3. Hyperbolinen, jos b²-4ac > 0.

Jaottelulla ei sinänsä suurta merkitystä numeriikan kannalta, mutta:

• Hyperboliset ja paraboliset yhtälöt kuvaavat usein jonkin ajallista kehitystä, ja ovat siten alkuarvotehtäviä
• Elliptiset yhtälöt kuvaavat usein "staattista" tilannetta, jolla on tietyt reuna-arvot.

Tärkeimmät ratkaisumenetelmät:

• differenssimenetelmät
• elementtimenetelmät

Yksinkertainen osittaisdifferentiaaliyhtälö on esimerkiksi Laplacen yhtälö, joka suorakulmaisessa koordinaatistossa on

\pd²u / \pd²x + \pd²u / \pd²y = 0.

Laplacen yhtälössä a=c=1 ja b=0, joten b²-4ac = -4, ja kyseessä on siis elliptinen yhtälö.

Jos oikea puoli ei ole nolla, saadaan Poissonin yhtälö

\pd²u / \pd²x + \pd²u / \pd²y = f(x,y).

Esimerkiksi gravitaatiopotentiaali V toteuttaa Poissonin yhtälön

\pd²V / \pd²x + \pd²V / \pd²y + \pd²V / \pd²z = \nabla² V = -4\pi G\rho(x,y,z),

missä \rho on tiheys ja G gravitaatiovakio.

Elliptiset yhtälöt kuvaavat usein ajasta riippumatonta potentiaalia tai tasapainossa olevaa systeemiä.

Differenssimenetelmä / elliptiset yhtälöt

  \pd u(xi,yj) / \pd x = (u(xi+1,yj)-u(xi-1,yj)) / 2\Delta x,
  \pd u(xi,yj) / \pd y = (u(xi,yj+1)-u(xi,yj-1)) / 2\Delta y,
  \pd2 u(xi,yj) / \pd2 x = (u(xi+1,yj)-2u(xi,yj)+u(xi-1,yj)) / (\Delta x)2,
  \pd2 u(xi,yj) / \pd2 y = (u(xi,yj+1)-2u(xi,yj)+u(xi,yj-1)) / (\Delta y)2.

  \pd2  ui,j / \pd2 x = (ui+1,j-2ui,j+ui-1,j) / h2,
  \pd2 ui,j / \pd2 y = (ui,j+1-2ui,j+ui,j-1) / h2.

Laplacen yhtälö on nyt

\nabla² u_i,j = (u_i+1,j-2u_i,j+u_i-1,j) / h² + (u_i,j+1-2u_i,j+u_i,j-1) / h² = 0

eli

(1 / h²)(u_i+1,j+u_i-1,j+u_i,j+1+u_i,j-1-4u_i,j) = 0.

Laplacen operaattori voidaan esittää kaaviona

            1             1
 \nabla2 = ---        1  -4   1
           h2             1

            1        1     4    1
 \nabla2 = ---       4    -20   4
           6h2       1     4    1

            1              1    1
 \nabla2 = ---       1    -6   1
            h2       1     1

Tylsä vakioesimerkki: suorakaiteen muotoinen levy, jonka yksi reuna pidetään 100 asteen lämpötilassa ja muut reunat 0 asteessa.

Dirichlet'n reunaehto: ratkaistavan funktion arvo tunnetaan koko ratkaisualueen reunalla.

Käytetään vain kolmea hilapistettä:

          0    0    0 
     0    u1   u2   u3   100 
          0    0    0 

(1 / h²)(0+0+u₂+0-4u₁) = 0
(1 / h²)(u₁+0+u₃+0-4u₂) = 0
(1 / h²)(u₂+0+100+0-4u₃) = 0

eli

  | -4   1   0  |  | u1 |   |   0 |
  |  1  -4   1  |  | u2 | = |   0 |
  |  0   1  -4  |  | u3 |   | 100 |

  u1   u2   u3   u4   u5  u6 
  u7   u8   u9   u10  u11  u12
  u13  u14  u15   u16  u17  u18 

      ...           
      1            1 -4   1               1
            1        1   -4   1              1 
               1          1  -4  1              1 
      ...           

Eliminoinnissa matriisin vinorivien väliselle alueelle muodostuu nollasta poikkeavia alkioita.

Jos pisteitä on paljon, on edullista käyttää iteratiivisia ratkaisumenetelmiä, jolloin em. tapauksessa tarvitsee tallettaa vain matriisin viisi vinoriviä.

Liebmannin iteraatio

  | -4   1   0 |  | u1 |   |   0 |
  |  1  -4   1 |  | u2 | = |   0 |
  |  0   1  -4 |  | u3 |   | 100 |

u₁ = (u₂) / 4,
u₂ = (8u₁+u₃) / 4,
u₃ = (u₂+100) / 4.

Näillä iteraatiokaavoilla päivitetään ratkaisuvektoria u, kunnes ratkaisu ei enää muutu vaadittua tarkkuutta enempää. Iterointi lopetetaan esim., kun suurin u:n alkioiden muutoksista on riittävän pieni.

Laplace-operaattorin kaaviosta saadaan suoraan iteraatioon sopivat kaavat:

                      |      1    |
 \nabla2 u = (1 / h2) | 1   -4   1 |  uij=0,
                      |      1    |

u_ij = (1/4)(u_i,j+1+u_i,j-1+u_i-1,j+u_i+1,j).

Alkuarvo voidaan laskea esimerkiksi ratkaisemalla yhtälö hyvin karkealla hilalla. Ratkaisun avulla interpoloidaan likimääräiset arvot tiheämmän hilan pisteissä.

Liebmannin iteraation suppeneminen voi olla hidasta. Suppenemista voi kiihdyttää ylirelaksaatiomenetelmällä.

Merkitään yläindeksiksi iteraatiokierros, jolla arvo on laskettu:

u_ij^k+1 = (1/4) (u_i,j+1^k+u_i,j-1^k+1+u_i-1,j ^k+1+u_i+1,j^k).

Lisätään ja vähennetään u_ij^k:

u_ij^k+1 = u_ij^k + ( {u_i,j+1^k+u_i,j-1^k+1+u_i-1,j ^k+1+u_i+1,j^k - 4u_ij^k) / 4.

Sulkulauseke on korjaus, jolla u_ij:tä päivitetään. Kun ratkaisu on löytynyt, sulkulauseke on nolla.

Korjaus on usein liian pieni, joten u:n arvot lähestyvät hitaasti oikeaa ratkaisua. Suppeneminen nopeutuu, jos korjaustermiä suurennetaan sopivasti:

u_ij^k+1 = u_ij^k + \omega ( u_i,j+1^k+u_i,j-1^k+1+u_i-1,j^k+1+u_i+1,j^k - 4u_ij^k) / 4.

Vakio \omega on välillä [1, 2]. Kun pisteiden määrä kasvaa hyvin suureksi, optimaalinen arvo -> 2.

(Suorakulmaisen hilan ja Dirichlet'n reunaehtojen tapauksessa optimaalinen arvo on

\omega = 4 / (2+\sqrt(4-c²)),

missä

c = cos \pi/n + cos \pi/m,

ja n ja m ovat hilapisteiden määrä x- ja y-suunnissa.)

Differenssimenetelmät/paraboliset yhtälöt

\pd u / \pd t = c \pd² u / \pd x²

alkuehdolla u(x,t=0)=g(x) ja reunaehdoilla u(x=0,t)=a, u(x=1,t)=b.

\pd²u \pd x² = (1/h²) (u_i+1^k - 2u_i^k + u_i-1^k),
\pd u / \pd t = (1/\Delta t) (u_i^k+1-u_i^k).

Diskretoitu yhtälö

u_i^k+1 = u_i^k + (c \Delta t / h²) (u_i+1^k - 2u_i^k + u_i-1^k).

Oikean puolen termi u_i^k häviää, kun 2c\Delta t/ h² = 1. Voidaan osoittaa, että tästä saadaan raja stabiilisuudelle. Jotta yhtälö olisi stabiili, on oltava 2c\Delta t/h² < 1. Jos siis pituus\-askelta lyhennetään, on samalla lyhennettävä aika-askelta.

Menetelmä on eksplisiittinen: kullakin aika-askelella tarvitaan vain edellisellä aika-askelella laskettuja arvoja. Virhe on O(\Delta t)+O(h²).

Implisiittiset menetelmät ovat ajan suhteen tarkempia. Esimerkiksi Crank-Nicholsonin menetelmä:

u_i^k+1 = u_i^k + (c\Delta t / 2h²) (u_i+1^k+1 - 2u_i^k+1+u_i-1^k+1+u_i+1^k-2u_i^k + u_i-1^k).

Virhe on O((\Delta t)²)+O(h²), ja menetelmä on stabiili kaikilla \Delta t.

Haittapuoli on, että kutakin aika-askelta kohti joudutaan ratkaisemaan tridiagonaalinen yhtälöryhmä.

Differenssimenetelmät/hyperboliset yhtälöt

\pd²u / \pd t² = c² \pd²u / \pd x²

Diskretoimalla saadaan

(u_i^k+1-2u_i^k+u_i^k-1) / (\Delta t)² = c² (u_i+1^k-2u_i^k+u_i-1^k\over h²)

eli

u_i^k+1 = 2u_i^k - u_i^k-1 + (c²(\Delta t)² / h²) (u_i+1^k-2u_i^k+u_i-1^k).

Jos aika-askel valitaan siten, että c\Delta t = h, saadaan

u_i^k+1 = -u_i^k-1 + u_i+1^k+u_i-1^k.

Yhtälö voi kuvata esimerkiksi värähtelevää kieltä. Alkuarvoina tarvitaan poikkeamat alkuhetkellä, u_i⁰, ja nopeus, josta voidaan laskea u_i¹.

Elementtimenetelmät

Käytetään erityisesti rakenneanalyysissä ja virtausmekaniikassa.

• Jaetaan käsiteltävä alue pieniin elementteihin. (Verkon generointi)
• Kunkin elementin alueella kuvataan ratkaisua yksinkertaisen kantafunktion (esim. polynomin) avulla. Kukin kantafunktio on nollasta poikkeava vain pienellä alueella.
• Yhtälön ratkaisu saadaan kantafunktioiden lineaarikombinaationa.
• Kantafunktiot muodostetaan niin, että kahden elementin rajalla molemmat kantafunktiot saavat samat arvot elementtien solmupisteissä (esimerkiksi kolmikulmaisen elementin kärjissä). Lisäksi voidaan vaatia, että myös derivaatoilla on samat arvot solmupisteissä.
• Kantafunktioiden kertoimet saadaan optimointitehtävän ratkaisuna. Ratkaisun laskeminen johtaa (yleensä suureen) lineaariseen yhtälöryhmään.

Elementtien muodot ja koot voidaan valita melko vapaasti tehtävän mukaan. Menetelmää on helppo soveltaa myös tilanteisiin, joiden geometria on mutkikasta.

Elementtien kokoa voidaan helposti pienentää kriittisissä kohdissa.

Kysymyksessä on minimointitehtävä, joka on yleensä aina ratkaistavissa. Menetelmä on stabiilimpi kuin differenssimenetelmät.

Funktionaali

F[v] = \int v(x,y,y') dx

saa ääriarvon, kun v toteuttaa Lagrangen yhtälön

\pd v / \pd y - (d/dx) (\pd v /\pd y') = 0.

Valitaan

v = (y'(x))² + 2f(x)y(x),

jolloin

\pd v / \pd y = 2f(x),
\pd v / \pd y' = 2y'
(d/dx) (\pd v / \pd y') = 2y''.

Lagrangen yhtälö on siis

2f(x)-2y''(x)=0

eli

y''=f(x).

Differentiaaliyhtälö

y''=f(x)

voidaan siis ratkaista myös variaatiotehtävänä etsimällä minimi funktionaalille

F[v] = \int y'²+2f(x)y dx

Vastaavasti esimerkiksi Poissonin yhtälö

\nabla² u = f(x,y,z)

voidaan ratkaista minimoimalla funktionaalia

F[u] = \int ( (\pd u/ \pd x)² + (\pd u/ \pd y)² + (\pd u/\pd z)² +2f(x,y,z) u) dx dy dz.

Approksimoidaan ratkaisua kantafunktioiden \phi_i lineaarikombinaationa:

v = \sum_j a_j \phi_j,

jolloin minimoitava funktionaali on

F[v] = \int ( (\sum_j a_j\phi'_j)²+2f(x)\sum_j a_j \phi_j) dx.

Minimi saadaan yhtälöistä

\pd F[v] / \pd a_i = 0.

\pd F[v] / \pd a_i =
\int ( 2\phi'_i (\sum_j a_j\phi'_j) + 2f(x) \phi_i ) dx
= 2 \sum_j a_j \int \phi'_i\phi'_j\,dx + 2 \int f(x)\phi_i dx,

josta

\sum_j a_j \int\phi_i'\phi_j'dx = - \int f(x)\phi_i.

Saadaan yhtälöryhmä

Fa = b,

missä

F_ij = \int\phi_i'\phi_j' dx ,
b_i = -\int f(x) \phi_i dx.

Esimerkki: yritetään ratkaista yhtälö y''=x reunaehdoilla y(0)=y(1)=0.

Valitaan kantafunktioiksi

         1+(x-xi)/h, kun xi-1 <= x <= xi
 \phii = 1-(x-xi)/h, kun xi <= x <= xi+1
         0  muuten

Kerroinmatriisin alkiot ovat

                       -1/h, kun i=j +- 1,
 \int \phi'i\phi'j dx = 2/h, kun i=j,
                        0 muuten

 \int01 f\phi1\,dx = 
  \int0h x (1 + (x-h)/h) dx +
  \inth2h x (1 - (x-h)/h) dx 
 = h2/3 + 2h2/3 = h2,
 \int01 f\phi2 dx = 2h2,
 \int01 f\phi3 dx = 3h2.

Valitaan yksinkertaisuuden vuoksi h=0.25, jolloin kantafunktioita on vain kolme.

Minimointitehtävästä saadaan yhtälöryhmä

         |  2  -1   0 | | a1 |           | 1 |
 1/0.25  | -1   2  -1 | | a2 | =  -0.252 | 2 |
         |  0  -1   2 | | a3 |           | 3 |

     | -0.0391 | 
 a = | -0.0625 |
     | -0.0547 |

Yhtälön tarkka ratkaisu on y=(x³-x)/6, josta y(0.5)=-0.0625.

Ominaisarvot

Ax = \lambda x

eli

(A - \lambda I) x = 0.

Jotta yhtälöllä olisi ei-triviaali ratkaisu, on oltava

det (A - \lambda I) = 0.

Kun determinantti kirjoitetaan auki, saadaan karakteristinen polynomi, josta ominaisarvot voidaan ratkaista.

Jos matriisi on reaaliarvoinen ja symmetrinen, ominaisarvot ovat reaalisia.

Esimerkiksi matriisin

   |  1   4 |
   |  1   1 |

       |  1-\lambda    4 |
   det |                 | =   = \lambda2 - 2\lambda - 3,
       |  1     1-\lambda|

QR-hajotelma

A = QR,

missä Q on ortogonaalimatriisi ja R yläkolmiomatriisi.

Muunnos voidaan tehdä eri tavoilla:

1. Householderin muunnos
2. Givensin rotaatiot
3. Gram-Schmidt-ortogonalisointi

Householderin muunnos

     
      |  3  2  1  |
 A  = |  1  1  2  |
      |  2  1  3  |

               | 3 |
 x1 = a(:,1) = | 1 |
               | 2 |

              | 1 |   | -0.7416574 |
  u1=x1-||x1|| | 0 | = | 1          |          
              | 0 |   |  2         |

P1=I - 2 (u1u1^T / || u1 ||2) =
  | 0.8017837   0.2672612   0.5345225 |
  | 0.2672612   0.6396433  -0.7207135 |
  | 0.5345225  -0.7207135  -0.4414270 |

A1= P1 A =
 |  3.7416574   2.4053512    2.9398737 |
 |  0           0.4534522   -0.6155927 |
 |  0          -0.0930955   -2.2311854 |

                 |  0.4534522 |
 x2 = a(2:3,1) = |            |
                 | -0.0930955 |

             |  1  |   | -0.0094578 |
 u2=x2-||x||2 |     | = |            |
             |  0  |   | -0.0930955 |

P2=I - 2u2u2T / || u2 ||2 = 
      |  1   0           0         |
      |  0   0.9795688  -0.2011093 |
      |  0  -0.2011093  -0.9795688 |

A2=P2 A1 = 
 | 3.7416574    2.4053512    2.9398737 |
 | 0            0.4629100   -0.1543033 |
 | 0            0            2.3094011 |

 Q = P1 P2  =
 | 0.8017837    0.1543033  -0.5773503 | 
 | 0.2672612    0.7715167   0.5773503 |
 | 0.5345225   -0.6172134   0.5773503 |

 R = P2 P1 A =
 | 3.7416574    2.4053512   2.9398737 |
 | 0            0.4629100  -0.1543033 |
 | 0            0           2.3094011 |

  
  QR = 
   |  3   2   1  |
   |  1   1   2  |
   |  2   1   3  |

 Q QT = 
   |  1   0   0  |
   |  0   1   0  |
   |  0   0   1  |

  Q = P1 P2 ... Pn, 
  R = Pn ... P2 P1 A.

Ominaisarvojen laskeminen

Neliömatriisi on lohkokolmiomatriisi, jos se on muotoa

 | T11  T12  T13  ... T1n |
 | 0    T22  T23 ...  T2n |
 |                       |
 | 0    0    0   ... Tnn |

Voidaan osoittaa, että tällaisen matriisin ominaisarvot ovat lävistäjälohkojen T_ii ominaisarvot.

=> jos matriisi on lävistäjä- tai kolmiomatriisi, ominaisarvot ovat lävistäjällä olevat alkiot.

Ominaisarvot eivät muutu similariteettimuunnoksissa

A -> S^-1 A S.

Muunnetaan matriisi aluksi johonkin yksinkertaisempaan muotoon similariteettimuunnoksilla.

Matriisin muuntaminen Householderin muunnoksilla

 A  =
 |  4    3    2    1 |
 |  3    4   -1    2 |
 |  2   -1    1   -2 |
 |  1    2   -2    2 |

      |  3  |
 x1 = |  2  |
      |  1  |

                  | 1 |   | -0.7416574 |
  u1 = x1 - ||x1|| | 0 | = | 2          |
                  | 0 |   |  1         |

  p1 = I - 2u1u1T/ ||u1|| = 
  | 0.8017837    0.5345225    0.2672612 |
  | 0.5345225   -0.4414270   -0.7207135 |
  | 0.2672612   -0.7207135    0.6396433 |

 P1  =
 | 1    0           0            0         |
 | 0    0.8017837   0.5345225    0.2672612 |
 | 0    0.5345225  -0.4414270   -0.7207135 |
 | 0    0.2672612  -0.7207135    0.6396433 |

 A1 = P1 A P1 =
  |  4           3.7416574   0           0         |
  |  3.7416574   2.4285714   1.2977396   2.1188066 |
  |  0           1.2977396   0.0349563   0.2952113 |
  |  0           2.1188066   0.2952113   4.5364723 |

       |  1.2977396 |
 x2  = |            |
       |  2.1188066 |

                    | -1.1869072 |
  u2 = x2 - ||x2|| = |            |
                    |  2.1188066 |

 p2  =
  |  0.5223034    0.8527597 |
  |  0.8527597   -0.5223034 |

  P2 =
  |  1    0    0           0         |
  |  0    1    0           0         |
  |  0    0    0.5223034   0.8527597 |
  |  0    0    0.8527597  -0.5223034 |

 
 A2 = P2 A1  P2 =
 |  4           3.7416574    0            0         |
 |  3.7416574   2.4285714    2.4846467    0         |
 |  0           2.4846467    3.5714286   -1.8708287 |
 |  0           0           -1.8708287    1         |

Muutetaan sitten alkuperäinen matriisi epäsymmetriseksi:

  A  =
   | 4    2    3    1  |
   | 3    4   -2    1  |
   | 2   -1    1   -2  |
   | 1    2   -2    2  |

A2  =
 |  4           3.4743961  -1.3039935  -0.4776738 |
 |  3.7416574   1.7857143   2.9123481   1.1216168 |
 |  0           2.0934787   4.2422252  -1.0307759 |
 |  0           0.         -1.2980371   0.9720605 |

H =  
   | x  x  x  x  x |
   | x  x  x  x  x |
   | 0  x  x  x  x |
   | 0  0  x  x  x |
   | 0  0  0  x  x |

QR-algoritmi

Valitaan aluksi A₁=H. Sitten toistetaan seuraavia askelia:

-- Lasketaan QR-hajotelma Q_iR_i = A_i.

-- Lasketaan uusi matriisi A_i+1 = R_i Q_i.

QR-hajotelman matriisi Q on ortogonaalinen ja A=QR, R=Q^TA, joten

RQ = Q^T A Q

on similariteettimuunnos ja säilyttää ominaisarvot.

Matriisien A_i jono suppenee kohti yläkolmio- tai lohkomatriisia, josta ominaisarvot voidaan poimia.

Viimeisen esimerkin tapauksessa raja-arvo on 50 iteraatiokierroksen jälkeen

 | 7.0363389  -0.7523758  -0.7356716  -0.3802631 |
 | 4.265E-08   4.9650342  -0.8892339  -0.4538061 |
 | 0           0          -1.9732687  -1.3234202 |
 | 0           0           0           0.9718955 |