A numerikus módszerek néhány fejezete

Jeney, András

Új Széchenyi Terv logó.

Miskolci Egyetem

Kelet-Magyarországi Informatika Tananyag Tárház

Kivonat

Nemzeti Fejlesztési Ügynökség http://ujszechenyiterv.gov.hu/ 06 40 638-638

Lektor

Dr. Kálovics Ferenc

Miskolci Egyetem, Analízis Tanszék, ny. egyetemi docens, a matematikai tudomány kandidátusa

A tananyagfejlesztés az Európai Unió támogatásával és az Európai Szociális Alap társfinanszírozásával a TÁMOP-4.1.2-08/1/A-2009-0046 számú Kelet-Magyarországi Informatika Tananyag Tárház projekt keretében valósult meg.


Tartalom

1. Bevezetés
2. Mátrixok, vektorok
2.1. Mátrix- és vektorműveletek
2.2. Speciális mátrixok, vektorok
2.3. Mátrixok particionálása
2.4. Lineáris függetlenség, mátrixok rangja
2.5. Mátrixok determinánsa és inverze
2.6. Vektor- és mátrixnormák, mátrixok kondíciószáma
2.7. Feladatok
3. Klasszikus hibaszámítás
3.1. Alapfogalmak
3.2. Az alapműveletek abszolút és relatív hibái
3.3. Függvények elsőrendű hibabecslése
4. Lebegőpontos hibaszámítás
4.1. A lebegőpontos aritmetika modellje
4.2. Kerekítési hibák és becslésük
4.3. A kerekítési hibák halmozódásának kompenzálása
4.4. A lebegőpontos aritmetikai szabvány
5. Érzékenység, numerikus stabilitás
5.1. Függvények kondíciószáma
5.2. Direkt és inverz hibák
5.3. Feladatok
6. Lineáris egyenletrendszerek megoldása
6.1. Direkt módszerek
6.1.1. Háromszögmátrixú egyenletrendszerek
6.1.2. A Gauss-módszer
6.1.3. A Gauss-módszer műveletigénye
6.1.4. A főelemkiválasztásos Gauss-módszer
6.1.5. A Gauss-Jordan elimináció, mátrixinvertálás
6.1.6. Az LU- és a Cholesky-módszer
6.2. Iteratív módszerek
6.2.1. A Jacobi-iteráció
6.2.2. A Seidel-iteráció
6.3. Feladatok
7. Egyenletrendszerek hibaanalízise
7.1. Érzékenységvizsgálat
7.1.1. Inverz hiba a jobboldali vektorban
7.1.2. Inverz hiba az együtthatómátrixban
7.1.3. Inverz hiba az együtthatómátrixban és a jobboldalon is
7.2. A relatív hiba becslése Wilkinson tételével
7.3. Utólagos hibabecslés
7.4. Iteratív javítás
7.5. Feladatok
8. A MATLAB-ról, röviden
8.1. Mátrixok és vektorok a MATLAB nyelvben
8.2. Utasítások, függvények és eljárások a MATLAB nyelvben
8.2.1. Utasítások
8.2.2. Függvények
8.2.3. M-adatállományok, eljárások
9. Előadásvázlatok fóliákon
9.1. Fólia
9.2. Fólia
9.3. Fólia
9.4. Fólia
9.5. Fólia
9.6. Fólia
9.7. Fólia
9.8. Fólia
9.9. Fólia
9.10. Fólia
9.11. Fólia
9.12. Fólia
9.13. Fólia
9.14. Fólia
9.15. Fólia
9.16. Fólia
9.17. Fólia
9.18. Fólia
9.19. Fólia
9.20. Fólia
9.21. Fólia
9.22. Fólia
9.23. Fólia
9.24. Fólia
9.25. Fólia
9.26. Fólia
9.27. Fólia
9.28. Fólia
9.29. Fólia
Irodalomjegyzék

Az ábrák listája

2.1. Háromszögmátrixok geometriai sémája
2.2. Sávmátrix geometriai sémája
2.3. Indukált norma geometriai jelentése
5.1. Numerikusan stabil és instabil függvények
5.2. A direkt és az inverzhiba kapcsolata
5.3. Nem érzékeny és érzékeny feladat
5.4. Kerekítési hibák hatása érzékeny feladatnál
6.1. A Gauss elimináció sémája
6.2. Felső sávszélesség növekedése elimináció során
6.3. Alsó háromszögmátrixból kettős tükrözéssel felső
6.4. Iteráció 1. lépése
6.5. Iteráció 2. lépése
6.6. Iteráció 3. lépése
9.1. Direkt hiba, inverz hiba (Direkt hiba, inverz hiba (dx=\Delta x , dy=\Delta y ), Direkt hiba, inverz hiba (dx=\Delta x , dy=\Delta y ))
9.2. Kerekítési hibák halmozódása, I.
9.3. Kerekítési hibák halmozódása, II.
9.4. Kerekítési hibák halmozódása, III.
9.5. Kerekítési hibák halmozódása, IV.
9.6. Kerekítési hibák halmozódása, V.
9.7. Kerekítési hibák halmozódása, VI.
9.8. Alsó háromszögmátrixú egyenletrendszer sémája
9.9. Felső háromszögmátrixú egyenletrendszer sémája
9.10. Az Az LU -felbontás sémája-felbontás sémája

1. fejezet - Bevezetés

A matematikát bizonyos szempontok alapján fel szokták osztani tisztán elméleti matematikára és alkalmazott matematikára. Meg kell azonban jegyezni, hogy élesen nem lehet elválasztani a két területet, hiszen a gyakorlatban közvetlenül hasznosítható matematikai eredmények mögött elméleti háttér van, ugyanakkor a gyakorlat is felvet olyan kérdéseket, amelyek elméleti kutatásokat indíthatnak el. Azt sem állíthatjuk, hogy a ma még tisztán elméletinek tűnő eredmények soha nem kerülnek be a gyakorlatba. Mindenesetre ha van tudományterület, amelynek gyakorlati alkalmazhatóságát senki nem vonhatja kétségbe, a numerikus módszerek területe bizonyára az. (Azért megjegyzendő, hogy 200 évvel ezelőtt, mondjuk, egy 10000 ismeretlent tartalmazó lineáris egyenletrendszer megoldására irányuló kutatást is feltehetőleg haszontalannak ítélték volna sokan.) Elméleti alapjainak fejlesztése mellett pontosan a gyakorlat által felvetett feladatok inspirálják, adnak újabb és újabb irányt a kutatásoknak. Ugyanakkor a numerikus módszereket a véges sok lépésben véges sok számjeggyel való számítások miatt szokták a kerekítés művészetének is nevezni. A numerikus módszereket több helyen is „közelítő számítások“ címszó alatt tárgyalják, amin csak azok csodálkoznak, akik nem veszik figyelembe, hogy számok pontos leírására – akármilyen számrendszerben – legtöbbször végtelen sok számjegyet kellene alkalmazni, ami – számítógép ide vagy oda – lehetetlen.

Egy konkrét gyakorlati feladat megoldását – kissé elnagyoltan – a következő fázisokra lehet bontani:

(i) A feladat megfogalmazása, célkitűzés.

(ii) Szakmai modell felállítása.

(iii) Matematikai modell kialakítása.

(iv) Algoritmus készítése és hibaelemzés.

(v) Számítógépi program készítése, a számítások végrehajtása.

(vi) Az eredmények értékelése és esetleg visszacsatolás: valamelyik fázistól módosítás, majd onnan a folyamat újraindítása.

A numerikus módszerek alapvető feladata a (iv) fázis, azaz kész matematikai modellek megoldására elfogadható hibahatáron belüli számszerű erdményeket produkáló eljárások kidolgozása. Ugyanakkor nem választható el az (v) fázistól sem élesen – sőt, egyre szélesedik a vele való átfedés –, hiszen figyelembe kell venni a végrehajtáshoz szükséges időt, a rendelkezésre álló eszközöket: a számítógépeket, az igénybe vehető szoftvereket. A matematikai modellalkotáshoz is köze van, mert egy-egy algoritmus gyakran az eredeti modellt is módosítja.

Könyvünk anyagának összeállítását elsősorban gyakorlati szempontok vezérelték, ám a szükséges definíciók, tételek kimondását sem mellőzhettük. Szemelvényszerűen néhány bizonyítást is közlünk a tárgyalt tétel jobb megértése céljából, egyben bemutatva az alkalmazható eszközök egyikét-másikát.

A könyv a numerikus módszereknek csupán két nagyobb területére ad betekintést: a klasszikus és a lebegőpontos hibaanalízis területére (ami a számítások végeredményében jelentkező hibák keletkezésének vizsgálatát és nagyságának becslését jelenti), valamint a lineáris egyenletrendszerek leggyakrabban alkalmazott megoldási módszereibe és azok hibaanalízisébe. Követelmény, hogy a hallgató a tárgy felvételekor alapvető lineáris algebrai ismeretekkel rendelkezzen, bár a mátrixokkal kapcsolatos néhány fogalmat szükségesnek tartottunk itt is megemlíteni. Feltételezzük továbbá, hogy a hallgató némi általános számítástechnikai tudás birtokában van. A gyakorlatokat számítógép mellett javasoljuk megtartani, ezért (is) igyekeztünk bemutatni az egyes részek tárgyalásánál a számítógépen való végrehajtást: több algoritmusnál megadtuk azok egy lehetséges MATLAB kódját. Az algoritmus jobb megértését segítendő céllal akkor is, ha létezik rá beépített MATLAB függvény. (Különben is: vannak, akik az algoritmus fogalmát úgy értelmezik, hogy az egy számítógépi kód.) A példákat a 6.5. verzióval oldottuk meg. Tapasztalat szerint érdemes a közölt programok futtatása során bizonyos részeredményeket is kinyomtatni, ami szintén a jobb megértést szolgálja. Megjegyezzük, hogy a szükséges, alapvető MATLAB ismeretek a tárgy keretében is megtaníthatók, kevés idő ráfordítással; könyvünk utolsó fejezete ezt kívánja megkönnyíteni.

Számos kidolgozott típus- és magyarázó példát is beépítettünk, valamint gyakorlásra szánt bőséges feladatgyűjteményt is kínálunk.

2. fejezet - Mátrixok, vektorok

Bár – mint ahogy a bevezetésben is utaltunk rá – feltételezzük a mátrixok, vektorok fogalmának és a velük kapcsolatos műveletek ismeretét, szükségesnek tartjuk rövid áttekintésüket és a könyvben alkalmazott jelölések rögzítését. Mivel könyvünkben a rövidebb és egységesebb tárgyalásmód érdekében a vektorokat speciális mátrixoknak tekintjük, ezért – sok könyvtől eltérően – először a mátrixok definíciójával kezdjük.

2.1. Definíció. Az típusú (méretű) valós mátrixon valós számok alábbi táblázatát értjük:

Az és az értelemszerűen pozitív egész számok. Beszélünk a mátrix -edik soráról és -edik oszlopáról. A sorok és oszlopok metszéspontjában vannak a mátrix elemei (). Komplex számokból (sőt, más – absztrakt – struktúrák elemeiből) is felépíthetünk mátrixokat.

Az típusú valós mátrixok halmazát jelöli, ennek megfelelően például az azt jelenti, hogy egy típusú (méretű) valós mátrix.

2.2. Definíció. Az egyetlen sorból vagy egyetlen oszlopból álló mátrixot vektornak nevezzük.

A sorvektor szokásos megadási módja: . Az a vektor egy eleme, másképpen komponense. Az oszlopvektorokat az

formában szoktuk megadni, ahol az komponensű oszlopvektorok halmaza (tulajdonképpen , bár itt inkább izomorfiáról beszélhetnénk).

2.3. Megjegyzés. Ha egyszerűen csak vektort említünk és nem tesszük hozzá, hogy az sor- vagy oszlopvektor-e, akkor azon mindig oszlopvektort értünk. A komponensek kettős indexelése egyrészt felesleges (az egyik mindig 1), másrészt ez is mutatja, hogy valójában a mátrixoktól függetlenül is értelmezhető matematikai objektum.

A mátrixok nevét (azonosítóját) a könyvben a latin nagybetűk közül választjuk (pl. , ); elemeinek azonosítására pedig általában a megfelelő kisbetűket használjuk, feltüntetve a sor- és oszlopindexeket (pl. , ). Bizonyos esetekben a mátrixelemeket is nagybetűvel adjuk meg, ilyenkor az indexeket zárójelek közé írjuk, pl. , .

A definíciókban látható szögletes zárójelek helyett szokás kerek zárójeleket is használni, valamint gyakran célszerű a következő, tömörebb jelölést választani:

Ha a sorok és oszlopok száma megegyezik, akkor a mátrixot négyzetesnek nevezzük, tömörebb jelölése is némileg egyszerűsödik. Például az méretű mátrix esetén:

2.1. Mátrix- és vektorműveletek

Ebben az alfejezetben röviden emlékeztetünk a mátrixok, vektorok körében bevezetett fontosabb műveletekre és azok tulajdonságaira. Általában mátrixokról beszélünk, és mivel a vektorokat is (speciális méretű) mátrixoknak tekintjük, ahol nem lényeges, ott nem említjük azokat külön is.

I. Összeadás. Csak azonos méretű mátrixok között értelmezzük, mégpedig elemenként. Azaz, ha , akkor

Az összeadás fontos tulajdonságai: kommutatív és asszociatív (a tagok sorrendje felcserélhető és csoportosítható), invertálható (adott mátrixokhoz létezik , hogy teljesül), továbbá van zéruselem (lásd a speciális mátrixoknál).

II. Számmal való szorzás. Legyen és legyen valós szám (azaz ). Ekkor

Nyilvánvaló, hogy , továbbá a fenti két művelet értelmezéséből következik az alábbi két – disztributívitást kimondó – szabály:

2.4. Megjegyzés. Az halmazt (vagy az halmazt) a fenti két művelettel algebrai struktúrának tekinthetjük és mint ilyen, rendelkezik mindazon tulajdonságokkal, amelyek a lineáris teret (más szóval – itt ugyan szokatlanul hangzik – vektorteret) definiálják; például van zéruselem, az összeadás invertálható, stb. Ezért nem lesz meglepetés, ha az általunk érzékelt háromdimenziós térben definiált vektorok bizonyos jellemzőit (pl. hosszúság) általánosítjuk többdimenzióban is, sőt ezt az általánosítást kiterjesztjük mátrixokra is.

III. Transzponálás (tükrözés). Az transzponáltját jelölje , amit a következőképpen definiálunk:

Úgy is fogalmazhatjuk: a sorokat és az oszlopokat felcseréljük. Fontos tudni a definícióból adódó szabályt.

Négyzetes mátrixok esetén a transzponálás a főátlóra való tükrözést jelent. (A főátlót azon elemek alkotják, melyek sor- és oszlopindexe megegyezik, azaz az elemek.)

A transzponálás felhasználásával az oszlopvektorokat meg lehet adni még az , a sorvektorokat pedig az

formában is.

IV. Szorzás. Ha , és , akkor a szorzatukat a

előírással definiáljuk.

Látható, hogy amint az összeadásnál, úgy itt is lényeges a műveletben szereplő operandusok (tényezők) mérete, csak itt más a követelmény: az első tényező oszlopainak száma és a második tényező sorainak a száma kell, hogy egyenlő legyen; az eredmény mérete a tényezőkéből következik. Megemlítjük, hogy esetén van egységelem (lásd a speciális mátrixoknál), a szorzás asszociatív és az összeadással összekapcsolva disztributív; de nem kommutatív, tehát általában (a két mátrix között többnyire nem is értelmezhető mindkét sorrendben a szorzás) és általában nem is invertálható. Az összeg transzponálásához analógiát mutat az szabály.

Ha akkor az szorzás neve skalár szorzás. A skalár szorzat eredménye egy skalár szám (-es mátrix), szokták külön is definiálni.

2.5. Definíció. Az skaláris szorzata:

A skalár szorzás definícióját (2.1)-gyel összevetve látható, hogy a szorzatmátrix bármelyik indexű elemét úgy kapjuk, mint egy skalár szorzatot: az első tényező -edik sorát (mint sorvektort) szorozzuk a második tényező -edik oszlopával, azaz

A továbbiakban a mátrix és mátrix-vektor műveletek felírásánál feltesszük, hogy az ott szereplő mátrixok, ill. vektorok méretei olyanok, amelyek lehetővé teszik az adott művelet elvégzését.

2.2. Speciális mátrixok, vektorok

Mátrixokat sokféle szempontból osztályozhatunk. Itt most néhány olyan mátrixot, vektort sorolunk fel, melyeknek meghatározóan sok elemük azonos (gyakran éppen zérus) és ezek szabályosan helyezkednek el. Az ebből eredő különleges tulajdonságaikat a későbbiekben többször fogjuk használni.

2.6. Definíció. Az mátrix szimmetrikus, ha .

Nyilvánvaló, hogy csak négyzetes mátrix lehet szimmetrikus és ekkor ().

2.7. Definíció. Ha az mátrix összes eleme zérus, akkor -t zérusmátrixnak nevezzük:

A zérusmátrix az összeadásra nézve a zéruselem, és szorzásra nézve is az „elvárt” tulajdonságú: azaz megfelelő méretű mátrixok esetén

2.8. Definíció. Az mátrix egységmátrix, ha

Az egységmátrix (amit gyakran a magyar elnevezésének kezdőbetűjével, -vel is jelölünk) a szorzásra nézve egységelem, azaz minden esetén

2.9. Megjegyzés. Szokás az egységmátrix fogalmát a nem négyzetes mátrixokra is kiterjeszteni, a definícióját úgy adva meg, hogy az azonos indexpárú elemei 1-esek, a többi zérus (például ). Ilyenkor persze az egység elnevezés már nem jogos, félrevezető. Mindenesetre nagyon hasznos kiterjesztés, több matematikai szoftver (a MATLAB is) él vele.

Az egységmátrix mellett fontos fogalom az egységvektor.

2.10. Definíció. Az vektort (-edik) egységvektornak nevezzük, ha az -edik komponense 1-es, a többi pedig zérus.

A transzponáltját felírva, tehát:

2.11. Definíció. A mátrixot diagonálmátrixnak nevezzük, ha esetén .

A diagonális elemek indexei azonosak, ezért egyetlen indexszel hivatkozhatunk rájuk: . Magát a diagonálmátrixot gyakran a vagy () formában is jelöljük.

2.12. Megjegyzés. Szokták a diagonális fogalmat is kiterjeszteni, úgy tekintve, mint azon mátrixelemek összességét, melyeknél az oszlopindex és a sorindex különbsége állandó, és ezt a számot nevezzük a diagonális sorszámának. Így például a főátló (az ún. fődiagonális), a -adik az e fölötti átlók sorszáma pozitív, az alattiaké pedig negatív. Példa rá az alábbi -es mátrix, ahol a -edik és az -edik átló elemeit írtuk be, a többi elemet -gal jelöltük:

2.13. Definíció. A mátrix permutációmátrix, ha minden sorában és oszlopában pontosan egy darab 1-es van és a többi elem zérus.

Például az alábbi mátrix egy -es permutációmátrix:

2.14. Definíció. Az mátrix alsó háromszög alakú, ha minden esetén és felső háromszög alakú, ha minden esetén .

Az alsó háromszögmátrixok alakja sematikusan a következő

a felső háromszögmátrixoké pedig

Gyakran csak ábrával jelöljük, hogy (alsó vagy felső) háromszögmátrixról van szó:

2.1. ábra - Háromszögmátrixok geometriai sémája

Háromszögmátrixok geometriai sémája

2.15. Megjegyzés. Beszélni szoktunk valamely négyzetes mátrix alsó és felső háromszögrészéről, amelyekre a , illetve jelöléssel hivatkozunk. Ezeket értelemszerűen úgy kapjuk, hogy az eredeti mátrix főátlója feletti, illetve alatti elemeit kicseréljük zérusokra. Például az mátrix esetén

(Az az angol lower, pedig az upper szó kezdőbetűje. A szakirodalomban több helyen találkozhatunk az angol elnevezésből eredő jelöléssel; ilyen például az egységmátrixra használt jelölés is.)

2.16. Definíció. Az mátrix sávmátrix alsó sávszélességgel és felső sávszélességgel, ha teljesül, hogy

Más szóval: a ()-edik diagonálisuk alatti és -adik diagonálisuk feletti elemeik zérusak. A sávot, amelynek elemei lehetnek nem zérusok is, azon elemek alkotják, amelyek indexeire teljesül, hogy , vagy ekvivalens módon . Sematikusan, illetve ábrán is megmutatva:

2.2. ábra - Sávmátrix geometriai sémája

Sávmátrix geometriai sémája

Például a diagonális mátrix egy alsó és felső sávszélességű mátrix. Végső soron minden mátrix felfogható lenne sávmátrixként is, de nyilván csak akkor érdekesek, ha és jóval kisebb mint . Tárolásukra is speciális, helytakarékos módok léteznek, hiszen a zérusokat nem kell tárolni, ha tudjuk, hogy azok hol szerepelnek.

2.17. Megjegyzés. Általában ritka mátrixoknak nevezzük azokat a mátrixokat, amelyek viszonylag sok, ismert pozíciójú zérust tartalmaznak. Ilyenek pédául a sávmátrixok, vagy lehetnek szabálytalan (de ismert, rögzített) elhelyezkedésű sok zérust tartalmazó mátrixok is. Az ilyen mátrixok tárolása a zérusok figyelmen kívül hagyásával helytakarékosan oldható meg, és a velük való különböző manipulációkat végrehajtó programok is gazdaságosan, művelettakarékosan írhatók meg. (Persze, tudni kell, hogy közben nem változnak-e meg a zérusok, illetve azt, hogy legfeljebb hol változhatnak.) Ebből a szempontból ritkának tekinthetők a háromszögmátrixok is, a nemzérus elemek tárolása történhet vektorban, például oszlopfolytonos sorrendben. A szimmetria is tekintetbe vehető; ilyen esetben elég csak az alsó (vagy a felső) háromszögrészt tárolni. (Itt is ügyelni kell a programozás során, hogy hol romlik el esetleg az algoritmus végrehajtása során a szimmetria.)

2.3. Mátrixok particionálása

A particionálás részekre bontást jelent. Geometriai analógiával: ha egy négyzetrácsos („kockás”) papírra írnánk fel az mátrixot, akkor annak ollóval vízszintes és függőleges, széltől-szélig (ún. guillotine vágásokkal) való szétszabdalását jelenti. Egy papírdarabkára eső része a mátrixnak az egy partíciója, amit -vel jelölünk. Például, ha darab vízszintes és darab függőleges csíkra vágjuk a mátrixot, akkor

és ha történetesen az -edik vízszintes csík az eredeti mátrix sorindexű, valamint a -edik függőleges csík a oszlopindexű részét foglalja magában, akkor

ahol és .

A partíciót blokknak is szoktuk nevezni. Az értelmezésből következik, hogy az azonos sorban álló blokkok sorainak száma azonos. Hasonlóképpen, az azonos oszlopban álló blokkok oszlopainak száma azonos. Továbbá értelemszerűen fennállanak az

összefüggések.

A fenti részletes jelölés helyett használhatjuk a tömörebb alakot, azaz a sor- és oszlopindex helyére a megfelelő „tól-ig“ indexpárokat írjuk. Ha egyetlen sorból vagy oszlopból áll a partíció, akkor csak azt az egyet írjuk, illetve, ha teljes sor- vagy oszlophosszúságú a partíció, azt a kettőspont önálló használatával is jelölhetjük. Pédául:

(A legutóbbi a teljes mátrixot jelenti.)

Azonosan particionált mátrixok összegzését és skalárral való szorzását blokkonként végezhetjük, úgy mintha a blokkok számok lennének. Particionált mátrixok szorzását is végezhetjük blokkok szerint, de csak akkor, ha az első mátrix függőleges particionálása megegyezik a második tényező vízszintes particionálásával, hogy a kívánt blokkszorzatok egyáltalán kiszámíthatóak legyenek. Például az

particionált mátrixok szorzata a particionálás megtartásával akkor lehetséges, ha esetén . Ekkor (a példánkban) nyilván a többi megkívánt blokkszorzás is elvégezhető, mert például és . Így

Leggyakrabban sorok, illetve oszlopok szerint particionálunk, tehát, ha , akkor

Szokás még az mátrix -edik sorvektorát -vel, a -edik oszlopvektorát -vel jelölni, ezekkel

Mindenesetre ügyeljünk arra, hogy ez a jelölés félreértésre adhat okot. Tudniillik itt az és az nem egymás transzponáltja; valójában semmi közük egymáshoz, hanem

A sorok, illetve oszlopok szerinti particionálások felhasználásával az és mátrixok szorzata felírható

alakban is. Tehát a sorok és oszlopok skalárszorzataiból álló mátrix. Az mátrixszorzatot felírhatjuk még a következőképpen is:

Az egységvektorokkal való szorzás eredménye pedig:

azaz a másik tényező megfelelő sor-, illetve oszlopvektorát kapjuk eredményként.

Ezek után a diagonálmátrixszal való szorzás eredménye is egyszerűen belátható: a sorok, illetve az oszlopok szorzódnak a megfelelő diagonális elemmel. Legyenek az mátrixok sorok, illetve oszlopok szerint particionálva, ekkor tehát:

Minthogy a permutációmátrix sorai transzponált egységvektorok, oszlopai pedig maguk az egységvektorok, permutációmátrixszal való szorzás a másik tényező sorainak, illetve oszlopainak permutálását eredményezi (ez magyarázza is a permutációmátrix elnevezését):

ahol az és a az , illetve az számok egy-egy permutációja.

A permutációmátrix néhány további tulajdonságával kapcsolatban a 2.7. alfejezet 6. feladatát ajánljuk figyelmükbe.

2.18. Megjegyzés. A MATLAB-ban a partíciók (részmátrixok) fogalmánál sokkal általánosabb módon is hivatkozhatunk az adott mátrix bizonyos elemeire, illetve azokból felépített másik mátrixra. Ezeket a lehetőségeket a 8.1. alfejezetben tárgyaljuk. Itt most két példát mutatunk be; az első:

illetve a második: ha fel kell cserélni egy mátrix -edik és -adik sorát, akkor ezt a

utasítással elérhetjük. Ezt a lehetőséget mi is fogjuk használni az algoritmusok leírásánál.

2.4. Lineáris függetlenség, mátrixok rangja

A későbbiekben szükségünk lesz az itt következő néhány fogalomra. Tekintsük a darab, egyenként komponensű vektort. Gyakran használjuk az jelölést is. Legyenek továbbá valós számok.

2.19. Definíció. A kifejezést az vektorok lineáris kombinációjának nevezzük.

2.20. Definíció. Az vektorok lineárisan függetlenek, ha a csak a triviális esetén teljesül.

Megemlítjük, hogy esetén ezzel ekvivalens: lineárisan függetlenek, ha egyik vektor sem állítható elő a többi lineáris kombinációjaként.

2.21. Definíció. Az vektorok közül kiválasztható maximális számú lineárisan független vektorok számát a vektorrendszer rangjának nevezzük.

2.22. Definíció. Egy mátrix rangja az oszlopvektoraiból (vagy sorvektoraiból) álló vektorrendszer rangja.

Nem magától értetődő, hogy a mátrix oszlop-, illetve sorvektorrendszer rangja megegyezik, de igaz. A rang jelölésére a kifejezést fogjuk használni.

Végezetül egy fontos megjegyzéssel zárjuk ezt a részt.

2.23. Megjegyzés. Az vektorok összes lineáris kombinációja lineáris teret (vektorteret) alkot (lásd a 2.4. megjegyzést is). A vektortér dimenziója éppen a közülük való maximálisan sok lineárisan független vektorok darabszáma és egy ilyen független vektorsokaság a tér egy bázisa. Maga a lineáris tér fogalma sokkal általánosabb: az összeadás és a skalár számmal való szorzás bevezetésével és annak tulajdonságaival kapcsolatos. Lineáris teret alkotnak például (a megfelelő műveletekkel) az méretű mátrixok vagy éppen az intervallumon folytonos függvények is.

2.5. Mátrixok determinánsa és inverze

Jelölje azt az -es mátrixot, amelyet az első sora és a -edik oszlopa elhagyásával kapunk, tehát particionálva:

2.24. Definíció. A

előírásokkal számított valós számot a négyzetes mátrix determinánsának nevezzük.

Megjegyezzük, hogy többnyire másként szokták a determináns fogalmát bevezetni, itt most ekvivalens módon, az ún. kifejtési tétellel (egyfajta kiszámítási szabállyal) definiáltuk.

2.25. Definíció. Az mátrixot a négyzetes mátrix inverzének nevezzük, ha , ahol az egységmátrix.

Ha az inverz mátrix létezik, akkor egyértelmű. Az inverz mátrix jelölése . Az inverz mátrixra fennállnak az alábbi tulajdonságok:

Azokat a mátrixokat, melyeknek létezik inverze, nemszinguláris mátrixoknak nevezzük.

2.26. Tétel. Az mátrixnak akkor és csak akkor van inverze, ha .

Értelemszerűen, ha , akkor a mátrixot szingulárisnak mondjuk.

2.6. Vektor- és mátrixnormák, mátrixok kondíciószáma

Mint már korábban említettük, az akárhány komponensű vektorok (de az összeadás és a skalár számmal való szorzás bevezetésével értelmezett lineáris térben akár a mátrixok is) az általunk érzékelhető, geometriailag irányított egyenes szakasszal ábrázolható vektorok általánosításaként foghatók fel. Ezért bizonyos fogalmak is általánosíthatók, a leglényegesebb tulajdonságaik megtartásával. Ilyen fogalom a hosszúság, amelyet a lineáris térben normának nevezünk.

2.27. Definíció. Az függvényt mátrixnormának (vagy vektornormának, ha ) nevezzük, ha

A norma szokásos jelölése: .

Normát nyilván nagyon sokféleképpen adhatunk meg. Bár a vektorokat is (speciális méretű) mátrixokkal azonosítottuk, vannak szempontok, amelyek alapján mégis csak különböző objektumokról van szó, ezért konkrét normákat különböző módon vezetünk be. Vektorok esetén a normák fontos osztályát alkotják az ún. hatványnormák:

ahol egész szám. Ezek közül is a leggyakrabban használt vektornormák a következők:

A legutóbbit a határátmenettel kapjuk. (Az norma megfogalmazásánál természetesen a jobboldalon az abszolút érték jel elhagyható valós vektorok esetén, de akkor nem, ha az elemek komplexek.)

2.1. Példa. Legyen . Ekkor

,

, (kerekítve),

.

A háromkomponensű vektorok által bezárt szög skalár szorzás segítségével való kiszámításának szabályát kiterjesztve, értelmezhetjük az akárhány komponensű vektorok közötti szöget is.

2.28. Definíció. Az () vektorok szöge , amelynek koszinuszát a

összefüggés definiálja.

Legyen . A leggyakrabban használt mátrixnormák a következők:

(A Frobenius-normánál itt is csak komplex elemek esetén kell a jobboldalon az abszolút érték jel.)

2.29. Megjegyzés. A spektrálnormát más fogalom segítségével is szokták definiálni, itt csak a legismertebbet, a sajátértékkel valót közöltük. (Egy mátrix sajátértékének azt a valós vagy komplex számot nevezzük, melyhez létezik vektor, hogy fennáll.) Ennek a normának sok helyen különösen fontos szerepe van, ugyanakkor a meghatározása rendkívül nehéz feladat. A Frobenius-norma egyik nagy jelentősége, hogy felülről becsli a spektrálnormát (a hétköznapi életben is gyakran segít, ha egy távolság nagyságának legalább valamilyen felső korlátját ismerjük).

2.2. Példa. Legyen .

Ezen mátrix normái:

,

,

.

(Az .)

Az egyes mátrixnormák lábindexelését indokolhatjuk formálisan is: ha az objektumot -beli vektornak tekintjük és úgy számoljuk az említett normáit, ugyanakkor -beli mátrixként is meghatározzuk a megfelelő normákat, akkor azok megegyeznek. A mátrixnormák alábbi, fontos osztálya azonban a jelölésnek is mélyebb magyarázatát adja. (Itt és a továbbiakban is négyzetes mátrixokra szorítkozunk.)

2.30. Definíció. A mátrixnormát a vektornorma által indukált mátrixnormának nevezzük, ha

Az indukált mátrixnorma geometriai jelentése: az egységnormájú vektorok megnyújtásának maximális mértéke:

2.3. ábra - Indukált norma geometriai jelentése

Indukált norma geometriai jelentése

Könnyen igazolható, hogy az , az , pedig az vektornorma által indukált mátrixnorma.

2.3. Példa. Felhasználva az indukált mátrixnorma definícióját, igazoljuk, hogy esetén .

Megoldás. Az értelmezés szerint

Tehát a feltételes szélsőérték feladatot kell megoldanunk (). Analitikus eszközökkel könnyen előállítható a megoldás: . Eredményünket az egyenlőséglánc jobboldalába helyettesítve megkapjuk a példa állítását.

2.31. Tétel. Indukált mátrixnormában és ( és ).

2.32. Megjegyzés. Az állítás nem minden mátrixnormára igaz. Ha csak a 2.31. tételben szereplő első egyenlőtlenség teljesül, akkor azt mondjuk, hogy a két norma kompatibilis. A definíciók alapján világos, hogy az indukált mátrixnormák az őket indukáló vektornormákkal kompatibilisek, de a Frobenius norma – bár nem indukált norma – is kompatibilis az vektornormával, azaz

Végül megemlítjük a mátrixok kondíciószámát.

2.33. Definíció. Legyen nemszinguláris. A mennyiséget az mátrix kondíciószámának nevezzük.

A kondíciószám normafüggő, nyilván csak nemszinguláris mátrixok esetén értelmezhető. Bevezetésének motivációját a 7.1. alfejezetben adjuk meg.

2.7. Feladatok

1. Legyen és . Rögzített értékekhez adjuk meg -t és -et úgy, hogy mind az , mind a szorzás elvégezhető legyen.

2. Melyik mátrix szimmetrikus bármely mellett a következők közül?

(i) ,

(ii) .

3. Igazoljuk a következőket:

(i) két alsó (vagy felső) háromszögmátrix szorzata is alsó (vagy felső) háromszögmátrix,

(ii) egy háromszögmátrix determinánsa a főátlóbeli elemek szorzata,

(iii) alsó (vagy felső) háromszögmátrix inverze is alsó (vagy felső) háromszögmátrix (ha létezik).

4. Legyen . Igazoljuk, hogy , .

5. Legyen diagonálmátrix. Tudva azt, hogy bármely diagonálmátrix sajátértékei megegyeznek a fődiagonális elemeivel, mutassuk meg, hogy

6. Legyen tetszőleges permutációmátrix.

(i) Mutassuk ki, hogy determinánsának abszolút értéke ,

(ii) igazoljuk, hogy az inverze megegyezik a transzponáltjával,

(iii) adjuk meg valamennyi tanult normáját,

(iv) igazoljuk, hogy ha nem az egységmátrix, akkor létezik olyan részmátrixa, amelynek determinánsa (),

(v) igazoljuk, hogy bármely két (azonos méretű) permutációmátrix szorzata is permutációmátrix.

7. Tekintsük az teret. Ábrázoljuk az , , és az „egységsugarú köröket” az kordinátarendszerben.

8. Igazoljuk, hogy mátrixnorma.

Mutassuk meg továbbá, hogy az mátrixokra .

9. Igazoljuk, hogy .

10. Felhasználva, hogy a -es mátrixnorma a -es vektornorma által indukált norma, igazoljuk a (2.13) egyenlőtlenséget.

11. A (2.3) mintájára írjunk egyetlen utasítást, amely az mátrixban elvégzi az , cseréket.

3. fejezet - Klasszikus hibaszámítás

3.1. Alapfogalmak

Valós számmal kifejezhető minden adatnak van egy elméleti – pontos – értéke; pédául a tér két adott pontja közötti távolság valamilyen mértékegységben egy elméleti érték. Ugyanígy a is létezik elméletileg, de hozhatnánk példákat az élet bármely területéről. Ezeket a pontos értékeket rendszerint nem ismerjük, csupán valamekkora hibát tartalmazó közelítéseiket (becsléseiket). Előző példáink mellett maradva, mondhatjuk, hogy Miskolc és Budapest két kijelölt pontja között a távolság légvonalban kb. 160 km. vagy, hogy . A közelítések akkor érnek valamit, ha a tévedésnek valamilyen felső korlátját ismerjük. Ilyen korlátok rendszerint megadhatók. A hibát és a hibakorlátot a következőképpen definiáljuk.

3.1. Definíció. Legyen pontos érték, pedig annak valamilyen közelítése. A mennyiséget az közelítés hibájának nevezzük, a számot pedig az abszolút hibájának. Azt a értéket pedig, amelyre fennáll, hogy , az abszolút hibakorlátjának mondjuk.

3.2. Megjegyzés. A előjele nem lényeges, gyakran (ha levezetésekben, bizonyításokban úgy kényelmesebb) fordított előjellel használjuk, a végcél úgyis mindig az abszolút értékének becslése. Felhívjuk továbbá a figyelmet arra, hogy a későbbiekben abszolút hiba helyett egyszerűen csak hibát mondunk, sőt rendszerint a hibakorlátot értjük rajta.

A definíció értelmében használjuk az hivatkozást is, ami annyit jelent, hogy . Nyilván annál jobb a közelítés, más szóval annál élesebb a becslés (és erre törekedni kell), minél kisebb a .

Az abszolút hiba sok esetben semmitmondó, vagy éppen félrevezető. Például ugyanaz a abszolút hibakorlátú közelítés egész más pontosságot jelent, ha egy elméletileg -res nagyságrendű érték közelítéséről van szó, mint ha a becsült érték nagyságrendje . A közelítés jóságát ezért az abszolút hiba és az abszolút hibának a pontos érték egységére eső része – a relatív hiba – együtt jellemzi.

3.3. Definíció. Az szám valamely közelítő értékének relatív hibája a mennyiség.

Minthogy az pontos érték általában nem ismeretes, ezért a helyett a közelítést használjuk. Az így elkövetett hiba elhanyagolható, ha értelmes becslésről van szó, azaz és lényegesen nagyobb a másodrendben kicsiny mennyiségnél.

Szokás a relatív hiba helyett annak százalékos értékét megadni, azaz .

Egy matematikai modell megoldása (például egy integrál kiszámítása) közben számítások egész sorát végezzük el: a bemenő (input) adatokból egy-egy bonyolult algoritmus sok-sok aritmetikai művelet végrehajtásával állítja elő a végeredményt, az ún. output adatokat. A végeredmény az eredetileg kitűzött probléma elméleti megoldásához képest hibát tartalmaz. Ezen hibákat keletkezésük oka szerint a következőképpen osztályozhatjuk:

analitikus hibák: az eredeti matematikai modell és a numerikusan ténylegesen megoldott feladat közötti eltérésből adódó hiba. Ilyen hiba keletkezik például, ha nem az eredeti függvény integrálját, hanem valami egyszerűbben számolható közelítő függvényt integrálunk; vagy ha egy végtelen sort az első néhány tagjával közelítünk. Ez utóbbit csonkolási hibának is szokás nevezni.

öröklött hibák: az input adatok hibáinak az eredményre pontosan végzett számítások mellett gyakorolt hatásai.

kerekítési hibák: a számításokat véges sok tizedesjeggyel (vagy más – például -es – számrendszerbeli jeggyel) végezzük, ezért közben kerekíteni kell. Ezek a kerekítések a végeredményben szükségszerűen hibákat okoznak.

Az analitikus hibákat az egyes eljárások kidolgozása során (mint azoknak a módszer használhatóságát jellemző tulajdonságát) kell vizsgálni. Az öröklött hibákkal elsősorban a klasszikus hibaanalízis, a kerekítési hibákkal pedig a lebegőpontos hibaanalízis foglalkozik. A kettő azonban nem függetleníthető egymástól.

A klasszikus hibaszámítás keretein belül tehát azt vizsgáljuk, hogy amennyiben az elméleti értékek helyett az adott hibakorlátú közelítésekkel – de azokkal pontosan – végeznénk a számításokat, úgy a kapott eredmény (legfeljebb) mennyivel térne el az ismeretlen elméleti végeredménytől. Mivel minden aritmetikai számítás végső soron valós számok között a négy alapművelet valamilyen sorozatának a végrehajtása (akkor is, ha azokat előttünk „rejtett“ , a számoló-, számítógépbe beégetett algoritmusok végzik; például a kiszámítását), ezért alapvető ezen műveletek eredményeinek öröklött hibáit ismerni.

A következő jelöléseket és elnevezéseket használjuk: pontos értékek, és a közelítő értékeik és hibakorlátokkal, azaz és .

3.2. Az alapműveletek abszolút és relatív hibái

Jelölje a , , , műveletek bármelyikét. Az művelet eredményét az elméleti eredmény közelítésének tekintjük és a

illetve a

becsléseket keressük, ahol a művelet hibáját, pedig abszolút hibakorlátját jelöli. Az additív műveletek (összeadás, kivonás) hibaszámítás szempontjából egymás között hasonlóságot mutatnak, ezért egyetlen tételben adjuk meg a megfelelő hibakorlátokat.

3.4. Tétel. Az additív műveletek abszolút hibakorlátjai a következők:

Bizonyítás.

amiből a fenti állításunk következik.

Mivel mindkét művelet esetén ugyanazt az eredményt kaptuk, valójában az előjelükre semmilyen kikötést nem kellett tenni. Az eredmény akárhány, tetszőleges előjelű tagra kiterjeszthető. Tekintsük a

összegzést. Könnyen belátható, hogy . Természetesen ez az esetek nagy részében jelentősen túlbecsli a tényleges abszolút hibát, hiszen azt tételezi fel, hogy az egyes tagok hibáinak előjele a legkedvezőtlenebbül alakul. Valószínűségszámítási eszközökkel élesebb becslés is adható, jó megbízhatósággal.

3.5. Tétel. A multiplikatív műveletek abszolút hibakorlátjai a következők:

Bizonyítás. A szorzat abszolút hibakorlátjára kapjuk, hogy

Ha és , akkor a másodrendű hibatagot elhanyagolhatjuk és azzal éppen az állításunkat kapjuk.

Az osztás esetén természetesen feltesszük, hogy a nevező nem zérus és azt kapjuk, hogy

Itt pedig hasonló meggondolással a tagot hanyagolhatjuk el az mellett, amivel állításunk kiadódik.

3.6. Megjegyzés. Az osztás abszolút hibakorlátja -hoz közeli esetén rendkívül nagy lehet, ezért algoritmusainkat úgy alakítsuk, hogy lehetőleg minél nagyobb abszolút értékű számmal kelljen osztani!

Rátérve a relatív hibákra, a következő állítást igazolhatjuk, kikötve, hogy a nevező sehol sem lehet zérus, az additív műveleteknél pedig az operandusok megegyező előjelét is előírva.

3.7. Tétel. Az aritmetikai műveletek relatív hibakorlátjai a következők:

Bizonyítás. Tulajdonképpen csak az összeadással kell foglalkoznunk, a kivonásra adott összefüggés megegyezik a definícióval, a szorzás és osztás pedig behelyettesítés után azonnal adódik.

Az utolsó egyenlőség az és azonos előjeléből következik.

Az additív műveletekhez hasonlóan, itt a multiplikatív műveletekre lehet több tényező esetén egyszerű relatív hibakorlátot adni ():

Megemlítjük, hogy a kivonásra is le lehet vezetni olyan relatív hibakorlátot, amelyben az operandusoknak is csak a relatív hibái szerepelnek, de annak árát nem érdemes megfizetni: amellett, hogy túl bonyolult, jelentősen durvul a korlát (az összeadásnál is durvítottunk kicsit, de az még – úgymond – megéri).

3.8. Megjegyzés. A kivonás relatív hibája, amennyiben az eredmény zérushoz közeli, igen nagy lehet, ezért kerüljük el az ilyen kivonást! Rendkívül „alattomos“ hatása lehet: az ábrázolt értékes jegyek száma csökkenhet, ami nem tűnik fel (szemben a zérushoz közeli számmal való osztással, ahol már esetleg figyelmeztetést is kapunk, gyakran le is áll a program), de kerekítési hibaként halmozódva teljesen hamis végeredményt produkálhat.

Összefoglalva: Additív műveletek esetén az abszolút hibakorlátok, multiplikatív műveletek esetén a relatív hibakorlátok adódnak össze.

3.1. Példa. Két ellenállást műszerrel megmértünk és a következő értékeket kaptuk: , . A párhuzamos kapcsolással kapott eredő ellenállást az ismert képlettel számoljuk. Határozzuk meg az input adatok relatív hibakorlátjait és az eredő ellenállás közelítő értékét. Számítsuk ki a közelítő érték abszolút és relatív hibakorlátját kétféleképpen is:

(i) az input adatoknak csak az abszolút korlátjait használva a -t és ezután a relatív korlátot,

(ii) az input adatoknak csak a relatív korlátjait használva a relatív korlátot és ezután a abszolút korlátot.

Magyarázzuk meg, miért térnek el az eredmények.

Megoldás. , az eredő ellenállás közelítőleg .

(i) , a relatív hiba: .

(ii) .

Az eltérések a két számítás során alkalmazott különböző elhanyagolások, illetve eltérő becslések miatt adódnak.

3.2. Példa. Ismertek a és közelítő értékek, amelyek közös abszolút hibakorlátja , a közös relatív hibakorlát pedig . Számítsuk ki a közelítő értékét és annak relatív hibakorlátját.

Megoldás. A kivonás elvégzésével kapjuk, hogy , amelynek relatív hibakorlátja az általános formulából

azaz (a kiinduló adatok relatív hibáinak közel -szerese). Most lehetőségünk van az elméleti relatív hiba kiszámolására is, ami „csak“ . Ez a valóságos hiba is jelentős mértékű, a kiinduló adatokéihoz képest kb. -szoros. A különbség képzését elkerülhetjük a

átalakítással. A számláló pontos érték. A nevező abszolút hibája , a hányados relatív hibája pedig , azaz . Ez összhangban van az input relatív hibáival és lényegesen kisebb, mint amit a közvetlen kivonásnál kaptunk.

Hasonló fogásokat lehet alkalmazni más esetekben is.

3.3. Függvények elsőrendű hibabecslése

Külön foglalkozunk az egy- és a többváltozós esetekkel. Legyen legalább kétszer folytonosan differenciálható függvény, . Az helyett -t számoljuk. Az

másodrendű Taylor-formulából kapjuk, hogy

ahol (). A másodrendű tagot elhanyagolva kapjuk, hogy a függvénybehelyettesítés abszolút hibája

Rátérve a többváltozós esetre, legyen legalább kétszer folytonosan differenciálható függvény és , valamint , ahol . A többváltozós Taylor-formulából az egyváltozós esethez hasonlóan a másodrendű tagot elhanyagolva, kapjuk:

ahol . Ebből pedig adódik a

becslés.

3.9. Megjegyzés. Ezzel az (ún. lineáris) hibabecsléssel azonban óvatosan kell bánni, a másodrendű tag (az , ahol a Hesse mátrix valamely közbülső helyen) elhanyagolása gyakran nem jogos.

Függvények relatív hibája értelemszerűen

3.3. Példa. Oldjuk meg a 3.1. példát harmadikféleképpen is, az eredő ellenállást a két ellenállás (kétváltozós) függvényének tekintve.

Megoldás. , . Ezeket a közelítő értékekkel kiszámolva, majd behelyettesítve (3.10)-be, kapjuk, hogy . A relatív hiba pedig: .

Feltűnő ezen eredménynek a többitől való kissé nagyobb mértékű eltérése, ami demonstrálja, hogy az elsőrendű (másképpen: lineáris) hibabecsléssel szemben gyakran fenntartással kell élni (bár jelen példánkban az egyező nagyságrendek miatt éppen el is fogadhatjuk az eredményt).

4. fejezet - Lebegőpontos hibaszámítás

4.1. A lebegőpontos aritmetika modellje

A gyakorlatban csak véges sok számjeggyel dolgozhatunk, így a számítógépek is csupán egy véges számhalmazt ábrázolnak és az aritmetikai műveleteket ezekkel a számokkal végzik, mégpedig a tudományos-műszaki számítások zömét ún. lebegőpontos aritmetikában. Ennek legáltalánosabban elfogadott modellje a következő.

4.1. Definíció. A lebegőpontos számok halmaza

ahol

A -át az feltétel miatt kellett külön a halmazhoz adni.

A három leggyakrabban használt számrendszer a következő:

Minthogy a mantissza , egy -beli számot felírhatunk az

alakban is, ahol , de (az feltétel miatt). Az ilyen számrendszereket normalizáltaknak nevezzük. A normalizált alak mantisszájának jegyei az értékes számjegyek.

Az halmaz elemeinek száma könnyen meghatározható: az említett feltételek miatt féle mantissza és féle kitevő lehet, egy szám pozitív és negatív előjelet is viselhet, azonkívül külön számoljuk a -át; tehát összesen szám alkotja -et.

A mantissza egész része zérus (az miatt); ezt a -át és a tizedes- (bináris, stb.) pontot értelemszerűen nem szokás ábrázolni. Ha , akkor az első jegy csak az lehet, amelyet szintén felesleges ábrázolni.

Az halmaz elemei nem egyenletesen helyezkednek el a számegyenesen.

4.1. Példa. Legyen , , és . Adjuk meg az elemeit.

Megoldás. A mantissza lehetséges értékei (binárisan): és vagy tízes számrendszerben (racionális törtként felírva): és ; a tényező pedig az és számok valamelyike. Elegendő a -val kiegészített, összesen elemű halmaz pozitív részét megadni (a könnyebb összehasonlíthatóság érdekében mindegyiket nyolcadokban kifejezve):

A példa is mutatja, hogy az azonos kitevőjű szomszédos elemek távolsága , tehát növekedésével exponenciálisan nő ez a távolság: a legnagyobb távolság , a legkisebb pedig .

A mantissza legkisebb értéke (mint már többször is említettük) , a legnagyobb pedig akkor lesz, ha minden számjegye a lehető legnagyobb, azaz . Ha ehhez -t hozzáadunk, könnyen látható, hogy az eredmény éppen lesz, a legnagyobb mantissza tehát . Jelölje és a nem zérus elem abszolút értékének lehetséges legkisebbikét, illetve legnagyobbikát. Előbbiekből következik, hogy

4.2. Kerekítési hibák és becslésük

Legyen nem feltétlenül -beli, de . Az helyett a lebegőpontos aritmetikában az -hez legközelebbi -beli elem ábrázolása történik; jelölje ezt az elemet . Ílymódon az ábrázolást leképezésnek tekinthetjük és kerekítésnek nevezzük. Például ötjegyű decimális aritmetika esetén a számot a számhoz kerekítjük. Ha két szomszédos -beli szám az -től egyenlő távol van, akkor általában a nagyobbik számhoz kerekítünk.

Az alapműveleteket tekintve, legyen és jelölje a négy aritmetikai művelet bármelyikét. A következő esetek lehetségesek:

(pontos eredmény),

(aritmetikai túlcsordulás),

(aritmetikai alulcsordulás),

, (nem ábrázolható eredmény).

Az utolsó két esetben a lebegőpontos aritmetika az eredményhez hozzárendeli a legközelebbi -beli számot.

4.2. Megjegyzés. A fentebb említett eljárást helyes kerekítésnek nevezzük (vannak másfélék is), a helyesen kerekített szám számjegyeit pedig helyes jegyeknek. Az értékes jegyek egyben helyes jegyek is, ezért az a jó gyakorlat, ha csak helyes jegyeket írunk le. Előbbi számról mondhatjuk, hogy körülbelül , de nem mondhatjuk, hogy körülbelül , mert ezzel azt sugallanánk, hogy négy tizedesjegyre helyesen van kerekítve, pedig ez nem igaz.

A kerekítés során keletkezett hiba a kerekítési hiba. A számítások végeredményét nagy mértékben befolyásolják a kerekítési hibák, ezért nagyon fontos az abszolút és a relatív kerekítési hibák nagysága. Mint láttuk korábban, két szomszédos szám távolsága , azaz a közöttük elhelyezkedő számhoz rendelt szám abszolút kerekítési hibakorlátja . Ez bizony igen nagy lehet, de a relatív kerekítési hiba a változó kitevőtől független, a -n kívül csak az aritmetikára nézve konstans -től függ és annak növekedésével exponenciálisan csökken.

4.3. Tétel. Legyen . Az kerekítés relatív hibájára teljesül

A tétel tulajdonképpen azt mondja ki, hogy a lebegőpontos aritmetikában a kerekítés relatív hibája korlátos és ez a korlát , az egységnyi kerekítés mértéke. Az egységnyi kerekítés elnevezést az indokolja, hogy éppen az -hez való kerekítés legnagyobb abszolút (és a alapján egyben relatív) hibája.

Az aritmetika pontosságának jellemzésére az kétszeresét, az értéket szokás használni, amit gépi epszilonnak nevezünk. Az az és a hozzá legközelebbi -nél nagyobb szám távolsága. Bináris alap esetén a következő algoritmussal határozhatjuk meg az értékét

MATLAB program a gépi epszilon meghatározására:

  1    function [epsM]=gepieps() 
  2    x=1; 
  3    while 1+x>1; 
  4       x=x/2; 
  5    end 
  6    epsM=2*x; 
  7    return 

Duplapontosságú szabványos lebegőpontos aritmetikában: .

Két, egyébként -beli és közötti művelet eredményére már nem feltétlenül igaz a 4.3. tétel állítása, hanem csak akkor, ha a kivonás esetén rendelkezésre áll egy tartalék jegy, az ún. ellenőrző jegy. Ez a legtöbb számítógépnél meg is van, de nincs például a CRAY szuperszámítógépeknél, és jó néhány zsebkalkulátornál. Követve a szabványos modellt, mindenesetre feltesszük a lebegőpontos aritmetikai műveletek eredményére vonatkozóan a következőt:

A feltevés fontos következménye, hogy esetén a műveletek relatív hibájára ugyancsak teljesül, hogy

Tehát az aritmetikai műveletek relatív kerekítési hibája kicsi.

Lebegőpontos aritmetikában a műveletekre vonatkozó algebrai azonosságok a kerekítések miatt általában nem állnak fenn.

4.2. Példa. Számoljuk ki a algebrai egyenlőség két oldalát számjegyű, -es számrendszerű aritmetikában.

Megoldás. A műveletekben szereplő valamennyi tag pontosan ábrázolható, azokat tehát külön nem kell kerekíteni. Ugyanakkor , írható tehát, hogy

és

Megjegyezzük, hogy az utóbbi adta a helyesebb eredményt, sőt, az adott aritmetikában attól pontosabb eredmény nem is érhető el. (A duplapontos szabványos aritmetikát megvalósító MATLAB 6.5-ös rendszerben ugyancsak az utóbbi eredmény adódik.)

4.3. Példa. Írjunk MATLAB programot az

összeg kiszámítására a felírt és a fordított, tehát

sorrendben is a természetes

     
    for 
       
    end

rekurzív algoritmussal. Hasonlítsuk össze az eredményeket az esetén.

Megoldás.

MATLAB program a (4.3) és (4.4) számítására:

  1    function [s,sf]=szumma(n) 
  2    s=1; 
  3    for i=1:n 
  4       s=s+1/(i*(i+1)); 
  5    end 
  6    sf=0; 
  7    for i=n:-1:1 
  8       s=s+1/(i*(i+1)); 
  9    end 
 10    sf=sf+1; 
 11    return 

A fentebb említett MATLAB rendszerben a (4.3) formulával (azaz az összegzést a tagok nagyság szerint csökkenő sorrendjében végezve) az , míg a (4.4) összefüggéssel (vagyis növekvő sorrendben összegezve) az eredményt kapjuk. Tehát az utóbbi adta az általános formulára igazolható eredmény -re pontos értékét.

Amikor az összegzést a kisebb tagokkal kezdjük, akkor ezek összegei értékes jegyeket érnek a végső eredményben. Ez a magyarázata annak, hogy az előző két példában a tagok növekvő sorrendben való összegzése adta a jobb eredményt.

4.3. A kerekítési hibák halmozódásának kompenzálása

Nagy mennyiségű, előjelben és nagyságrendben eltérő szám nagy pontosságú összeadása nem egyszerű feladat. Külön válogatva az azonos előjelűeket és azok megfelelő (azaz abszolút értékben növekvő) sorrendű összegzése, majd a két részösszeg kivonása után érhetnénk el a legkedvezőbb eredményt, ez azonban a rendezés miatt nagyon időigényes (költséges) eljárás. Több módszer ismeretes, amelyek megtakarítják a rendezés költségeit, bár nyilván azoknál is növekedik a műveletek száma. Ha nem is adják az optimális végeredményt, de elfogadhatóan közelítik azt. Az egyik legérdekesebb, ilyen célra kifejlesztett eljárás, W. Kahantól származik.

MATLAB program Kahan algoritmusára:

  1    function [s]=kompenzaltsum(x) 
  2    n=length(x); 
  3    s=0; 
  4    kvazi0=0; 
  5    for i=1:n 
  6       temp=s; 
  7       modxi=x(i)+kvazi0; 
  8       s=temp+modxi; 
  9       kvazi0=(temp-s)+modxi; 
 10    end 
 11    return 

Mint látjuk, az algoritmus az algebrailag mindenkor zérus (ténylegesen viszont nem szükségképpen nullához kerekített, azaz aktuálisan csak közel zérus értékűnek ábrázolt) változónak az -hez való hozzáadásával csökkenti a kerekítési hibákat. Megjegyezzük hogy a MATLAB rendszerben a fenti eljárás a (4.3) sorrendű összegzésre is ugyanazt az eredményt adja (vagyis a helyesebbet), mint a (4.4).

4.4. A lebegőpontos aritmetikai szabvány

Az ANSI/IEEE Std 754-1985 bináris () lebegőpontos aritmetikai szabványt 1985-ben hozták nyilvánosságra. A szabvány specifikálja az alapvető lebegőpontos műveleteket, összehasonlításokat, kerekítési módokat, az aritmetikai kivételeket és kezelésüket, valamint a különböző aritmetikai formák közti konverziót. A négyzetgyökvonás az alapvető műveletek közé tartozik. A szabvány nem mond semmit az exponenciális és transzcendens függvényekről.

A szabvány két fő lebegőpontos formátumot ismer: az egyszeres és a dupla pontosságút:

Mindkét formátumban egy bitet az előjelnek tartanak fenn. Minthogy a lebegőpontos számok normalizálva vannak és az első jegy mindig (bináris alap!), ez a jegy nincs tárolva. A mantissza számjegyeinek számában szereplő ezt a rejtett bitet jelzi.

Az IEEE aritmetika zárt rendszer. Minden aritmetikai műveletnek van matematikailag értelmes vagy értelmetlen eredménye. A kivételes műveleteknél jelzést ad, amely után a számításokat előírásszerűen folytatja. Kivételes művelet például a vagy az (ahol véges nem zérus). Előbbi eredménye NaN (Not a Number), utóbbié . Az IEEE aritmetikai szabvány kielégíti a (4.2) modellt.

Egyszeres pontosság esetén a mantissza hossza kb. 7 értékes jegynek felel meg a tizes számrendszerbe átszámolva, míg dupla pontosság esetén kb. 16 értékes jegynek. (Tudniillik és ). Létezik még egy 80 biten ábrázolt, ún. kiterjesztett pontosság is, ahol , a kitevő pedig bites.

5. fejezet - Érzékenység, numerikus stabilitás

Egy feladat érzékenységét az jellemzi, hogy pontos számítások mellett mennyire érzékeny a végeredmény az input adatok megváltozására (perturbációjára). Gyakorlati feladatok input adatai rendszerint mérési eredmények, hibával terheltek; ráadásul ugyanazt többször mérve, más és más eredményt kaphatunk. Az érzékenység fogalma helyett használhatjuk vele ellentétesen a feladat stabilitása fogalmát is. Nyilván a kevésbé érzékeny feladat a stabilabb.

E fogalom egyik jellemzője a kondíciószám, amely a relatív hibákat hasonlítja össze.

5.1. Függvények kondíciószáma

Egy függvényérték kiszámítása rendszerint egy algoritmussal történik (tulajdonképpen már egyetlen aritmetikai művelet is annak tekinthető), ezért érdekes megvizsgálni, hogy az input adat relatív hibáját az algoritmus hányszorosra „nagyítja“ fel, amit a

mennyiség fejez ki. Egyszerű átalakításokkal adódik, hogy

5.1. Definíció. A

mennyiséget az függvény pontbeli kondíciószámának nevezzük.

Egy függvényt numerikusan instabilnak, vagy rosszul kondicionáltnak nevezünk, ha nagy a kondíciószáma. A függvény stabil, vagy jól kondicionált, ha a kondíciószám kicsi. Természetesen a kicsi és nagy jelző relatív. Ezek a relatív jelzők adott feladat esetén a rendelkezésre álló számítógép aritmetikájától és a közelítés megkövetelt pontosságától függenek.

5.1. Példa. Vizsgáljuk az függvényt. Ennek kondíciószáma , amely esetén nagy. Tehát az értékekre a relatív direkt hiba nagy lesz.

5.2. Példa. Az és . Ekkor

ami tetszőlegesen nagy lehet, ha elég közel van -hez. Ezért a példa függvénye numerikusan instabil. Ha bevezetjük az új változót, akkor kapjuk, hogy . Ennek a függvénynek a helyen vett kondíciószáma

Ha , azaz , akkor a kondíciószám kicsi marad. Tehát stabilizáltuk a számítást egy egyszerű átalakítással.

5.3. Példa. Tegyük fel, hogy -t hibakorláttal tudjuk kiszámítani. Melyik kifejezést lehet kisebb relatív hibával kiszámítani és hányszor kisebbel az alábbi, elméletileg egyenlő két kifejezés közül:

(i) ,

(ii) .

Magyarázzuk is meg, miért.

Megoldás. Legyen , , , . , . Az (i) kifejezést tehát kb. 5577-szer kisebb relatív hibával számíthatjuk ki.

Magyarázni azzal is lehet, hogy az függvénynek jóval nagyobb a kondíciószáma az -nél, mint a -nak. Ugyanis , , behelyettesítve az értéket, majd alkalmazva az (5.1) összefüggést, előbbi kb. , utóbbi pedig kb. . A két kondíciószám aránya éppen az előbbi 5577.

A kondíciószámot értelmezhetjük az többváltozós (ún. vektor-vektor) függvényre is (lásd az 5.4. példát). A levezetést mellőzve az (5.1) összefüggésből formálisan is adódik, hogy

ahol az ún. Jacobi-mátrix.

A kondíciószám normafüggő.

5.4. Példa. Mutassuk meg, hogy az kifejezés numerikusan instabil esetén.

Számítsuk ki a kifejezés értékét esetén. Alakítsuk át a kifejezést numerikusan stabillá.

Megoldás. Az és a számítása során keletkezett kerekítési hibák egymástól függetlennek tekinthetők, tehát az függvény viselkedését vizsgálhatjuk az esetén.

, , így euklídeszi normája , és mivel az függvény határértéke az pontban , ezért . Ezeket behelyettesítve (5.2)-be, megkapjuk a kondíciószámot: . Tekintettel arra, hogy , ez a kondíciószám tetszőlegesen nagy lehet, az (és így az eredeti is) tehát igen rosszul kondicionált az adott pontban. Zérushoz közeli -ekre , és ezekkel számolva a kondíciószám pl. -nél kb. . MATLAB-ban -nél (a kb. helyett).

Az azonosság másik oldala viszont numerikusan stabil a megadott környezetben. A kondíciószám meghatározását akár a (ahol ), akár a (ahol ) függvénynél az olvasóra bízzuk.

Az alábbi két ábra jól szemlélteti a viszonyokat.

5.1. ábra - Numerikusan stabil és instabil függvények

Numerikusan stabil és instabil függvények

Kékkel (a baloldali ábrán csak néhány pontot megadva) a helyes függvényt ábrázoltuk, pirossal pedig az eredeti -et. A jobboldalon a környezetében kinagyítottuk az ábrát.

5.2. Direkt és inverz hibák

A függvényértékek számítása során – mint már említettük – hiba következhet be. Jelölje és a pontos értékeket és legyen pontosan , a ténylegesen számított behelyettesítési érték pedig . Az eltérést, azaz a értéket direkt hibának nevezzük. Amennyiben az -ra valamely értékkel pontosan fennáll, hogy , akkor a értéket inverz hibának mondjuk.

5.2. ábra - A direkt és az inverzhiba kapcsolata

A direkt és az inverzhiba kapcsolata

Az ábra a kétféle hibát szemlélteti (, ). Kék színnel itt is a pontos számítást (a folytonos vonal a pontos értékkel, a szaggatott a megváltoztatott inputtal történőt jelöli) mutattuk be, piros pont-vonallal pedig a tényleges számítást (tehát aminek a végén csak valamilyen közelítő érték jelenik meg). Az és az megváltozást (vagy megváltoztatást) perturbációnak is szoktuk említeni. Az inverz hiba elemzését és becslését inverz hibaanalízisnek nevezzük. Ha több inverz hiba is létezik, akkor a (valamilyen normában) legkisebb inverz hiba meghatározása az érdekes. (Gondoljunk például arra, hogy ha és , akkor többféle is szolgáltathatja ugyanazt az eredményt.)

A direkt és az inverz hiba kapcsolatának vizsgálatához tegyük fel, hogy kétszer folytonosan differenciálható. Ekkor tehát felírható a következő Taylor-polinom: , ahol . Így a számított megoldás abszolút hibája

A relatív hiba pedig

Innen kapjuk az (5.1) figyelembe vételével az alábbi, hibaszámítási „ökölszabálynak“ is nevezett

közelítő egyenlőtlenséget, amely szóban kifejezve a következő:

Az egyenlőtlenség azt mutatja, hogy egy rosszul kondicionált probléma számított megoldásának nagy lehet a (relatív) direkt hibája. Egy értéket számító algoritmust direkt stabilnak nevezünk, ha a direkt hiba kicsi és inverz stabilnak nevezzük, ha bármely értékre olyan számított értéket ad, amelyre a inverz hiba kicsi. A „kicsi“ jelző környezetfüggő. Egy direkt stabil módszer nem feltétlenül inverz stabil. Ha az inverz hiba és a kondíciószám kicsi, akkor az algoritmus direkt stabil.

5.2. Megjegyzés. A gyakorlatban természetesen a számítás végeredményének a hibája, a direkt hiba a fontos. Az inverz hibaanalízis jelentősége abban áll, hogy sokszor az inverz hibát tudjuk becsülni. Az alkalmazott számítógép számábrázolási pontossága rendszerint ismert, gyakran annak mérőszámát, a gépi epszilont vagy az azzal arányos mennyiséget tekinthetjük inverz hibának. Az arányossági tényező megállapítása tapasztalatok alapján történik, szakkönyvek is ajánlanak értékeket. Jól kondicionált feladat esetén pedig az inverz hibából következtethetünk a direkt hibára.

Példaként említjük az aritmetikai műveleteket. Az összeadás és a szorzás stabil művelet bármilyen két adattal végezve is. Ugyanakkor a kivonás vagy az osztás bizonyos input adatok mellett érzékeny, instabil. A műveletek érzékenysége kapcsán figyelmükbe ajánljuk az 5.3. alfejezet 14. feladatát.

Rosszul kondicionált (érzékeny, instabil) feladat megoldását általános célú algoritmussal rendszerint nem tudjuk elfogadható hibahatáron belül elérni, az valamilyen speciális eljárást igényel, sokszor a feladat átalakítását. Ilyesmit láttunk a 3.2. vagy az 5.4. példában is. Most egy más jellegű példát mutatunk be.

5.5. Példa. Tekintsük az , , és

sorozatot és vizsgáljuk meg a feladat (pl. az kiszámítása) érzékenységét, ha input adatoknak az és értékét tekintjük (az együtthatókat pontos állandóknak tartjuk).

Megoldás. Az (5.5) képlet egy ún. rekurzív előírás a sorozat elemeire (magát az egyenletet szokás differenciaegyenletnek is nevezni, amelynek az a független változója és annak differenciája ). A rekurzió lineáris, a régóta ismert elmélettel megoldható és a sorozat általános tagjára az explicit összefüggést kapjuk, ahol és tetszőleges konstansok. Figyelembe véve az input adatokat (vagyis az , előírást), a és adódik. Tehát a konkrét feladat megoldása:

Ha a megadott kezdeti értékekhez képest kicsit más input adatokkal oldjuk meg a feladatot (például mérési hibák miatt), akkor a és értéke is kicsit változik, de akármilyen kicsi is a , a zérustól különböző értéke miatt a tag előbb-utóbb hatalmas változást okoz (5.6)-hoz képest. A következő ábrán szemléltetjük az elméleti (5.6) megoldást és mellette az (5.5) direkt alkalmazásával írt MATLAB program eredményeit.

5.3. ábra - Nem érzékeny és érzékeny feladat

Nem érzékeny és érzékeny feladat

Szemléltessük egy ábrasorozaton külön-külön is az eseteket (kékkel itt is a pontosnak tekinthető (5.6), pirossal a rekurzív (5.5) formulával számolt értékek szerepelnek):

5.4. ábra - Kerekítési hibák hatása érzékeny feladatnál

Kerekítési hibák hatása érzékeny feladatnál

Az ábrákból is látható, hogy a rekurzióval számolva -től kezdve az elméleti megoldáshoz képest lényeges és egyre növekvő eltérés jelentkezik. Ráadásul a monotonitás is elvész, sőt a nemnegatívitás is. A sorozat úgy viselkedik, mintha a tag is zérustól különböző együtthatóval szerepelne. A problémát valójában az okozta, hogy a számítások során kerekítési hibák léptek fel. Ha a kerekített eredményt tekintenénk pontosnak, az kerekítések nélkül más, megváltozott (perturbált) input adatokból származna. Jelentkeznek tehát a és a inverz hibák; mintha az input adatok valamely és lennének. Például, ha a pontos , helyett a MATLAB számítás során az (5.5) képlettel kapott (kerekítési hibákkal terhelt) , értékeket vesszük pontos adatoknak, akkor a és adódik; tehát nagyon pici inverz hiba hatalmas direkt hibát okozott. Ugyanakkor az ebből adódó és értékek alig változtak az eredeti , illetve értékekhez képest.

Itt meg is érkeztünk egy másik problémához, az algoritmusok stabilitásához (itt is használatos a jellemzésre az ellentett fogalom, az érzékenység). Az algoritmusoknak a kerekítési hibákra való érzékenységét (vagy éppenséggel kevésbé érzékenységét, robusztusságát) értjük rajta. Éppen ezért szokták a numerikus jelzővel összekapcsolni és numerikus stabilitásnak nevezni. Használatos itt is a jól (vagy gyengén) kondicionált algoritmus kifejezés is. Kapcsolat van a feladat érzékenységével (előbb láttuk hogy az instabil feladat végrehajtása közben a kerekítések hatása olyan, mint az input adatok perturbációja), de valójában más fogalom. Említettük, hogy az összeadás stabil művelet: akárhány azonos tag összeadásának elméleti eredménye relatíve nem változik jelentősen, ha az összeadandókba kis relatív hibák kerülnek. Ugyanakkor, ha azt az algoritmust találnánk ki, hogy olyan sorrendben adjuk össze az adatokat, hogy a meglévő részösszeghez a hátralévők közül – mindaddig, amíg lehetséges – azt adjuk, amelyik egyirányú kerekítést (például mindig lefelé kerekítést) okoz – és lehetőleg minél nagyobbat –, akkor sok szám összeadása után valószínűleg csodálkoznánk a végeredményen. Ez az algoritmus instabil. Ellentétben a Kahán-algoritmussal, amelyik éppen kompenzálni akarja a kerekítési hibákat, és numerikusan stabilnak mondható.

Fontos szabály, hogy a gyakorlatban csak stabil (jól kondicionált) algoritmusokat használunk. Instabil (inkorrekt kitűzésű), vagy rosszul kondicionált feladatot általános célú algoritmusokkal általában nem tudunk megoldani.

5.3. Feladatok

1. A fény sebessége légüres térben km/s. Mekkora a relatív hibakorlátja? Miért helytelen, ha a fény sebességét légüres térben durván km/s-nak adjuk meg? Hogyan kell helyesen megadni?

2. A normálfogazású fogaskerék modulusza: , ahol az osztást jelenti mm-ben (két szomszédos fog távolsága az ún. osztókörön mérve). A moduluszok szabványosítva vannak. Az osztókör ívén értékét mm-nek mértük. Mekkora a relatív hiba, ha tudjuk, hogy a modulusz pontosan ?

3. Egy gömbalakú kondenzátor kapacitása

ahol és a fegyverzetek külső és belső sugara, pedig az anyag dielektromos állandója, melyet pontos értéknek tekintünk.

A mérési eredmények: cm, cm. Számítsuk ki relatív hibáját.

4. Az cm hosszú inga lengésidejét mp-nek mértük. Határozzuk meg a

képletből közelítő értékét a mérés helyén. Adjunk becslést a közelítés abszolút és relatív hibájára a közelítésének abszolút hibája függvényében. Konkrét hibakorlátot is adjunk, ha -t 6 tizedesjegy pontossággal ismerjük.

5. Az és pontok magasság különbségét műszerrel állapítottuk meg. A mért adatokat táblázatba foglaltuk:

A keresett szintkülönbséget a könnyen belátható geometriai összefüggésből számítjuk. Adjuk meg közelítő értékét és becsüljük meg a hibáját, ha értékét 4 tizedesjegy pontossággal ismerjük.

6. A erővel húzott hosszúságú és átmérőjű acéldrót megnyúlása

ahol a -t és az rugalmassági modulust is magába foglaló konstans helyes jeggyel megadott értéke . Feltéve, hogy a erő pontosan , mekkora közös abszolút hibakorláttal kell megadnunk a kb. m nagyságú és a kb.  mm nagyságú értékét, hogy a képletből a megnyúlás  mm pontossággal meghatározható legyen?

7. Az háromszögről a következő adatokat ismerjük:

 m

 m

Legfeljebb mekkora lehet abszolút hibája, ha a -t a koszinusz-tétel alkalmazásával

(i)  m,

(ii)  m

hibakorláttal akarjuk kiszámítani? Adjunk magyarázatot a két eredmény jelentős eltérésére.

8. Egy bója függőleges állásában magasan van a vízszint felett. A felfüggesztő zsinór a vízszintet az pontban döfi. A feltámadó szél megdönti úgy, hogy a felfüggesztő zsinór továbbra is egyenes marad, de a bója éppen a víz felszínére kerül, az ponttól távolságra. Állapítsuk meg a víz mélységét az adatokból és becsüljük meg az eredmény hibáját. (.)

9. A földi gravitációs gyorsulás (a tengerszint feletti magasság függvényében)

ahol a gyorsulás magasságban, a gyorsulás a tenger szintjén, pedig a föld középsugara. Ha , akkor használhatjuk a

közelítő képletet is. Adjuk meg függvényében a abszolút és relatív hibáját.

10. Számítsuk ki az oldalú háromszög területét kétféleképpen:

(i) , ahol ,

(ii) .

Melyik képlet ad pontosabb eredményt , esetén és miért?

11. Vizsgáljuk meg kiszámításának alábbi módjait:

Legyen és , valamint és a hozzájuk tartozó relatív hibakorlátok. Adjuk meg az és relatív hibakorlátait és függvényében. Melyik eljárás mikor jó?

12. Ha elég kicsi, akkor az alábbi (i)-(iii) feladatokban alkalmazhatjuk az ott megadott közelítéseket. Igazoljuk, hogy a közelítő összefüggésből kapott értékek relatív hibakorlátjai az egyes feladatokban százalékosan megadott értékek.

(i) , a relatív hibakorlát:

(ii) , a relatív hibakorlát:

(iii) , a relatív hibakorlát:

13. A

konvergens sor összege . Ha csak az első tagot (de azokat pontosan) vesszük figyelembe, akkor a maradék sor elhagyásából származó hiba kisebb, mint .

Milyen pontossággal kell az egyes tagokat kiszámítani, hogy ezen tizedik részletösszeggel való közelítés hibakorlátja legyen?

14. Tekintsük az alapműveleteket kétváltozós függvényeknek, azaz , ahol a műveletek valamelyike.

(a) Vezessük le az alapműveletek hibakorlátjait, mint kétváltozós függvényekét.

(b) Adjuk meg ezen függvények kondíciószámát. Mikor rosszul kondicionáltak ezek a függvények?

(c) Vezessünk le hibakorlátot a hatványozásra; úgy is, hogy a kitevőt pontos értéknek tekintjük és úgy is, hogy mind az alap, mind a kitevő hibával terhelt.

(d) Legyen és . Adjuk meg függvényében az maximális hibakorlátját úgy, hogy az közelítéssel számolva relatív hibakorlátja legfeljebb legyen.

15. Számítsuk ki az rekurzió első 50 tagját, ha és . Ábrázoljuk a értékeket és vessük össze az eredményt az elméleti megoldással.

16. Számítsuk ki az rekurzió első 70 tagját, ha , és . Ábrázoljuk a értékeket és vessük össze az eredményt az elméleti megoldással. Stabil-e a rekurzió a kezdeti értékekre nézve?

17. Vizsgáljuk meg kiszámításának alábbi két módját:

(i) ,

(ii) .

A közbülső és eredmények relatív hibával vannak kerekítve. Melyik változat a jobbik?

18. Számítsuk ki kétféleképpen az függvény értékét a intervallumon az pontokban :

(i) ,

(ii) .

Ábrázoljuk az eredményeket. Mi az eredmények magyarázata?

19. Számítsuk ki az függvény értékét az intervallum darab ekvidisztans pontjában háromféleképpen is:

(i) a binomiális tétel kifejtésével kapott

polinom tagjait összegezve,

(ii) a polinomot a Horner séma szerint számolva:

(iii) közvetlenül hatványozva: .

Ábrázoljuk az eredményeket. Mi az eredmények magyarázata?

20. Ismert, hogy nemnegatív és mellett a számtani közép nem kisebb a mértaninál, azaz

és az egyenlőség csak esetén áll fenn. Így van-e ez numerikusan is? Ellenőrizzük kísérletileg az állítást, ha és nagyságrendje azonosan igen kicsiny, illetve igen nagy, valamint ha jelentősen eltérő.

21. Mutassuk meg, hogy a következő kifejezések numerikusan instabilak esetén:

(i) ,

(ii) ,

(iii) ,

(iv) .

22. Legyen lebegőpontos szám az IEEE duplapontos rendszerben és teljesüljön, hogy . Mutassuk meg, hogy értéke vagy vagy ( a gépi epszilon).

23. Legyen . Ezt az értéket az esetleges túlcsordulást elkerülendő a formulával számíthatjuk, ahol . Igazoljuk, hogy az alulcsordulás nem okoz számottevő hibát. (Másképpen: ha számítása során alulcsordulás lép is fel, az eredmény nem változik lényegesen ahhoz képest, mintha a korrekt értékkel számolnánk és a végén kerekítenénk.) Igaz-e ez az állítás akkor is, ha ?

24. Melyik állítás és mikor igaz az alábbiak közül:

(i) ; ,

(ii) ,

(iii) ,

(iv) ,

(v) ,

(vi) , ha .

6. fejezet - Lineáris egyenletrendszerek megoldása

A lineáris egyenletrendszerek általános alakja egyenlet és ismeretlen esetén:

Az egyenletrendszert megadhatjuk a tömörebb

formában is, ahol

Ha , akkor alulhatározott, ha , akkor túlhatározott egyenletrendszerről beszélünk. Az esetben pedig az egyenletrendszert négyzetesnek nevezzük. Az egyenletrendszerek geometriai tartalmát a következőképpen írhatjuk le.

6.1. Definíció. Az euklideszi tér normálvektorú és ponton átmenő hipersíkját az

egyenletet kielégítő pontok határozzák meg.

A hipersík egyenlete más formában kifejezve:

Felhasználva az mátrix sorok szerinti

particionálását, ahol , az egyenletrendszert felírhatjuk az ekvivalens

alakban. Innen jól láthatjuk, hogy a lineáris egyenletrendszer megoldása hipersík közös része. Ennek megfelelően három eset lehetséges:

    (i) az egyenletrendszernek nincs megoldása,

    (ii) az egyenletrendszernek pontosan egy megoldása van,

    (iii) az egyenletrendszernek végtelen sok megoldása van.

6.2. Definíció. Ha az lineáris egyenletrendszernek legalább egy megoldása van, akkor az egyenletrendszert konzisztensnek nevezzük. Ha az egyenletrendszernek nincs megoldása, akkor az egyenletrendszer inkonzisztens.

Például a , egyenletrendszer inkonzisztens.

Az egyenletrendszert felírhatjuk az ekvivalens

alakban is, ahol az mátrix -edik oszlopát jelöli ( pedig skalár: a megoldásvektor -edik komponense). A összeg az vektorok lineáris kombinációja (lásd a 2.20. definíciót). Az egyenletrendszer akkor és csak akkor oldható meg, ha kifejezhető az oszlopvektorainak lineáris kombinációjaként.

A megoldhatóságot megállapíthatjuk a rang fogalmának (lásd a 2.22. definíciót) segítségével is:

az egyenletrendszernek akkor és csak akkor van megoldása, ha . Ha , akkor az egyenletrendszernek pontosan egy megoldása van.

A továbbiakban csak négyzetes egyenletrendszerekkel foglalkozunk. Feltesszük tehát, hogy . Ismert a következő tétel.

6.3. Tétel. Az egyenletrendszernek akkor és csak akkor van pontosan egy megoldása, ha létezik . Ekkor a megoldás .

Fenti tételnek egyenes következménye az alábbi állítás, amelyet fontossága miatt külön tételben is meg szoktak fogalmazni.

6.4. Tétel. Az homogén lineáris egyenletrendszernek akkor és csak akkor van nemtriviális megoldása, ha .

Az lineáris egyenletrendszer megoldására jól ismert az ún. Cramer-szabály, mely szerint , ahol (oszlopok szerint particionált formában felírva): . (Az tehát -ból úgy származtatható, hogy az -edik oszlopát kicseréljük a jobboldali vektorral.) Az darab determináns kiszámítása viszont rendkívül számításigényes, még akkor is, ha a determinánsok kiszámítására nem a kifejtési szabályt (lásd a 2.24. definíciót), hanem a determináns tulajdonságait kihasználó, más, kisebb műveletszámú algoritmust használunk.

Kézenfekvő az – egyébként a középiskolában is tanított – következő eljárás: „egyik“ egyenletből fejezzük ki explicite „valamelyik“ ismeretlent; az arra kapott összefüggést helyettesítsük be a maradék többibe, így azok eggyel kevesebb ismeretlent tartalmazó egyenletrendszert alkotnak. Erre az () egyenletre ismételjük meg az eljárást, majd újra, egész addig, amíg csupán egyetlen ismeretlen marad. Ezt nevezhetjük első fázisnak. A folytatásban kiszámoljuk az említett egyetlen ismeretlent, ezután a már ismertté vált értékét visszahelyettesítve a többibe, újabb ismeretlent számolhatunk ki, és így tovább. A számítások ezen részét nevezhetjük második fázisnak.

Végeredményben el is mondtuk a – specialitást nem mutató – egyenletrendszerek megoldására ma is a legáltalánosabban alkalmazott Gauss-módszert. Csupán az „egyik“ , illetve „valamelyik“ kiválasztására kell számítógépre programozható és a kerekítési hibákat kordában tartó szabályt felállítani.

Mire eljutunk az egyetlen ismeretlent tartalmazó egyenletig, akkorra az eredeti egyenletrendszert egy vele ekvivalens (a megoldást őrző) háromszögmátrixú egyenletrendszerré alakítottuk. Ezért először az ilyen egyenletrendszerek megoldásával foglalkozunk.

6.1. Direkt módszerek

6.1.1. Háromszögmátrixú egyenletrendszerek

Tekintsük az egyenletrendszert, ahol felső háromszögmátrix (lásd a 2.14. definíciót). Az egyenletrendszert részletesen is írjuk ki:

Az egyenletrendszer akkor és csak akkor oldható meg egyértelműen, ha , azaz (lásd a 2.7. alfejezet 3. feladatát). Ezen egyenletrendszer megoldását adja a (szavakban végül is már elmondott) következő algoritmus:

Az kiszámítását végző sorban szereplő szumma tulajdonképpen két vektor (az mátrix és az vektor egy-egy partíciója) skaláris szorzata, így felírhatjuk a a részmátrixokra bevezetett jelöléssel (lásd a 2.3. alfejezetet) is:

persze azzal a feltétellel, hogy az oszlopvektor. Ezt kihasználva írhatjuk meg a megfelelő MATLAB programot.

MATLAB program felső háromszögmátrix esetén:

  1    function [x]=fhmatrixszal(A,b) 
  2    n=length(A); 
  3    x=zeros(n,1); 
  4    for i=n:-1:1 
  5       x(i)=(b(i)-A(i,i+1:n)*x(i+1:n))/A(i,i); 
  6    end 
  7    return 

A 3. sor tulajdonképpen méret deklaráció (értékadásként felesleges lenne); azért volt szükség rá, hogy az vektor oszlopvektor legyen, mert csak így értelmezhető az 5. sorban lévő partíciók szorzása. Kihasználva, hogy az üres vektor, ha , az kiszámítására sem kellett külön utasítást adni, az beépíthető a ciklusba.

Az

alsó háromszögmátrixú egyenletrendszert megoldó algoritmus pedig értelemszerűen:

A MATLAB program megírását gyakorlásképpen az olvasóra bízzuk. (Később látni fogjuk, hogy külön programra valójában nincs is szükség; a fenti program alkalmas hívásával alsó háromszögmátrixú egyenletrendszert is meg lehet oldani.)

6.1.2. A Gauss-módszer

A Gauss-módszer, mint már elmondtuk, két fázisból áll.

I. Azonos (a megoldást őrző) átalakításokkal az egyenletrendszert felső háromszög alakúra hozzuk:

6.1. ábra - A Gauss elimináció sémája

A Gauss elimináció sémája

II. A kapott felső háromszögmátrixú egyenletrendszert a (6.5) algoritmussal megoldjuk.

Az azonos átalakítást itt most úgy végezzük el, hogy egyik egyenletből kivonjuk a másik egyenlet alkalmasan megállapított konstansszorosát, így ott a szóban forgó ismeretlen együtthatója zérussá válik. Kiiktatjuk – idegen szóval: elimináljuk – az ismeretlent, ezért is nevezzük a módszert eliminációs eljárásnak. Az első lépésben az

alakot kell tehát kapnunk, amit ha (6.1)-ben , akkor elérhetünk úgy, hogy az -edik sorból () kivonjuk az első sor -szeresét:

Az feltételből kapjuk, hogy . Ha ezt a értéket behelyettesítjük a (6.7) egyenletbe, akkor könnyen ellenőrizhető, hogy itt tényleg arról van szó, hogy az első egyenletből kifejezzük az -et, és a kapott kifejezést behelyettesítjük a maradék többibe. A számítások során nem kell az szimbólumokat magunkkal cipelni, elég, ha az együtthatómátrix elemein hajtjuk végre a megfelelő módosítást. Így könnyen programozható a kinullázás alábbi algoritmusa:

Az algoritmus felülírja az mátrix indexű és a vektor indexű elemeit. A futhatna éppen -től is, de felesleges az első oszlop alatti elemeit számolni, hisz tudjuk, hogy azok zérusak lesznek. Beírni is felesleges ezeket a -kat, hiszen láttuk, hogy a visszahelyettesítő (6.5) algoritmus nem hivatkozik rájuk (illetve eleve a egyenletrendszert oldja meg). Az első egyenletet nem tekintve, most már egy ismeretlent tartalmazó kisebb, -es egyenletrendszerhez jutottunk. Erre alkalmazhatjuk az előbb ismertetett kiiktatást, majd folytathatjuk ugyanígy a még kisebb méretű egyenletrendszeren. Tegyük fel, hogy az első oszlopban a kinullázást már elvégeztük és az

egyenletrendszert kaptuk. (Itt megjegyezzük, hogy az együtthatók már nem az eredeti mátrix értékeit tartalmazzák, hiszen azok az előző lépésekben módosultak. Valójában egy harmadik – például felső – indexszel kellene ellátnunk őket, de ez itt inkább csak bonyolítana; nélküle is érthetők a viszonyok.) Ha , akkor kinullázzuk -nak az alatti együtthatóit, azaz az -edik sorból a -adik sor -szeresét kivonjuk:

Az feltételből kapjuk, hogy . A -adik oszlop alatti kinullázását tehát a

algoritmussal végezhetjük el. Ez szóról-szóra azt jelenti, hogy a -adik egyenletből kifejezzük a -adik ismeretlent és behelyettesítjük a hátralévő egyenletekbe, csak éppen a nullák helyén fizikailag otthagyjuk a régi számot.

A kinullázást mindaddig folytathatjuk, amíg az ; egészen a -ig (az -edik egyenletben már nincs a főátló alatt elem). Ha sikerül az mátrixot felső háromszög alakra hozni (úgy, hogy a végén is fennáll), akkor a (6.5) algoritmus alkalmazásával eljutunk a megoldáshoz.

A főátlóbeli elemeknek kulcsszerepük van az eliminációs folyamatban: pivot elemnek nevezzük az éppen aktuálisat. A megállapítására egyszerűen megjegyezhető verbális szabályt adhatunk: egy tört, melynek számlálója a kinullázandó elem, nevezője pedig a pivot elem. Tehát az algoritmus:

A Gauss-módszer algoritmusa

I. (eliminációs) fázis:

II. (visszahelyettesítő) fázis:

Az algoritmust végrehajtó MATLAB program megírásánál itt is felhasználjuk a partícionálás lehetőségét, valamint azt, hogy az -edik sor módosításánál kihagyhatjuk a sor első elemét (hiszen azok már zérusok és azok is maradnak). Sőt az -edik módosításától is eltekinthetünk, mivel a II. fázis egyáltalán nem hivatkozik a főátló alatti elemekre, azok konkrét aktuális értéke közömbös, nem kötelező, hogy zérus legyen, ott maradhat a megelőző. Azaz csupán az partícióból kell a -adik sor megfelelő partíciójának -szorosát kivonni.

MATLAB program a Gauss-módszerre:

  1    function [x]=Gauss(A,b) 
  2    n=length(b); 
  3    % 
  4    %eliminációs fázis: 
  5    % 
  6    for k=1:n-1 
  7       for i=k+1:n 
  8          g=A(i,k)/A(k,k); 
  9          A(i,k+1:n)=A(i,k+1:n)-g*A(k,k+1:n); 
 10          b(i)=b(i)-g*b(k); 
 11       end%i 
 12    end%k; 
 13    % 
 14    %vége az eliminációs fázisnak 
 15    % 
 16    [x]=fhmatrixszal(A,b) 
 17    return 

Itt a 18. sorban a (6.5) algoritmus MATLAB programját hívjuk.

6.1. Példa. Oldjuk meg a következő egyenletrendszert.

Megoldás. Hajtsuk végre szó szerint a fenti eljárást, de zárójelbe téve azért jelezzük az elvileg -kat.

Determináns kiszámítása Gauss-módszerrel. Ismeretes, hogy a determináns értéke nem változik, ha bármely sorához hozzáadjuk egy másik sor akárhányszorosát. Mivel a fenti eljárás közben csak ilyen átalakításokat hajtunk végre a mátrixon (és így a determinánsán is), a determináns marad az eredeti. Azt tehát a végén a főátlóbeli elemek szorzata adja (lásd a 2.7. alfejezet 3. (ii) feladatát). Előbbi példánkban tehát: .

6.1.3. A Gauss-módszer műveletigénye

A Gauss-módszer véges sok lépésben, véges sok aritmetikai alapművelet elvégzése után megadja az () egyenletrendszer megoldását. A szükséges aritmetikai műveletszám (műveletigény) az egyenletrendszer-megoldó eljárások fontos minőségi jellemzője, mert az ilyen algoritmusok számítógépideje nagyjából arányos az aritmetikai műveletigénnyel (a kerekítési hibák halmozódása is – bár ez konkrét esetekben nem törvényszerű – kedvezőtlenebb lehet nagyobb műveletszám esetén). A műveletszámot az elegendő információt szolgáltató nagyságrendjével szoktuk megadni. Az itt használatos nagyságrend-fogalom a következő:

6.5. Definíció. Az mennyiség nagyságrendű, ha van olyan szám, hogy ().

6.6. Megjegyzés. A műveletigényt szokás az algoritmus költségének is nevezni, ennek megfelelően beszélhetünk olcsó, illetve drága eljárásokról.

Ha már a fejezet elején a Cramer-szabály nagy műveletigényét említettük, számláljuk össze a Gauss-módszer által igényelt additív és multiplikatív műveleteket. Jelöljük az additív műveleteket -val, a multiplikatívakat pedig -mel. (A hagyományos felépítésű számítógépeknél az összeadás és a kivonás nagyjából egyező időigényű, tőlük nagyobb, de egymás között szintén azonosnak vehető a szorzás és az osztás időigénye; indokolt tehát a műveletek ilyen szempontból két csoportba való sorolása.)

Az I. fázis -adik lépésének utasításai, valamint a bennük szereplő additív (azaz: ) és multiplikatív (azaz: ) műveletek és a műveletszámok (illetve költségek) a következők:

Fentieket összeadva és -val szorozva (hisz a ciklusmag -szor fut le), kapjuk:

Minthogy ez a belső ( változójú) ciklus a külső ciklus lépéseire fut le, ki kell számolni a

összegeket. Ismert azonosságok felhasználásával adódik az

eredmény. A II. fázis műveletigénye hasonló számítással:

Összegezve az lépések műveletigényét, az helyett a változó bevezetésével kapjuk, hogy a II. fázis összköltsége

Az I. és II. fázis költségét összeadva kapjuk a Gauss-módszer számítási összköltségét:

Említettük már, hogy aritmetikai műveleteket használó akármely véges algoritmus lényeges jellemzője a végrehajtásához szükséges aritmetikai műveletek száma, természetesen a feladat paramétereinek (pl. ismeretlenek száma, együttható mátrix mérete, stb.) függvényében. A régebbi számítógépek esetén a multiplikatív műveletek végrehajtási ideje lényegesen nagyobb volt, mint az additívaké. Ezért ezeket külön határozták meg. Később megfigyelték, hogy a lineáris algebra számítási eljárásaiban az additív és a multiplikatív műveletek száma nagyon gyakran közel azonos. Ezért alkották meg a számítási igény mérésére a flop fogalmát.

6.7. Definíció. 1 régi flop az a számítási munka, amely az művelet (1 összeadás + 1 szorzás) elvégzéséhez kell. 1 új flop pedig az a számítási munka, amely egyetlen (mindegy, hogy additív vagy multiplikatív) aritmetikai művelet elvégzéséhez szükséges.

Az új flop bevezetését az indokolja, hogy az újabb számítógépeken a multiplikatív és az additív műveletek ideje azonosnak tekinthető. Tehát egy régi flop 2 új floppal azonos a mai számítógépeken.

Például az skalárszorzat () kiszámítása az

algoritmussal régi flop, de új flop műveletigényű.

Figyelembe véve mindezeket, valamint, hogy a (6.8) formulában nagy -ekre a legmagasabb fokú tag válik dominánssá, kimondhatjuk a Gauss-módszer műveletigényéről szóló tételt.

6.8. Tétel. A Gauss-módszer additív és ugyanennyi multiplikatív műveletet igényel, azaz összességében a műveletigénye régi és új flop.

Igazolni lehet, hogy ha csak sor- és oszlopműveleteket (sor vagy oszlop számmal való szorzása; sorok vagy oszlopok cseréje, valamint sor vagy oszlop számszorosának másik sorhoz vagy oszlophoz való hozzáadása) engedünk meg, akkor nem is lehet régi flopnál kevesebb művelettel az lineáris egyenletrendszert megoldani.

Kutatások folynak ezen műveletszám lejjebb szorítására, de jelenleg nem ismerünk numerikusan stabil, a Gauss-módszernél olcsóbb eljárásokat.

6.1.4. A főelemkiválasztásos Gauss-módszer

A Gauss-módszer I. fázisában előfordulhat, mondjuk a -adik lépésben, hogy az elem zérus. Például a

rendszernél . Ilyen esetekben a sorok, vagy az oszlopok felcserélésével megkísérelhetjük elérni, hogy az helyére zérustól különböző elem kerüljön. A fenti esetben például az első és harmadik sor felcserélésével kapjuk, hogy

Az első és második oszlop cseréjével pedig azt kapjuk, hogy

A sorok cseréjénél az egyenletek sorrendje, az oszlopok cseréjénél pedig a változók sorrendje változik meg. Általában, így az előző példában is, több választási lehetőségünk is van sor-, vagy oszlopcserére. Ha azonban az elem alatt a többi együttható is zérus, akkor az részmátrix oszlopai lineárisan összefüggők, szinguláris és az eliminációs eljárás sorcserével sem folytatható. Hasonló a helyzet, ha is és a sorában, tőle jobbra, minden együttható is zérus: ekkor is szinguláris.

A sorok felcserélését úgy, hogy az új pivot elem zérustól különböző legyen, pivotálásnak vagy főelemkiválasztásnak nevezzük. A pivot elem megválasztása nagymértékben befolyásolja az eredmények numerikus stabilitását (lásd a 7.2. alfejezetet).

6.2. Példa. Oldjuk meg a következő egyenletrendszert -es alapú, jegyű aritmetikában Gauss-módszerrel, pivotálás nélkül és pivotálással is. (Mint már korábban említettük, ez az aritmetika a pontosságát tekintve az esetek nagy részében jól modellezi a szabványos lebegőpontos aritmetikát.)

Megoldás. Mindenek előtt megemlítjük, hogy az elméleti megoldás: , . Mindkettő ábrázolásához -nál több jegyre van szükség, előbbi egy nagyon kicsit nagyobb -nél, utóbbi ugyancsak nagyon kicsivel kisebb: mindkettő kerekített értéke . Tehát az értéknél jobb közelítésű megoldás az adott aritmetikában nem is lehetséges.

Pivotálás nélkül az elimináció után kapott háromszögmátrixú egyenletrendszer

A szükségszerű kerekítés utáni, ténylegesen ábrázolt együtthatójú egyenletrendszer pedig

Ennek megoldása az elméleti megoldáshoz képest katasztrófálisan rossz .

Ha az elimináció előtt főelemkiválasztást alkalmazunk (pivotálálunk), akkor a háromszögmátrixú egyenletrendszer

Itt a nyíl után látható a ténylegesen ábrázolt (kerekített) együtthatójú egyenletrendszer. Ennek megoldása: ; azaz nagyon jó, sőt, az elérhető legjobb közelítése az elméleti eredménynek.

Megjegyezzük, hogy valóságos számítógépen, MATLAB rendszerben ugyanezt tapasztalhatjuk.

Általában is igaz, hogy a közelítő megoldás pontosságát nagymértékben javítja a helyesen megválasztott pivotálás. Pivot elemnek nagy abszolút értékű elemet kell választani. Két alapvető pivotálási stratégiát használunk.

Részleges főelemkiválasztás: A -adik lépésben a -adik oszlop () elemei közül kiválasztjuk a maximális abszolút értékűt. Ha ennek indexe , akkor a -adik és a -edik sort felcseréljük. A pivotálás után teljesül, hogy

A részleges elnevezést az indokolja, hogy csak az aktuális oszlopbeli, szóba jöhető elemek közül (azaz az partícióból) választjuk ki a legnagyobb abszolút értékűt.

Teljes főelemkiválasztás: A -adik lépésben az () mátrixelemek közül (azaz az partícióból) választjuk ki a maximális abszolút értékűt. Ha ennek indexpárja , akkor a -adik és a -edik sort, valamint a -adik és az -edik oszlopot felcseréljük. A pivotálás után teljesül, hogy

Ismét felhívjuk a figyelmet arra, hogy oszlopcsere esetén változócsere is történik.

6.9. Megjegyzés. Ha determinánst számolunk pivotálással, akkor figyelemmel kell követnünk a sor- és az oszlopcserék számát, hiszen egy sor- vagy oszlopcsere a determináns előjelét megváltoztatja. Tehát a végső háromszögmátrix főátlóbeli elemeinek szorzatát -gyel szorozni kell, amennyiben páratlan számú cserét hajtottunk végre.

A főelemkiválasztásos Gauss-módszer esetén az I. fázis minden lépésében pivotálást hajtunk végre. A teljes főelemkiválasztást biztonságos (numerikusan stabil) stratégiának tekinthetjük. A részleges főelemkiválasztás egyéb technikákkal kiegészítve ugyancsak biztonságosnak tekinthető. Éppen ezért a teljes főelemkiválasztás árát nem éri meg megfizetni. (A nagyság szerinti összehasonlítás is műveletet igényel; rendszerint különbségképzést és az eredmény előjel-bitjének figyelését. Ehhez jön még a többszöri abszolút érték-képzés és a nagyobb adatmozgatás.) A részleges főelemkiválasztás lényegesen kevesebb pótlólagos aritmetikai műveletet igényel mint a teljes főelemkiválasztás. A gyakorlatban tehát legtöbbször csak részleges főelemkiválasztást alkalmazunk, de azt – később említett kivételektől eltekintve – mindig!

Figyeljük meg, hogy a jobboldali vektor elemeit ugyanúgy kell transzformálni az elimináció során, mint az elemeit. Megtehetjük tehát, hogy kibővítjük az mátrixot -edik oszlopként a vektorral és a rögzített soron („vízszintesen“ ) futó ciklusváltozót nem -ig, hanem -ig futtatjuk (azaz úgy, hogy a partíció utolsó oszlopindexe legyen):

MATLAB program részleges főelemkiválasztással:

  1    function [x]=Gauss_rfoelem(A,b) 
  2    n=length(b); 
  3    A=[A,b]; 
  4    %eliminációs fázis: 
  5    % 
  6    for k=1:n-1 
  7       for i=k+1:n 
  8          [amax,ind]=max(abs(A(k:n,k))); 
  9          t=ind+k-1; 
 10          %ind az A(k:n,k) particiobeli index, t pedig az A-beli 
 11          if t  k 
 12             %a k-adik es a t-edik sor csereje: 
 13             A([k,t], :)=A([t,k], :); 
 14          end % if 
 15          g=A(i,k)/A(k,k); 
 16          A(i,k+1:n+1)=A(i,k+1:n+1)-g*A(k,k+1:n+1); 
 17       end %i 
 18    end %k; 
 19    % 
 20    %vége az eliminációs fázisnak 
 21    A (6.5) algoritmus MATLAB programja: 
 22    [x]=fhmatrixszal(A(:, 1:n),A(:, n+1)) 
 23    return 

6.10. Megjegyzés. Ez a szokásos programozási eljárás. Ha darab, csupán a jobboldali vektorban különböző egyenletrendszert kell megoldanunk, akkor mindegyikkel kibővíthetjük az mátrixot, és a soron végigfutó ciklust -ig futtatva egyszerre történik mindegyiknél az I. fázis. (A II. fázist is lehet egyszerre végezni.) A kibővítés történhet az utasítással vagy az módon, ahol a vektorok a mátrix oszlopvektorai. Konkrét kivitelezést láthatunk majd a Gauss-Jordan eliminációnál. A MATLAB-ban egyébként az utasítás algebrailag az eredményt szolgáltatja (amennyiben ez létezik).

A fenti programmal kapcsolatban felhívjuk a figyelmet, hogy a már „kinullázódott“ elemeket nem kellene cserélni, a 13. sor lehetne A([k,t], k:n+1)=A([t,k], k:n+1); is. Később látni fogjuk, hogy ugyanez az algoritmus használható az ún. -felbontásnál, ott viszont a fizikailag tárolt elemekre szükség van. Ám sávmátrixok esetén kihasználhatjuk, hogy vannak biztosan zérus és annak is maradó elemek, így a műveletigényt jelentősen lecsökkenthetjük. A gyakorlatban jelentős szerepet játszó, ún. tridiagonális mátrixok () esetén ez a műveletigény mindössze flop. Arra viszont ügyelni kell, hogy a felső sávszélesség megnőhet -ra. A 6.2. ábra éppen azt az esetet szemlélteti, amikor a -adik sort a -edik sorral kell felcserélni:

6.2. ábra - Felső sávszélesség növekedése elimináció során

Felső sávszélesség növekedése elimináció során

Megemlítjük még a gyakorlatban többször előforduló, azon eseteket, amikor nem kell pivotálni: az eljárás főelemkiválasztás nélkül is stabilan viselkedik. Ezek a következők:

Az mátrix diagonálisan domináns a következő értelemben:

Az mátrix szimmetrikus és pozitív definit.

A pozitív definitség fogalma a következő:

6.11. Definíció. Az mátrix pozitív definit, ha , esetén .

Szimmetrikus és pozitív definit mátrix esetén a Gauss-eliminiáció egy speciális alakját, a Cholesky-módszert használjuk az egyenletrendszer megoldására (lásd később). Itt utalunk a 8. fejezetben megemlített MATLAB műveletre, amely az megoldását a Help szerint a következőképpen adja: gyorsteszttel megvizsgálja, hogy szimmetrikus, pozitív definit-e, és ha igen, akkor Cholesky-módszert alkalmaz, egyébként pedig részleges főelemkiválasztásos Gauss-módszert (bár, meglehet, ezután még valamilyen javító lépést is alkalmaz; erről hallgat a MATLAB Help).

A numerikus stabilitást nagymértékben befolyásolja a pivotelemek növekedési tényezője. Ennek definiálása érdekében itt mégiscsak jelezzük egy felső indexszel, hogy az eliminációs lépések során tulajdonképpen újabb és újabb együttható mátrixok keletkeznek; a régi értékek transzformálódnak.

Legyen és . A Gauss-módszer eliminációs fázisában egy

ekvivalens egyenletrendszerekből álló sorozatot képezünk, ahol

A végső állapotban

ahol a -adik fő-, vagy pivotelem. A pivotelemek növekedési tényezője:

Ez a növekedési tényező összefüggésben áll az eljárás numerikus stabilitásával.

Itt kell megjegyezni, hogy újabban találtak a gyakorlatban (differenciál- és integrálegyenletek közelítő megoldásánál) előforduló egyenletrendszereket, melyeknél a részleges főelemkiválasztáson alapuló Gauss-módszer növekedési tényezője exponenciálisan nő és ilyenkor mégiscsak a teljes főelemkiválasztást kell alkalmazni.

A főelemkiválasztás nélküli Gauss elimináció esetén a növekedési tényező nagyon nagy lehet. Meggyőző a 6.2. példa, ahol főelemkiválasztás nélkül , ugyanakkor pivotálással (az eredmények közti különbséget már láttuk). Egy másik – a gyakorlatból vett – példa

Itt a növekedési tényező pivotálás nélkül (MATLAB 6.5. verzióban), ugyanakkor részleges főelemkiválasztás esetén ez a szám: . (Az mátrixú egyenletrendszer megoldásában persze ez esetben nem olyan látványos az eltérés, csak , de adott esetben ekkora hiba is elfogadhatatlan.) Ezért ismét hangsúlyozzuk, hogy a Gauss-módszert főelemkiválasztás nélkül ne alkalmazzuk.

6.1.5. A Gauss-Jordan elimináció, mátrixinvertálás

Ugyanazzal a technikával, mint ahogy a -adik oszlopban az alatti elemeket kinulláztuk, a fölötte lévő elemeket is zérussá lehet tenni. Azaz az eliminációs fázisban minden értékére az ciklusváltozót nemcsak -től -ig, hanem -től -ig futtathatjuk, kivéve az esetet. (Ez annak felel meg, mintha az -nak a -adik egyenletből való kifejezése után azt nemcsak a -nál nagyobb sorszámú, hanem az összes többi egyenletbe behelyettesítenénk.) Az I. fázis végeredménye így egy diagonálmátrixú egyenletrendszer, vagyis a II. fázis ekkor csupán az () utasításokból áll (a sor -vel való osztását menet közben, egy-egy oszlop teljes kinullázása után – vagy még előtte – azonnal is megtehetjük).

Persze, szekvenciálisan végrehajtva ez a módszer nem előnyös, hiszen jelentősen megnő a műveletek száma. (Párhuzamosítva csökkenthetjük a végrehajtási időt.) Ha viszont csak azután kezdünk a főátló fölötti elemek nullázásával foglalkozni, miután kialakítottuk a felső háromszögmátrixot, és ezt a nullázást a , sorrendben végezzük (tehát az oszlopok szerint visszafelé haladva), akkor az mátrix elemeihez már nem kell hozzányúlni. Ugyanis az elvégzése során a -adik sorban az elem kivételével minden elem (elvileg) már . A -adik oszlopba sem kell a -ákat beírni (kiszámolni), mert (mint láttuk), teljesen mindegy, hogy fizikailag mi van azok helyén ábrázolva. A II. fázis úgy tekinti, hogy ott zérusok állnak. A főátló fölötti elemek nullázása tehát nem más, mint a már tárgyalt Gauss-módszer II. fázisa.

Fogalmazzuk meg most az algoritmust MATLAB programmal úgy, hogy több, de ugyanolyan együttható mátrixú (, ) egyenletrendszert oldjon meg. A jobboldali vektorokat a mátrix tartalmazza.

MATLAB program Gauss-Jordan eliminációra:

  1    function [X]=Gauss_Jordan(A,B) 
  2    [n,m]=size(B); 
  3    A=[A,B]; 
  4    % I. fazis: foatlo alatti eliminacio 
  5    % 
  6    for k=1:n-1 
  7       for i=k+1:n 
  8          [amax,ind]=max(abs(A(k:n,k))); 
  9          t=ind+k-1; 
 10          %ind az A(k:n, k) particiobeli index, t pedig az A-beli 
 11          if t  k 
 12             %a k-adik es a t-edik sor csereje: 
 13             A([k,t], :)=A([t,k], :); 
 14          end %if 
 15          g=A(i,k)/A(k,k); 
 16          A(i,k+1:n+m)=A(i,k+1:n+m)-g*A(k,k+1:n+m); 
 17       end %i 
 18    end %k; 
 19    II. fazis: foatlo folotti eliminacio 
 20    for k=n:-1:1 
 21       %A kovetkezo i ciklus k=1 eseten ures 
 22       for i=1:k-1 
 23          g=A(i,k)/A(k,k); 
 24          A(i,n+1:n+m)=A(i,n+1:n+m)-g*A(k,n+1:n+m); 
 25       end %i 
 26       X(k,1:m)=A(k,n+1:n+m)/A(k,k); 
 27    end %k; 
 28    % 
 29    %barmely j-re (j=1,2,... ,m) 
 30    %az Ax=B(:, j) egyenletrendszer megoldasa: X(:, j) 
 31    return 

A programban itt is kihasználtuk, hogy az elimináció során az mátrixban keletkező -ákat felesleges beírni (kiszámolni). Egyben szeretnénk eloszlatni azt a (több helyen is olvasható) tévhitet, hogy a Gauss-Jordan elimináció nagyobb költségű lenne (ezért hátrányosabb), mint a Gauss-módszer. Fenti program végrehajtásának műveletigénye ugyanakkora, mint a Gauss-módszeré (akkor is, ha a soroknak a főátló elemeivel való osztását az eliminációs lépés előtt hajtjuk végre; több helyen így tárgyalják a módszert). A műveletek – így a kerekítési hibák is – ugyanazok, csak a sorrendjük más (akárcsak a később tárgyalt -módszernél). A műveletigényt a 6.1.3. alfejezetben látott módon számolhatjuk össze. (Az ott megismert eredmény birtokában még egyszerűbben is, hisz csak a vektorral való műveleteket kell leszámlálni és nem -szer, hanem -szer figyelembe venni). A következő eredményt kapjuk:

esetén ez régi flop.

Fenti algoritmus alkalmas mátrixinvertálásra. Könnyen belátható ugyanis, hogy az egyenletrendszer megoldása éppen az inverz mátrix -edik oszlopvektora ( az -edik egységvektor). Ha az algoritmusban az egységmátrix, akkor a végeredmény: .

6.3. Példa. Invertáljuk Gauss-Jordan eliminációval az mátrixot, természetesen részleges főelemkiválasztást alkalmazva.

Megoldás. Keretezve jelöljük a pivotelemeket és zárójelbe zárjuk azokat, amelyek elméletileg zérusok, csak éppen nem számoltuk ki, hisz nem használjuk azokat. A kezdeti permutáció: .

A főátló alatti elemek kinullázása után következik a főátló fölötti elemek nullázása (a program 20. sora):

A bekeretezett elemekkel osztva a sorokat, a jobboldali partícióban megjelenik az inverz:

6.1.6. Az LU- és a Cholesky-módszer

A gyakorlatban sokszor könnyít a dolgunkon, ha egy mátrixot sikerül két, speciális tulajdonságú mátrix szorzatára bontani (faktorizálni).

6.12. Definíció. Az mátrix -felbontásán a mátrix szorzatra bontását értjük, ahol alsó, pedig felső háromszögmátrix.

Az -felbontás nem egyértelmű. Ugyanis a 2.3. alfejezetben leírtak alapján: ha egy nemszinguláris diagonális mátrix, akkor . Ha nemszinguláris, igazolható, hogy egyetlen -felbontásból az összes többi csak ilyen módon származtatható. Az is belátható, hogy amennyiben nemszinguláris és van egy szorzatra bontása, akkor olyan faktorizációt is találunk, melyben az vagy főátlóbeli elemeit tetszőlegesen, de -tól különböző számként előírjuk. Ilyenkor közbeszúrjuk azt a -et, amely elemeire , ha ott az főátlóbeli eleme van -nek megadva, és , ha a -edik helyen az főátlóbeli elemét írjuk elő. Ebből az is következik, hogy ha van -felbontása -nak, akkor egyértelműen van olyan is, ahol az (vagy az ) minden főátlóbeli eleme -es. Az ilyen háromszögmátrixokat egység (alsó vagy felső) háromszögmátrixoknak nevezzük.

Az egzisztenciát pedig a következő állítás fogalmazza meg.

6.13. Tétel. Egy nemszinguláris mátrixnak akkor és csak akkor létezik -felbontása, ha összes főminor mátrixa is nemszinguláris, azaz

Bizonyítás. Egyszerű számolással belátható, hogy a Gauss-elimináció során a -adik oszlop kinullázása úgy is történhet, hogy elvégezzük a szorzást, ahol

azaz csak a -adik oszlop különbözik az egységmátrixtól; abban az oszlopban pedig a -k vannak. Az első fázisban tehát a transzformáció-sorozatot hajtottuk végre, aminek a végeredménye egy felső háromszögmátrix. A mátrixok nemszinguláris alsó háromszögmátrixok. Kihasználva a mátrixszorzás asszociativitását és a szorzat inverzének tulajdonságát, írhatjuk, hogy . Tudva azt, hogy háromszögmátrixok inverze is ugyanolyan típusú (alsó- vagy felső) háromszögmátrix, továbbá azonos típusúak szorzata is az (lásd a 2.7. alfejezet 3. feladatát), már tudjuk, hogy . Egyszerű szorzással ellenőrizhető, hogy inverze a -ból úgy származtatható, hogy minden -ról levesszük a negatív előjelet. Az is könnyen igazolható (például indukcióval), hogy a szorzat úgy áll elő, hogy az egységmátrix feltöltődik a -kat tartalmazó oszlopokkal (mintha átlátszó papírra írva a mátrixokat egymásra raknánk). Tehát

Ha az -edik lépésben elakad az elimináció, akkor a (6.11) jelölését használva , de ez azt jelenti, hogy , ami a feltétel szerint nem lehet. Tehát állításunkat az „akkor“ irányban igazoltuk. A fordítottja pedig abból következik, hogy ha létezik az faktorizáció, akkor , ahol egyik tényező sem szinguláris; vagyis a baloldali minor sem az.

Az -faktorizáció tehát szoros kapcsolatban áll a Gauss-eliminációval, gyakorlati megvalósítása ugyanaz. Sőt, a szorzók és így az számára sem kell külön helyet lefoglalni, mivel azt lehet az eredeti mátrix alsó részén tárolni (az eleve az felső háromszög része lesz). Úgy történhet tehát az előállítása, hogy az eliminációs fázisban fizikailag nem írjuk be a főátló alá a -kat, hanem (mint ahogy az algoritmusunkat is fogalmaztuk) ott hagyjuk a régi értéket, és a végén oszloponként végigosztunk a főátlóbeli elemmel. Hiszen a (6.10) sorozat jelölésével a (6.11) mátrix ekkor fizikailag a következő:

és ().

Mivel az algoritmus lényegi része a Gauss-módszer eliminációs fázisa, így a MATLAB programját itt nem írjuk le mégegyszer, csupán megemlítjük, hogy ezen fázis után a tényezők kinyerése elegánsan a következő két utasítással történhet:

  U=triu(A); 
  L=tril(A)*diag(1./diag(A)); 

Vannak nemszinguláris mátrixok, amelyeknél a feltétel nem teljesül, ezért nincs -felbontásuk; például az mátrixnak sincs. Viszont korábban láttuk, hogy nemszinguláris mátrixok esetén a Gauss-módszer részleges főelemkiválasztással mindig sikeres. Mivel a sorcserék alkalmazása a sorok permutálásának felel meg, ezért kimondhatjuk a következő tételt.

6.14. Tétel. Minden nemszinguláris mátrixhoz létezik olyan permutációmátrix, hogy a mátrixnak van -felbontása.

Előbbi -es példánkban .

Az -felbontást végző MATLAB programban ilyenkor természetesen a sorcseréket is rögzítenünk kell. Helytakarékosság okán elég egy permutációs vektorban regisztrálni az eredetihez képesti pillanatnyi sorrendet (lásd a 6.4. példát). A megfelelő MATLAB rutin hívása során a visszatérési értékek között viszont maga a permutációmátrix szerepel: [L,U,P]=lu(A).

Egy alsó háromszögmátrix transzponáltja felső háromszögmátrix és ez fordítva is igaz. Felmerül a kérdés, hogy milyen mátrixoknak van olyan -felbontása, melyben a két tényező egymás transzponáltja.

6.15. Tétel. Ha az mátrix szimmetrikus és pozitív definit, akkor létezik szorzatra bontása, ahol alsó háromszögmátrix. Ezt a felbontást nevezzük Cholesky-felbontásnak.

6.16. Megjegyzés. Néhol a Cholesky-felbontást alakban írják és vannak szoftverek (például a MATLAB), amelyek az tényezőt állítják elő. A szoftverek használata előtt tehát ajánlatos tájékozódni.

A Cholesky-felbontást úgy is megkaphatjuk ezen mátrixoknál, hogy elkészítjük Gauss-eliminációval az -felbontást, majd megkeressük azt a diagonál mátrixot, amelynél az szorzatban az és az diagonális elemei egyenlők. Az ilyen eljárásnak azonban nincs értelme, csak a műveletszámot növeli jelentősen, hisz a Cholesky-felbontásban elég az egyik tényezőt kiszámítani, ami durván fele annyi költséggel jár. Az

egyenletből , és , azaz és .

Pszeudókóddal felírva a Cholesky-felbontás algoritmusát:

Észrevehetjük, hogy akármelyik eleme kiszámításához felhasznált -beli elemre az algoritmus később már nem hivatkozik, ezért helytakarékosan szokás implementálni: az indexes elemek helyett írható a megfelelő indexű elem, így az az alsó háromszögrészén alakul ki, míg a főátló fölött megmaradnak eredeti elemei.

A Cholesky-felbontás költsége a fenti implementálásban

A MATLAB program megírását gyakorló feladatként az olvasóra bízzuk. A beépített MATLAB rutin hívása: U=chol(A).

A megoldás LU- és Cholesky-módszerrel.

Tekintsük az megoldását, ahol nemszinguláris. A 6.14. tétel szerint van permutációmátrix, hogy létezik a faktorizáció és írható, hogy . Bevezetve az új változót, végrehajtható tehát a következő:

Az -módszer algoritmus

1. Határozzuk meg a felbontást.

2. Oldjuk meg az egyenletrendszert -ra.

3. Oldjuk meg az egyenletrendszert -re.

A 2. és a 3. lépésben háromszögmátrixú egyenletrendszert kell megoldani, tehát az összköltséget dominánsan az 1. lépés határozza meg. Szimmetrikus, pozitív definit esetén az 1. lépésben természetesen a Cholesky-felbontást alkalmazzuk.

6.17. Megjegyzés. A mátrixot nem kell előre ismernünk. Az 1. lépésben a faktorizációt a főelemkiválasztást alkalmazó Gauss-eliminációval végezzük el, és ott regisztráljuk a sorcseréket (és az oszlopcseréket is, teljes főelemkiválasztás esetén; ekkor egy faktorizációt végzünk és az megoldásnál vesszük figyelembe a -et). Így a fentebb leírt -módszer lényegében (beleértve a kerekítési hibákat is) ugyanaz, mint a főelemkiválasztásos Gauss-módszer. A különbség ott jelentkezik, hogy az -val nem párhuzamosan transzformáljuk a jobboldali vektort, hanem azt elhalasztjuk. A 2. lépés pontosan a transzformációját hajtja végre. A mátrix helyett hely- (és adatmozgatás) takarékossági okokból elég csak egy permutációvektort tárolni (lásd a 6.4. példát).

MATLAB program az LU-módszerre:

  1    function [x]=linermegoldluval(A,b) 
  2    n=size(b); 
  3    [L,U,P]=lu(A); 
  4    Itt a beepitett rutint hasznaltuk a faktorizaciora 
  5    b=P*b 
  6    y(n:-1:1)=fhmatrixszal(L(n:-1:1, n:-1:1), b(n:-1,1)); 
  7    x=fhmatrixszal(U, y); 
  8    return 

A 6.-ik sorban kihasználtuk, hogy amennyiben egy alsó háromszögmátrixot tükrözünk a vízszintes tengelyre, majd a függőleges tengelyre is (más szóval felcseréljük a sorainak és az oszlopainak sorrendjét is), akkor egy felső háromszögmátrixot kapunk, az egyenletrendszer így olyan felső háromszögmátrixú lesz, ahol az egyenletrendszer skalár egyenleteinek sorrendje fordítottja az eredetinek, és a megoldás komponenseit is fordított sorrendben kapjuk vissza. A tükrözéseket ábrán is megmutatva:

6.3. ábra - Alsó háromszögmátrixból kettős tükrözéssel felső

Alsó háromszögmátrixból kettős tükrözéssel felső

6.4. Példa. Oldjuk meg ismét a 6.1. példában már látott egyenletrendszert, de most -módszerrel, részleges főelemkiválasztást alkalmazva.

Megoldás. Az együtthatómátrix megegyezik a 6.3. példában invertált mátrixszal. A főátló alatti elemek nullázásához vezető lépések ezért megegyeznek az ott már közöltekkel. Mindazonáltal hasznosnak tartjuk, ha egy-egy feladat megoldását mindig teljes egészében látjuk, ezért azt a részt közöljük itt is. A kezdeti permutáció: .

Az -módszert különösen akkor előnyös használni, ha egynél több, előre nem ismert jobboldalú

alakú egyenletrendszert kell megoldani. Ekkor elég az mátrix -felbontását egyszer meghatározni, majd rendre az , (), összesen darab háromszögmátrixú egyenletrendszert megoldani. Alkalmazhatjuk például mátrix invertálására, az előre rögzített egységvektorokkal (). Ekkor az eljárás egyébként teljesen ekvivalens a Gauss-Jordan elimináció 6.1.5. alfejezetben ismertetett módon való alkalmazásával.

6.2. Iteratív módszerek

6.2.1. A Jacobi-iteráció

Jelöljük most az egyenletrendszer pontos megoldását -gal. Ebben az alfejezetben egy olyan módszer két változatát ismertetjük, amelyek az -hoz konvergáló sorozatokat állítanak elő, ahol

Szükségünk lesz a komponensek feltüntetésére, ezért a sorozatbeli sorszámot felső indexszel jelöljük. Konvergencián a komponensenkénti konvergenciát értjük, ami véges ismeretlenszám esetén egyenértékű azzal, hogy valamilyen normában fennáll az .

Az ismertetni kívánt eljárás az

alakú egyenletrendszerek megoldására alkalmas, ezért először ilyen alakra kell hozni. Egy szokásos átalakítást később bemutatunk.

Maga az eljárás egyszerű. Vegyünk fel egy kezdeti vektort. A sorozat többi elemét az alábbi rekurzióval állítjuk elő:

A következő ábrasorozat az első három iterált bemutatásával geometriailag szemlélteti a viszonyokat három dimenzióban konvergencia esetén (a zsúfoltság elkerülése érdekében a pontos megoldás helyét csak a kék színű karakterrel jeleztük):

6.4. ábra - Iteráció 1. lépése

Iteráció 1. lépése

6.5. ábra - Iteráció 2. lépése

Iteráció 2. lépése

6.6. ábra - Iteráció 3. lépése

Iteráció 3. lépése

Kérdés, hogy milyen feltételek mellett konvergál a sorozat az megoldáshoz és a -adik elemének mekkora a hibája.

6.18. Tétel. Ha , akkor a sorozat tetszőleges választása esetén az egyetlen megoldáshoz konvergál. Érvényes továbbá a következő hibabecslés:

6.19. Megjegyzés. A feltétel azt jelenti, hogy az leképezés a teljes téren kontrakció, abban az értelemben, hogy bármely két különböző pont képe közelebb van egymáshoz, mint az eredeti két pont. Tehát, ha , akkor

Az egyenletrendszer megoldása pedig az a pont, aminek a képe saját maga, azaz a leképezés fixpontja. Éppen ezért ezt a módszert szokták fixpontiterációnak is nevezni.

Ha van elképzelésünk az megoldásról, akkor azt érdemes figyelembe venni az megválasztásakor, egyébként az választás a szokásos. (Lehetne is, de ekkor előre tudjuk a következő iteráltat, ami éppen a ).

A sorozat elemeinek kiszámítása leggyakrabban komponensenként történik, vagyis

A fenti (6.18) formula egyébként sugallja is az egyenletrendszernek iterációra alkalmas egyfajta átalakítását, nevezetesen: fejezzük ki az -edik skaláregyenletből az -t. Ekkor

Könnyen látható, hogy ilyen átalakítással a 6.18. tétel feltétele akkor és csak akkor teljesül, ha diagonálisan domináns (lásd a (6.9) képletet). Elképzelhető, hogy az eredeti skalár egyenleteinek felcserélése után lesz diagonálisan domináns az együtthatómátrix, de sok egyenletnél az iteráció nem alkalmazható. Mindenesetre vannak a gyakorlatban is előforduló (például differenciálegyenletek numerikus megoldása során keletkezett) egyenletrendszerek, amelyeknél lehet és célszerű is az iterációval való közelítés.

MATLAB program a Jacobi-iterációra:

  1    function [x, itszam]=Jacobi(B, c, x0, q, tol, itmax) 
  2    p=q/(1-q); 
  3    % 
  4    d=inf; % hogy legalabb egy iteraciot vegrehajtson 
  5    itszam=0; 
  6    while (d>tol) & (itszam<itmax) 
  7       itszam=itszam+1; 
  8       x=B*x0+c; 
  9       d=p*norm(x-x0, inf); 
 10       x0=x; 
 11    end %while 
 12    return 

Iteratív ciklusoknál mindig fel kell készülni arra, hogy az elméleti kilépési feltétel numerikusan nem érhető el, ezért a lépésszámra „illik“ egy felső korlátot előírni (itt most ez az itmax változó).

6.5. Példa. Oldjuk meg iterációval (ha az konvergens) az alábbi egyenletrendszert pontossággal:

Megoldás. A pontos megoldás: , majd ezzel is összehasonlítjuk a közelítéseket. Látható, hogy az egyenleteket a sorrendbe írva, az együttható mátrix diagonálisan domináns lesz, ezért mindjárt ebben a sorrendben alkalmazva a (6.19) összefüggést, a (6.18) egyenletrendszer a következőképpen adódik:

, tehát valóban konvergens az iteráció. Induljunk ki az vektorból. A (6.18) előírás szerint, egymás után a következő közelítések adódnak:

A közelítések után zárójelben feltüntettük a (6.16) szigorúbb egyenlőtlensége szerint becsült hibát, mellette a tényleges hibát. Meglehetősen lassú a konvergencia, aminek a viszonylag nagy értéke az oka. Azt is észrevehetjük, hogy itt a (6.16) elég durva becslést ad, egy nagyságrenddel durvábbat, mint a tényleges hiba.

Az iteráció gyorsítására több lehetőség ismert. Ezek közül mutatunk be egyet.

6.2.2. A Seidel-iteráció

Az eljárás alapgondolata (a numerikus módszerek más területein is alkalmazott, ún. Seidel-féle relaxációs elv), hogy mindig a legutolsó információt használjuk fel. A (6.18) egyenletrendszer általános tagját írjuk fel a következőképpen:

Ha egymás után számoljuk az egyes komponensek újabb közelítését, akkor az első szummában a komponenseknek már a -adik iteráltja is megvan, oda azokat lehet behelyettesíteni. Tehát:

MATLAB program a Seidel-iterációra:

  1    function [x, itszam]=Seidel(B, c, x0, q, tol, itmax) 
  2    n=length(c); 
  3    p=q/(1-q); 
  4    %hogy legalabb egy iteraciot biztos vegrehajtson: 
  5    d=inf; 
  6    itszam=0; 
  7    x=x0; 
  8    while (d>tol) & (itszam<itmax) 
  9       itszam=itszam+1; 
 10       for i=1:n 
 11          x(i)=B(i,:)*x; 
 12       end % i
 13       d=p*norm(x-x0, inf); 
 14       x0=x; 
 15    end % while 
 16    return 

A Seidel-iterációra ugyanaz a konvergenciatétel érvényes, mint a korábban tárgyalt Jacobi-iterációra.

6.6. Példa. Oldjuk meg most Seidel-iterációval az előbbi egyenletrendszert, pontossággal:

Megoldás. Ugyanazon átalakítás után, ugyancsak az választással a következőt kapjuk:

Kevesebb lépésben értük el az előírt hibakorlátot, mint a korábban tárgyalt változattal, de nagyságrendi különbséget nem érzékelünk.

6.20. Megjegyzés. Sorozatok konvergenciasebességét nem a lépésszámmal szokták jellemezni. Itt most mellőzzük a konvergenciasebesség fogalmának részletezését, csupán érzékeltetni kívánjuk. Bizonyos feltételekhez kötött kategóriákat állítunk fel; hétköznapi hasonlattal: vannak gyalogos-, személyautó-, repülőgép-sebességűek, stb. Többféle konvergenciarend fogalmat ismerünk. A leggyakrabban használt konvergenciarend azzal jellemezhető, hogy hibakorlátja az hibájának hanyadik hatványával arányos. A (6.16) becslésből következik, hogy , tehát mindkét változatban a sorozat konvergenciasebessége elsőrendű vagy másképpen lineáris. De mint ahogy a személyautók között is van különbség, úgy az azonos rendben konvergens módszerek között is. A Seidel eljárás abban az értelemben gyorsabb, hogy azonos kezdővektorból ugyanazt a hibakorlátot rendszerint kevesebb lépésben éri el.

6.21. Megjegyzés. A 6.18. tétel a konvergenciára elégséges feltételt mond ki. Ha nem teljesül, az még nem zárja ki a konvergenciát, de a (6.16) hibabecslés ekkor konvergencia esetén sem alkalmazható. Szükséges és elégséges feltételt is ki lehet mondani, de az abban szereplő ún. spektrálsugár megállapítása nehéz feladat.

Végül pár szó arról, hogy miért van szükség közelítő módszerre, ha van a megoldást elvileg pontosan előállító eljárás is, például a Gauss-elimináció. Először is azért, mert az valóban csak elvileg ad pontos megoldást, a gyakorlatban a kerekítési hibák halmozódása nagyon nehezen megítélhető módon eltérítheti a pontos megoldástól. Szemben az iterációval, ahol a kerekítési hibák nem halmozódnak, és amennyiben a konvergenciatétel feltételei teljesülnek, megbízható és könnyen elvégezhető hibabecslést ad a kezünkbe. Azonkívül: a gyakorlatban mindig elegendő valamilyen közelítés, azt pedig adott esetben iterációval olcsóbban is elérhetjük, hiszen egy lépés költsége csupán flop.

6.3. Feladatok

1. Oldjuk meg Gauss-eliminációval, részleges és teljes főelemkiválasztással is az alábbi egyenletrendszereket:

(i)     ,

(ii)     .

2. Oldjuk meg -módszerrel, részleges főelemkiválasztással az alábbi egyenletrendszert:

Számítsuk ki az együttható mátrix inverzét is Gauss-Jordan eliminációval.

3. Oldjuk meg Jacobi- és Seidel-iterációval is az alábbi egyenletrendszert, végtelen normában hibakorláttal:

Hasonlítsuk össze az iteráció két változatával kapott eredményeket és az iteráció számokat.

Oldjuk meg Gauss-eliminációval is. Mekkora a tényleges hibája az iterációval kapott közelítéseknek?

4. Tekintsük az egyenletrendszert, ahol

Oldjuk meg kétféleképpen is az egyenletrendszert:

(i) Gauss-eliminációval, részleges főelemkiválasztást alkalmazva,

(ii) végigszorozva az egyenletet -vel és az egyenletnél Cholesky-módszert alkalmazva.

Figyelembe véve az -vel való szorzást is, melyik módszer a költségesebb és mennyivel?

5. Adott az

mátrix.

(i) Adjuk meg azt az -felbontását, ahol , , továbbá teljesül.

(ii) Adjuk meg azt az szorzatra bontást, amelyben egység alsó, egység felső háromszögmátrix, pedig diagonál mátrix.

7. fejezet - Egyenletrendszerek hibaanalízise

A vizsgálatban a direkt és az inverz hibákat elemezzük. Az egyenletrendszer megoldásával kapcsolatban a következő jelöléseket és fogalmakat használjuk. Az elméleti megoldást , a közelítő megoldásokat pedig jelöli. A közelítő megoldás direkt hibája . Az mennyiséget reziduális hibának nevezzük. Az elméleti megoldás esetén , a közelítő megoldás esetén pedig

Az inverz hiba meghatározásához különféle modelleket használunk. A legáltalánosabb esetben feltesszük, hogy az számított megoldás kielégíti az egyenletrendszert, ahol és . A és mennyiségeket inverz hibáknak nevezzük.

7.1. Érzékenységvizsgálat

Mint láttuk az 5. fejezetben, érzékenységen azt értjük, hogy milyen mértékben reagál az output adat az input adatok (a feladat bemenő paraméterei) kismértékű változására. Lineáris egyenletrendszer paraméterei az és a , ezért három modellt állíthatunk fel, annak megfelelően, hogy csak az egyik, csak a másik, vagy mind a két paraméter változik. A lineáris egyenletrendszerek érzékenysége szoros összefüggésben van az együtthatómátrix kondíciószámával, amit ezen mennyiség bevezetésének motivációja is érzékeltet:

Legyen nemszinguláris mátrix, és tekintsük az függvényt, ahol (azaz az lineáris egyenletrendszer megoldása a jobboldali vektor függvényében). Egyszerű számolással megmutatható, hogy ezen függvény kondíciószáma , és ez a becslés pontos. Joggal nevezhetjük tehát a kifejezést az mátrix kondíciószámának.

7.1. Megjegyzés. A legtöbb, változóiban lineáris feladat kapcsolatba hozható egy mátrixszal, amelynek kondíciószáma egyben a feladat kondíciószámára is rámutat. Kézenfekvő példa erre a mátrix kondíciószámának előbb megmutatott motivációja, de megemlíthetnénk például az 5.5. példát is. Ennél az (5.6) képletben szereplő és együtthatókat az

egyenletrendszer megoldása adja, ahol az együtthatómátrix kondíciószáma (például) esetén .

7.1.1. Inverz hiba a jobboldali vektorban

Tegyük fel, hogy az egyenlet helyett a perturbált

egyenletrendszert oldjuk meg. Legyen és vizsgáljuk a két megoldás eltérését. A relatív direkt hibáról szóló következő tétel bizonyítását is bemutatjuk, hogy minél erősebben demonstráljuk a mátrix kondíciószáma jelentőségét. Egyben arra is fel szeretnénk hívni a gyakorlati szakemberek figyelmét, hogy más feladatok esetén sem biztos, hogy a megoldás kis környezetében járunk, ha valamilyen egyenlet két oldala majdnem egyenlő. (Persze, sokszor nem tehetünk mást, csak bízhatunk benne, de mindig érdemes összevetni ilyenkor is az eredményt más úton – például méréssel – nyert tapasztalatokkal.)

7.2. Tétel. Ha nemszinguláris és , akkor

ahol az mátrix kondíciószáma.

Bizonyítás. Az egyenlőségből , ahonnan következik. Másrészt , ahonnan . A két egyenlőtlenséget összeszorozva kapjuk, hogy

ami éppen a bizonyítandó állítás.

A tételből látható, hogy az mátrix kondíciószáma erősen befolyásolhatja az perturbált megoldás relatív hibáját. Az egyenletrendszert jól kondicionáltnak nevezzük, ha kicsi és rosszul kondicionáltnak nevezünk, ha nagy. Értelemszerűen a nagy és kicsi jelzők relatívak és környezetfüggők. A kondíciószám függ a normától. Ennek megfelelően például . Ha a normától való függés lényeges, akkor ezt külön jelöljük.

A most vizsgált modellben az inverz hiba , a 7.2. tétel pedig a relatív direkt hibára ad becslést. Ez teljes összhangban van az (5.4) hibaszámítási ökölszabállyal. A tétel állítása az összefüggés miatt átírható az ekvivalens

alakba. Az egyenlőtlenség jelentése az, hogy az perturbált megoldás relatív hibája kicsi, ha kondíciószáma kicsi és az relatív reziduális hiba kicsi. Ha azonban a rendszer rosszul kondicionált, akkor ez nem szükségképpen igaz.

7.1. Példa. Tekintsük az alakú

egyenletrendszert, ahol kicsiny és legyen az . Határozzuk meg reziduális és relatív hibáját végtelen normában. Adjunk magyarázatot a tapasztaltakra.

Megoldás. Az egyenletrendszer pontos megoldása bármilyen esetén ( esetén ez a megoldás egyértelmű is). A reziduális hiba , és , ami -re vetítve a relatív hibája is (). Ugyanakkor , azaz .

A magyarázat természetesen az, hogy kondíciószáma rendkívül nagy. Például esetén . Ha az egyenleteket, mint egyenes egyenleteit ábrázoljuk az kordinátarendszerben, akkor a két egyenes közel párhuzamos. Kis változás az együtthatókban ugyan az egyeneseket is kicsit mozdítja, de a metszéspont nagyon távol kerülhet a perturbáció előttitől. Ez egyben a kondicionáltság egy lehetséges geometriai értelmezését is adja.

7.1.2. Inverz hiba az együtthatómátrixban

Tegyük most fel, hogy az egyenletrendszer helyett a perturbált

egyenletrendszert oldjuk meg. Legyen ismét és . A kérdésre, hogy adott esetén van-e egyáltalán inverz hiba úgy, hogy teljesüljön, igen a válasz Sőt, sokféle lehet, ilyenkor – mint korábban is említettük – a valamilyen normában legkisebb az érdekes. De kicsi kondíciószám mellett minden, szintén elegendően kicsiny inverz hibára igaz a relatív reziduális hibára vonatkozó következő állítás.

7.3. Tétel. Ha , , nemszinguláris, és , akkor

Ha rosszul kondicionált, akkor az állítás nyilvánvalóan nem igaz.

7.1.3. Inverz hiba az együtthatómátrixban és a jobboldalon is

Itt most az egyenlet mindkét oldalán lévő paraméterek egyidejű perturbációjának következményét vizsgáljuk. Legyen

7.4. Tétel. Ha nemszinguláris, és , akkor

A bizonyítást, akár az előző tételnél, itt is mellőzzük.

7.5. Megjegyzés. A tételből kiolvasható a következő „ökölszabály”: Tegyük fel, hogy . Ha és elemei decimális jegyre pontosak és , ahol , akkor a számított megoldás kb. jegyre lesz pontos.

7.2. Példa. Tekintsük a következő egyenletrendszert, amelynek együtthatóit a MATLAB pontosan ábrázolja:

A feladat pontos megoldása , . A MATLAB paranccsal kapott közelítő megoldás , , amelynek relatív hibája

Minthogy és az eredmény lényegében megfelel az ökölszabálynak.

Az említett tételek gyakorlati alkalmazhatóságát igen megnehezíti a és a , de leginkább a kondíciószám ismeretének hiánya. A kondíciószám olcsó becslésére vannak ajánlott, tapasztalati módszerek. Ezekkel könyvünkben nem foglalkozunk, csak utalunk a szabadon hozzáférhető LINPACK programcsomagra. (Számos egyéb, gyakori feladatokra írt, kipróbált programok is letölthetők.) Az inverz hibákra több helyen ajánlják a , választást, ahol az abszolútérték elemenként értendő, pedig a gépi epszilon. Ezek a közelítések azonban nem mindig jogosak, ezért is van szükség elméletileg megalapozottabb közelítésekre, illetve felső korlátokra. A mindig használható, becslésére pedig a következő alfejezetben térünk ki.

7.2. A relatív hiba becslése Wilkinson tételével

A becslésére több eredmény is született, itt most Wilkinson eredményét mutatjuk be, részben a pivotelemek növekedési tényezőjének demonstrálása céljából is.

7.6. Tétel. Az egyenletrendszer Gauss-módszerrel lebegőpontos aritmetikában kapott közelítő megoldása kielégíti az

egyenletrendszert, ahol

és a pivot elemek növekedési tényezője, pedig az aritmetika egységnyi kerekítése.

Pl. a (7.2) példában figyelembe véve, hogy , , pedig -re adódott, a eredményt kapjuk.

Minthogy a gyakorlatban kicsi (főelemkiválasztás esetén), a

relatív inverz hiba is az.

A tételből kapjuk, hogy

Kis kondíciószám esetén feltehetjük, hogy , így a 7.4. tételből esetén a relatív direkt hibára az alábbi közelítő becslést kapjuk:

Ez az összefüggés is rámutat, hogy a feladatérzékenység és a numerikus stabilitás egymástól független fogalmak: külön-külön befolyásolják a számított eredmény pontosságát. Lineáris egyenletrendszernél a feladatérzékenységet a , míg a Gauss-módszer numerikus stabilitását a jellemzi. A 6.2. példában sem a feladatérzékenység okozta a főelemkiválasztás nélküli esetben keletkezett durva hibát, hiszen ott , hanem a pivotelemek növekedési tényezője (korábban már megmutattuk, hogy az említett példában ).

A (7.12) a kondicionáltságnak a számított eredmény hibájára gyakorolt hatása mellett Gauss-módszer esetén a pivotelemek növekedése fontosságát igazolja elméletileg is.

7.3. Utólagos hibabecslés

A valamilyen módszerrel kapott közelítő megoldás hibájának utólagos becslésére azért van szükség, hogy elképzelésünk lehessen az eredmény megbizhatóságáról. Itt csupán a gyakorlatban jól használható Auchmuty-tételt ismertetjük, bizonyítás nélkül.

7.7. Tétel. Jelölje az egyenletrendszer valamilyen módon kiszámított közelítő megoldását. Ekkor igaz, hogy

ahol a (konstansnak nevezett) az -tól és az hibavektor irányától függ (a nagyságától nem), továbbá .

Az összefüggés alapján a reziduális hiba és az hányados meghatározásával nagyságrendileg helyesen becsülhetjük a közelítő megoldás abszolút hibáját. A számítógépes tapasztalatok azt mutatják, hogy nem túl nagy szám, gyakran . De ha kételyeink vannak és van időnk, költségünk az adott -hoz -nek megbízhatóbb becslésére, úgy magunk is kísérletezhetünk. Az adott mátrixnál több vektort felvéve meghatározzuk a jobboldali vektorokat, majd véletlenszerű perturbációkat választva kiszámoljuk az reziduális hibákat, és a (7.13)-ból kifejezve kapjuk -ket. Ezek közül nyilván a legnagyobbat választjuk, egy biztonsági tényezővel (pl. -zel) szorozva és ezzel felső becslését a tétel segítségével megbízhatóan kaphatjuk.

7.3. Példa. Legyen

(Ez a Lotkin-mátrix.) Írjunk MATLAB programot az -hoz tartozó meghatározásához, és futtassuk -re, kísérletet végezve.

Megoldás. Az mátrix előállítását az olvasóra bízzuk.

MATLAB program az Auchmuty-tételre:

  1    function c=auchm(A,k) 
  2    n=length(A); 
  3    c=zeros(n,1); 
  4    rand('state',0); 
  5    for i=1:k 
  6       x=rand(n,1); 
  7       b=A*x; 
  8       dx=2e-10*(2*rand(n,1)-1); 
  9       xp=x+dx; 
 10       r=A*xp-b; 
 11       c(i)=norm(dx)*norm(A'*r)/norm(r)^2; 
 12    end 
 13    return 

A futtatás végeredménye:

.

Itt a legnagyobb érték a , az említett biztonsági tényezővel, kerekítve a megbízhatóan választható a felső korlát becsléséhez.

7.8. Megjegyzés. Sok kísérletet persze nem érdemes végezni, hisz költséges lehet, márpedig általánosan is kimondhatjuk: hibabecsléshez csak az eredeti probléma megoldása költségétől nagyságrendben alacsonyabb költségű módszert illik alkalmazni.

7.4. Iteratív javítás

Az egyenletrendszer valamilyen módszerrel (például Gauss-módszerrel) kapott javítására is születtek ötletek, ezek egyikét, mint iteratív módszert ismertetjük. Alapgondolata a következő. Legyen a pontos megoldás és a közelítő , a köztük fennálló összefüggéssel. Ekkor az egyenletből a meghatározása után ismert lehetne a pontos megoldás. Csakhogy nyilván -nek is csak egy közelítését határozhatjuk meg, így azzal nem pontossá tehetjük, hanem reményeink szerint is csupán javíthatjuk az közelítést.

Az iteratív javítás algoritmusa

Legyen

for

    

    Számítsuk ki az egyenletrendszer közelítő megoldását

    

end

Az eljárás különféle változatai ismertek. Az egyenletrendszerek megoldására az -módszer előnyös, hisz több, azonos együtthatójú olyan rendszert kell megoldani, melyek jobboldalai előre nem ismertek. Kilépési feltételként a egy korlát alá kerülését ajánlják, ám a tapasztalatok szerint sem a , sem az nem csökken monoton, ezért aztán elég szigorú korlátot egész váratlanul vagy egyáltalán nem ér el az algoritmus. Néhány lépést érdemes elvégezni és ezek közül kiválasztani a legkisebb reziduális hibájút. Mindenesetre elméletileg igazolt, hogy az iteratív javítás még egyszeres pontosságú lebegőpontos aritmetikában is javíthatja a közelítő megoldás pontossága.

7.4. Példa. Írjunk MATLAB programot iteratív javításra lépésben, majd javítással oldjuk meg az egyenletrendszert, ahol

az ú.n. Faddeev-Sominskii mátrix, és

(Az mátrix előállítását itt is az olvasóra bízzuk.)

Megoldás.

MATLAB program az iteratív javításra:

  1    function [rnorm,xlegjobb]=itjav(A,b,k) 
  2    n=length(b); 
  3    rnorm=zeros(k+1,1); 
  4    %egy kozelito megoldas + k javitas 
  5    [L,U,P]=lu(A); 
  6    b1=P*b; 
  7    y(n:-1:1)=fhmatrixszal(L(n:-1:1, n:-1:1), b1(n:-1:1)); 
  8    x=fhmatrixszal(U, y); 
  9    xlegjobb=x; 
 10    r=A*x-b; 
 11    rlegjobb=norm(r,inf); 
 12    rnorm(1)=rlegjobb; 
 13    % 
 14       for i=1:k 
 15  15         h=P*r; 
 16  16         y(n:-1:1)=fhmatrixszal(L(n:-1:1, n:-1:1),h(n:-1:1)); 
 17  17         d=fhmatrixszal(U, y); 
 18  18         x=x-d; 
 19  19         r=A*x-b; 
 20  20         rnorm(i+1)=norm(r,inf); 
 21  21         if rlegjobb>rnorm(i+1) 
 22  22            rlegjobb=rnorm(i+1); 
 23  23            xlegjobb=x; 
 24  24         end % if 
 25  25      end % i 
 26  26   return 

(A programban többször is használt fhmatrixszal nevű függvénnyel kapcsolatban utalunk a 6.1.1. alfejezetben megírt programra és annak a 6.1.6. alfejezetbeli alkalmazására.)

Hívás: [rnorm,xlegjobb]=itjav(A,b,4). A kapott eredmények:

Az eredményekből látható, hogy itt az első lépés javítás helyett éppen rontott, a második javító lépés szolgáltatta a legkisebb reziduális hibát ( normában), és így azt az eredményt őrizte meg a program, mint a megoldás legjobb közelítését, majd ismét elkezdett a sorozat romlani.

7.5. Feladatok

1. Oldjuk meg az

egyenletrendszert és vizsgáljuk a megoldás hibáját.

2. Tekintsük az

lineáris egyenletrendszert. Ennek pontos megoldása 8 tizedesjegyre helyesen kerekítve:

(i) Oldjuk meg az egyenletrendszert 8 tizedesjegy hosszú mantisszájú lebegőpontos aritmetikában és magyarázzuk meg, hogy a kapott eredmény pontossága miért annyira rossz. (A részleges főelemkiválasztás itt semmitmondó, teljes főelemkiválasztást pedig ne alkalmazzunk.)

(ii) Oldjuk meg az egyenértékű

egyenletrendszert részleges főelemkiválasztással. Mit tapasztalunk?

(iii) Számoljuk ki mindkét feladat mátrixának kondíciószámát végtelen-normában.

3. Tekintsük az (i) és az (ii) egyenletrendszert, ahol

pedig kétkomponensű vektor.

A két egyenlet közül melyik érzékenyebb a perturbációjára? Az érzékenyebb egyenletnél a -t hányszor kisebb relatív hibával kell ismernünk, hogy ugyanolyan pontossággal határozhassuk meg a megoldást, mint a másiknál?

4. Legyen . Oldjuk meg az egyenletrendszereket a sima és a részleges főelemkiválasztásos Gauss módszerrel is az alábbi adatok esetén. Határozzuk meg a kondíciószámokat és vizsgáljuk meg a kapott megoldások hibáját a fejezet módszereivel iteratív javítással és anélkül.

(i)

(ii)

(iii)

Mi a kapcsolata az egyes egyenletrendszereknek?

5. Legyen , és , ().

(i) Adjuk meg kondíciószámát és a relatív inverz hibát függvényében.

(ii) Szintén az függvényében mekkora az relatív reziduális hibája (azaz az érték, ahol ), ha az

pontos megoldása?

6. Legyen az ismeretlenes lineáris egyenletrendszer pontos megoldása, pedig egy tetszőleges komponensű vektor. Legyen továbbá az reziduális hibája: .

Igazoljuk, hogy indukált normában fennállnak a következő egyenlőtlenségek:

Értelmezze az egyenlőtlenségeket.

7. Tekintsük az lineáris egyenletrendszert, ahol a negyedrendű Hilbert mátrix, vagyis . Legyen továbbá . Ismert, hogy rosszul kondicionált, ezért az inverzét approximáljuk -vel, ahol

Tehát az megoldás egy közelítése: . Ez nyilván nem a pontos eredmény, de még csak nem is elfogadható közelítés, mert tudjuk, hogy a pontos megoldásnak is csak egész komponensei vannak. Alkalmazzuk az iteratív javítást az elfogadható integer megoldás megkeresésére, az helyett annak közelítését használva.

8. Tekintsük az lineáris egyenletrendszert, ahol

Mivel szimmetrikus, pozitív definit, ezért a Cholesky módszerrel oldjuk meg. Adjuk meg a pontos felbontást és az egyenletrendszer pontos megoldását.

9. Oldjuk meg az egyenletrendszert, ahol

és

Számítsuk ki a mátrix determinánsát, kondíciószámát és becsüljük is meg a közelítés hibáját a reziduális hibabecslés alapján. A közelítő megoldás elfogadható relatív hibája: 5%.

8. fejezet - A MATLAB-ról, röviden

A MATLAB matematikai szoftver a nevét a MATrix LABoratory kifejezésből kapta. A neve is arra utal, hogy benne a mátrixszámítások, illetve általában a lineáris algebrai algoritmusok rendkívül egyszerűen programozhatók. Természetesen a MATLAB szoftver állandó fejlődésen megy keresztül, mai eszköztára lényegesen bővebb, felsorolásukra itt nem vállalkozhatunk, csupán utalunk az irodalomjegyzékre. Ezt a fejezetet azért építettük a könyvünkbe, hogy valami képet mégis kapjunk a konkrét lehetőségekről, a korábbi fejezetek megértését segítsük, valamint, hogy elindítsuk az olvasót a MATLAB-nak – akár a saját Help-jén keresztüli – tanulmányozására.

A MATLAB Help-jében a leírások témák szerinti csoportosításban, különböző alkönyvtárakban vannak tárolva text, illetve hipertext (HTML fájl) formában. A parancsablakban segítséget kérhetünk utasítások közvetlen begépelésével is, vagy a menüsorban a ?, illetve a Help pontokra kattintással. Három szintje van a segítség kérésnek, itt csak a legegyszerűbbet ismertetjük:

A help kiadása után a leírásokat tartalmazó alkönyvtárak listája jelenik meg. Egy részlet például:

                

matlab\elfun - Elementary math functions.

matlab\specfun - Specialized math functions.

matlab\matfun - Matrixfunctions-numerical linear algebra .

                

A help alkönyvtár parancs hatására megjelenik az ott található MATLAB kulcsszavak listája, témák szerint csoportosítva. Például a help matfun-ra adott válasz egy részlete:

                

chol - Cholesky factorization.

cholinc - Incomplete Cholesky factorization.

lu - LU factorization.

                

Végül a help kulcsszó hatására kapjuk meg a konkrét kulcsszó MATLAB jelentésének leírását. Például a help lu-ra adott válasz egy részlete

LU        LU factorization.

                

[L,U,P] = LU(X) returns unit lower triangular matrix L, upper

triangular matrix U, and permutation matrix P so that

P*X = L*U.

                

See also COLAMD, LUINC, QR, RREF, UMFPACK.

Persze, ha csak a megoldandó feladatot ismerjük, de ahhoz kapcsolódó kulcsszavakat nem, akkor más utat kell választani. Ilyenkor nagyon hasznos a lookfor (keresd!) parancs. Általunk megadott szövegrészletet keres a kulcsszavak leírásának első sorában. Például, ha az integrálási lehetőségekről akarunk tájékozódni, a lookfor integr parancs után bőséges listát kapunk:

                

POLYINT Integrate polynomial analytically.

DBLQUAD Numerically evaluate double integral.

                

8.1. Mátrixok és vektorok a MATLAB nyelvben

A MATLAB legfontosabb adattípusa a valós vagy komplex elemű mátrix. Mindazonáltal a skalárokra és vektorokra nem kell (bár lehet) két indexszel hivatkozni, továbbá az újabb verziókban a többindexes tömbök és egyéb objektumok is megjelentek.

A mátrixokat legegyszerűbben elemeik sorfolytonos felsorolásával adhatjuk meg. A sorokat pontosvesszővel vagy enter-rel választjuk el, ugyanazon mátrixsorban álló elemeket pedig helyközzel vagy vesszővel. Az egész felsorolást szögletes zárójelbe tesszük. Pl. az

mátrixot megadhatjuk a következőképpen:

              

Mindjárt az értékadásra is láttunk példát. A mátrixok deklarációja a rájuk, vagy egy elemükre vonatkozó első értékadó utasítással automatikusan megtörténik. Azonosítót a szokásos módon választhatunk. Vannak foglalt azonosítók, ezek a HELP parancs segítségével megismerhetők. Most csak az eps foglalt azonosítót említjük ami a konkrét számábrázolástól függő, ún. gépi epszilon (rendszerint ) értékét jelenti. Nehezen kideríthető hibát okozhat a jelentésétől eltérő használata.

Egy mátrixot kisebb mátrixokból (blokkokból, partíciókból) is felépíthetünk. Pl. a

utasítás az előbbi -val a

mátrixot adja.

A mátrix valamelyik elemére kerek zárójelek segítségével hivatkozunk: . Példáinknál maradva előbbi értéke 2, utóbbié . Sor- vagy oszlopvektor (azaz vagy méretű mátrix) elemére egyetlen indexszel is hivatkozhatunk. Bizonyos esetekben nincs jelentősége annak, hogy egy vektor sor- vagy oszlopvektor-e, sokszor viszont igen. Ha a vektor egyindexes elemeit külön-külön adjuk meg, akkor az sorvektor lesz. Skalár adatot pedig más programnyelvben megismert módon írhatunk, egy legáltalánosabb példa: (vagyis: ).

Ha k,h,m már definiált számértékek, akkor a vektorok felépítésére egy speciális lehetőség is rendelkezésünkre áll:

Itt egy sorvektort definiáltunk: elemei rendre lesznek. Az utolsó elem a , vagy a relációkból adódik, attól függően, hogy pozitív vagy negatív. A elhagyható, ha értéke 1. A vektor (vagy mátrix) üres is lehet.

Láttuk, hogy egy mátrixot részmátrixok segítségével is megadhatunk. A fordítottját is megtehetjük: egy mátrixból kiemelhetünk egy-egy blokkot (részmátrixot). A következő

utasítás a

mátrixot hozza létre. Ha vagy , akkor elég az egyiket megadni. Teljes sort vagy oszlopot pedig elég a (kettőspont)-tal jelölni.

Fontos alkalmazása a részmátrixoknak az

utasítás, amely az mátrix -edik sorához hozzáadja a -adik sor -szeresét.

Egy mátrix elemeiből a részmátrixhoz képest sokkal általánosabb módon is felépíthető egy új mátrix a MATLAB-ban. Nevezetesen: ha egy elemű, pedig egy elemű vektor (mindegy, hogy sor-, vagy oszlopvektor), akkor a

utasítás azt az méretű mátrixot hozza létre, amelynek elemeire teljesül. Az kifejezés értékadó utasítás jobboldalán is állhat; pl. egyetlen utasítással felcserélhetjük két sorát:

Speciális mátrixokat beépített függvényekkel is létrehozhatunk. Néhány függvényhívás:

Ezek a függvények rendre az méretű egységmátrixot, a csupa 1-esekből álló mátrixot, a zérusmátrixot, illetve véletlenmátrixot adnak. (Akármilyen méretű egységmátrixot a következőképpen értelmezünk: az azonos indexű elemei 1-esek, a többi eleme 0.) Ha , akkor a függvényeknél elég egy argumentumot írni, pl.: .

A RAND a -be eső (itt most egyenletes eloszlású) véletlenszámokat állít elő, a utasítás inicializálja a véletlenszám generátort.

Logikailag ide (kis erőltetéssel akár a részmátrixokhoz is) tartozik további három függvény. A TRIL(A) egy alsó, a TRIU(A) pedig felső háromszögmátrixot képez az mátrix megfelelő pozíciójában lévő elemeiből. A DIAG tárgyalása előtt meg kell említeni, hogy a mátrixok diagonálisaira sorszámukkal is hivatkozhatunk. Az -es mátrix fődiagonálisának sorszáma 0, a fölötte lévőké rendre , az alatta lévőké pedig . Precízebben szólva a -adik diagonálist azon elemek alkotják, amelyekre . Legyen egy hosszúságú vektor. A DIAG(v,k) azt a mátrixot hozza létre, melynek -adik diagonálisa a elemeit tartalmazza, többi eleme pedig 0. DIAG(A,k) pedig azt a oszlopvektort adja, amelynek elemei az mátrix -adik diagonálisában találhatók. Ha , akkor feltüntetése elmaradhat.

Mátrixok, vektorok, skalárok között a jelekkel az algebrában megismert műveleteket végezhetjük el és a várt eredményt kapjuk. Van egy kivétel: skalárt bármilyen méretű mátrixhoz hozzáadhatunk (kivonhatunk). Ilyenkor a skalár a mátrix minden eleméhez hozzáadódik. A ^ a hatványozás jele. Az művelet ()-szer szorozza -t önmagával, ha pozitív egész, illetve esetén az inverzet adja, ha az létezik.

A és műveletekhez hasonlóan a szorzás, osztás és a hatványozás is elemenként történik, ha a műveleti jelet a (pont) előzi meg. Pl.: eredménye lesz, míg a eredménye . Speciális mátrixművelet a transzponálás: . Alkalmas méretű mátrixok között és jelű műveleteket is használhatunk. Itt csak a legfontosabb esetet említjük. Az

eredménye az lineáris egyenletrendszer megoldása (ha az egyenletrendszer konzisztens).

8.2. Utasítások, függvények és eljárások a MATLAB nyelvben

A MATLAB rendszerben nagyon sok beépített függvény és eljárás segíti a programozói munkát. (Eljáráson itt a hagyományos függvénykiértékelésnek nem tekinthető műveletsort – pl. képernyőre rajzolás – értünk.) Mindezen alprogramok körét maga a programozó is bővítheti, így aztán egy MATLAB program általában kevés utasításból és sok-sok függvény, illetve eljárás behívásából áll.

8.2.1. Utasítások

Konkrét utasításból csak néhány fajta van. A már látott egyszerű értékadó utasításon kívül tulajdonképpen csak három további utasítást említhetünk: az IF utasítást, és két ciklusutasítást. Ezek más programnyelvből is jól ismert összetett, ún. blokkutasítások.

Az IF utasítás általános alakja:

if feltétel

    utasítások

    {elseif

        utasítások}

    <else

        utasítások>

end

A akárhányszor ismétlődhet (el is maradhat), elmaradhat.

A külön sorba írt részek egy sorba is írhatók. Ilyenkor vesszőt, vagy pontosvesszőt kell közéjük tenni. (Minden önálló utasítással is így kell eljárni.) Az IF utasítás szemantikája magyarázatot nem igényel, minimális programnyelvi ismerettel rá lehet jönni.

Feltételként logikai kifejezéseket használunk. Túlnyomórészt egyetlen reláció a kifejezés, de logikai műveletekkel össze is kapcsolhatjuk őket. Relációjelek:

Az utolsó kettő az „egyenlő“ , illetve „nemegyenlő“ jel. A logikai összeadás, szorzás és a negáció jelei rendre:

A taxatív ciklusutasítás legáltalánosabb alakja:

for

    utasítások

end

Az a ciklusváltozó, pedig általában egy vektor amely leggyakrabban (vagy ) alakú. Szemantika (legyen pl. elemszáma ): az ciklusváltozó rendre felveszi a értékeket és mindegyik mellett lefut a ciklusmag. (Ha történetesen mátrix, akkor rendre a mátrix oszlopvektorait jelenti.)

Az iteratív ciklus általános alakja:

while feltétel

    utasítások

end

A szemantika itt is ismerős.

8.2.2. Függvények

A MATLAB-ban való programozáshoz sokféle függvényt használhatunk. Egy-két függvénnyel már találkoztunk, amelyek mátrixot közvetlenül (a „semmiből“) hoznak létre, vagy egy mátrixot egyszerűen alakítanak át. Most néhány olyan függvényt sorolunk fel, amelyek már létező mátrixokból állítanak elő egy vagy több output értéket. Ezek inputjai is lehetnek természetesen skalárok. Bizonyos függvények a matematikában csak skalárokra vannak értelmezve, ezeknél a függvény a MATLAB-ban elemenként értendő. Pl. azt a -as mátrixot jelenti, amelynek minden eleme . A függvények gazdag tárházából csak a legfontosabbakat említjük meg (a HELP segít a többit is megismerni).

Elemi matematikai függvények: SIN, COS, stb. Ezek a szokásos értelmezésűek, a ROUND szabályosan kerekít.

DET(A) az determinánsát, COND(A) a mátrix kondíciószámát (''2''-es normában), RANK(A) pedig a rangját adja. Néhány további fontos függvény hívását mutatja a következő példa:

Itt az első utasítás az mátrix (vagy vektor) méretét (sorainak és oszlopainak számát) adja (a kettő közül a nem kisebbet LENGTH(A)-val közvetlenül megkaphatjuk), a következő kettő pedig a „2“ -es és „ ” normáját. Az „1“ -es és a Frobenius norma is számolható a NORM(A,1), illetve a NORM(A,'fro') hívással.

A következő függvények outputjai sorvektorok (adott esetben persze skalárak): MAX, MIN, SUM, PROD. Ha az input argumentum mátrix, de nem sorvektor, az eredményvektor az input oszlopainak a legnagyobb elemét, legkisebb elemét, összegét, szorzatát tartalmazza. Ha az input egy sorvektor, akkor a vektornak veszi a legnagyobb elemét, legkisebb elemét, összegét, és elemeinek szorzatát (az eredmény tehát itt mindenképpen skalár). Fontos tudni, hogy a MAX, és a MIN függvények két outputtal is hívhatók. Ilyenkor a második output a megfelelő elemek sorindexét, vektor esetén az elem sorszámát tartalmazza. A MAX, MIN függvények egy fontos alkalmazását érdemes példán bemutatni. Legyen

csak a vektort adná, míg a

utasításokkal a

eredményeket kapjuk. A

eredménye pedig a , , azaz legnagyobb eleme és annak oszlopindexe. Az maximális elemének sorindexét ezek után adja, ami most természetesen 2.

A beépített függvények között találunk néhány olyat is, amelyek a lineáris algebra egy-egy bonyolult algoritmusát hajtják végre. Ezek közül megemlítjük a legfontosabbakat.

Az első kettő az -felbontás két verziója, az utolsó a mátrix (aminek szimmetrikus, pozitív definitnek kell lennie) Cholesky-felbontását adja. Az eredményekre igaz a következő: (azaz ) alsó háromszögmátrixok, és felső háromszögmátrixok, permutációmátrix, továbbá

Itt az inverze (ha az létezik), ortogonális, felső háromszögmátrix és .

Az első utasítás sajátértékeit helyezi az vektorba. A második után a sajátértékek a diagonális mátrixban jelennek meg, a mátrix -edik oszlopa pedig a -hez tartozó sajátvektor (aminek „2“ -es normája 1) lesz.

Itt említjük meg a numerikus módszerek nemlineáris eljárásaira készült függvényeket is. Néhány közülük: az FSOLVE nemlineáris egyenletrendszerek megoldását végzi, a QUAD integrált számol adaptív Simpson-formulával, a SPLINE természetes interpolációt végez.

8.2.3. M-adatállományok, eljárások

MATLAB utasítások sorozata egy MATLAB programot alkot. A programot kiterjesztésű adatállományban (-adatállomány) kell elhelyeznünk. A behívása egyszerűen a nevével történik. Pl. a

utasítás betölti és lefuttatja a prog.m adatállományban lévő programot. Az egyszerű program (script-file) változói globális változók.

Lehetőségünk van egy -adatállományban a beépített függvényekkel egyező módon hívható saját függvények definiálására. Ekkor az utasítássorozatot

sorral kell kezdenünk. A kimenő és a bemenő paraméterek száma (elvben) tetszőleges. A függvény azonosítóját és az adatállomány nevét érdemes egyeztetni, azaz a függvényt az fnév.m állományban elhelyezni. Az állományban (script-file-ban is) bárhol elhelyezhetünk RETURN utasítást, ennek hatása a szokásos. A függvény definíciós programjában szereplő, nem paraméter változók lokálisak. Egy-egy változó globálissá tehető a GLOBAL utasítással. Bemenő paraméterként szerepelhet függvénynév is. Ilyenkor a function programban a FEVAL eljárást kell hívni. Pl. ha az paraméter egy kétváltozós függvény, akkor hibás, helyette a a helyes utasítás.

Ismerkedjünk meg néhány további lehetőséggel is.

A DISP(változó) vagy DISP('szöveg') eljárás a képernyőn megjeleníti az argumentumát.

A DIARY oda utasítás kiadása után minden képernyőüzenet oda (tetszőleges adatállomány lehet) is kiíródik. Ezen naplózást a utasítással fejezhetjük be, ott folytatja, ahol előzőleg abbahagyta.

Ha és egy-egy azonos hosszúságú vektor, akkor a hatására a képernyőn egy görbét kapunk az pontokból. A PLOT eljárást más formában is hívhatjuk.

Az időmérést egy példán keresztül mutatjuk be:

A CLOCK egy 6 elemű vektorban a dátum és idő adatokat tárolja, az ETIME az argumentumainak különbségét adja másodpercekben (század pontossággal), a CPUTIME pedig a MATLAB indítása óta eltelt időt méri.

9. fejezet - Előadásvázlatok fóliákon

9.1. Fólia

EGY FELADAT MEGOLDÁSÁNAK FÁZISAI

(i) Megfogalmazás, célkitűzés.

(ii) Szakmai modell felállítása.

(iii) Matematikai modell kialakítása.

(iv) Algoritmus készítés, hibaelemzés

(v) Program készítés, futtatás.

(vi) értékelés, esetleg visszacsatolás.

A numerikus módszerek alapvető feladata a (iv) fázis,

azaz már kész matematikai modellek megoldására

elfogadható hibahatáron belüli számszerű erdményeket

produkáló eljárások kidolgozása.

9.2. Fólia

MÁTRIX ÉS VEKTOR FOGALMA

Az típusú (méretű) valós mátrixon valós

számok alábbi táblázatát értjük:

Az típusú valós mátrixok halmazát jelöli,

ennek megfelelően például az azt jelenti,

hogy egy típusú (méretű) valós mátrix.

Az egyetlen sorból vagy egyetlen oszlopból álló mátrixot

vektornak nevezzük.

9.3. Fólia

TÖMÖREBB JELÖLÉSEK

Ha a sorok és oszlopok száma megegyezik, akkor a

mátrixot négyzetesnek nevezzük. Például az méretű

mátrix jelölése ekkor:

Ha egyszerűen csak vektort említünk és nem tesszük

hozzá, hogy az sor- vagy oszlopvektor-e, akkor azon

mindig oszlopvektort értünk. A komponensek kettős

indexelése ekkor felesleges (az egyik mindig 1).

Az vagy a mátrix-, illetve vektorelemekre való

hivatkozás történhet

formában is.

9.4. Fólia

MŰVELETEK

Összeadás.

Számmal való szorzás.

Transzponálás (tükrözés).

Szorzás.

Ha , akkor az szorzat neve:

Skalárszorzat. Ezzel az elemei

9.5. Fólia

MŰVELETEK FONTOS TULAJDONSÁGAI

9.6. Fólia

SPECIÁLIS MÁTRIXOK, VEKTOROK

Zérusmátrix ()

Egységmátrix, -edik egységvektor ()

Diagonálmátrix ()

Permutációmátrix: minden sorában és oszlopában

pontosan egy darab -es van és a többi elem zérus

9.7. Fólia

Háromszögmátrixok ()

Sávmátrixok (, alsó- és felső

sávszélességgel)

9.8. Fólia

MÁTRIXOK PARTICIONÁLÁSA ()

Figyelem: itt és nem egymás transzponáltja

9.9. Fólia

MÁTRIXOK DETERMINÁNSA ÉS INVERZE ()

Determináns:

Inverz mátrix:

9.10. Fólia

mátrix- (vagy vektornoma, ha ), ha

A norma szokásos jelölése:

Leggyakoribb vektornormák:

Leggyakoribb mátrixnormák ():

Mátrixok kondíciószáma:

9.11. Fólia

HIBASZÁMÍTÁSI ALAPFOGALMAK

: pontos érték,          : közelítő érték}

: az hibája,          : abszolút hibája

: abszolút hibakorlát, ha

: az relatív hibája;         

analitikus hibák: az eredeti matematikai modell és a numerikusan ténylegesen megoldott feladat közötti eltérésből adódó hiba.

öröklött hibák: az input adatok hibáinak az eredményre pontosan végzett számítások mellett gyakorolt hatásai.

kerekítési hibák: a számításokat véges sok számjeggyel végezzük, és ez a végeredményben szükségszerűen hibákat okoz.

9.12. Fólia

AZ ALAPMŰVELETEK ÖRÖKLÖTT HIBÁI

pontosak és a közelítések

és

9.13. Fólia

FÜGGVÉNYEK HIBÁI

Egyváltozós eset ()

Kondíciószám:

Többváltozós eset}

(, , , )

Kondíciószám ():

9.14. Fólia

LEBEGŐPONTOS HIBASZÁMÍTÁS

Ha , akkor

,         kerekítés:

az egységnyi kerekítés; a gépi epszilon:

A műveletekre is feltesszük, hogy

9.15. Fólia

ÉRZÉKENYSÉG, STABILITÁS

Direkt hibák, inverz hibák

Legyenek és a pontos, és pedig a közelítő értékek.

direkt hiba:

inverz hiba: , ha

9.1. ábra - Direkt hiba, inverz hiba (Direkt hiba, inverz hiba (dx=\Delta x , dy=\Delta y ), Direkt hiba, inverz hiba (dx=\Delta x , dy=\Delta y ))

Direkt hiba, inverz hiba (dx=\Delta x , dy=\Delta y )

Az inverz hiba elemzése: inverz hibaanalízis

Hibaszámítási „ökölszabály”

9.16. Fólia

FELADATÉRZÉKENYSÉG

Pontos számítás mellett mennyire érzékeny a végeredmény az input adatok változására(perturbációjára).

Egyik jellemzője a kondíciószám.

Algoritmusok numerikus stabilitása

A kerekítési hibákra való érzékenysége (vagy éppenséggel a kevésbé érzékenysége, robusztussága)

Példa:

Ha és , akkor az

rekurzióval definiált sorozatban és .

ám, ha

és

,

akkor és az .

A rekurziós képlettel, duplapontos aritmetikában az

eredeti és kezdőértékekkel nem a pontos, hanem éppen az és az eredményeket kapjuk.

Tehát

Itt a direkt hibák:

igen nagyok, pedig az input adatok megváltozásai, a

rendkívül kicsinyek.

Ha az , értékeket tekintjük pontosnak, akkor a , mennyiségek az inverz hibák.

Látható a feladatérzékenység és a kerekítési hibák kölcsönhatása: a kerekítési hibák hatása olyan, mintha a perturbált input adatokkal pontosan számolnánk.

A kerekítési hibák halmozódását a következő ábrasorozat mutatja (a vízszintes tengelyen az , a függőlegesen pedig az (kék színnel) és az (piros színnel) látható.

(Az algoritmusok stabilitására példát a lineáris egyenletrendszereknél mutatunk be.)

9.17. Fólia

9.2. ábra - Kerekítési hibák halmozódása, I.

Kerekítési hibák halmozódása, I.

9.3. ábra - Kerekítési hibák halmozódása, II.

Kerekítési hibák halmozódása, II.

9.4. ábra - Kerekítési hibák halmozódása, III.

Kerekítési hibák halmozódása, III.

9.18. Fólia

9.5. ábra - Kerekítési hibák halmozódása, IV.

Kerekítési hibák halmozódása, IV.

9.6. ábra - Kerekítési hibák halmozódása, V.

Kerekítési hibák halmozódása, V.

9.7. ábra - Kerekítési hibák halmozódása, VI.

Kerekítési hibák halmozódása, VI.

9.19. Fólia

LINEÁRIS EGYENLETRENDSZEREK

( egyenlet és ismeretlen)

Ha

akkor tömör alakban:

Ha , akkor az egyenletrendszer

Ha

akkor a megoldás az alábbi hipersíkok közös része:

Három eset lehetséges:

(i) nincs megoldás

(ii) pontosan egy megoldás van

(iii) végtelen sok megoldás van

A továbbiakban feltesszük hogy . Ekkor

Cramer-szabály:

A Cramer-szabály alkalmazása költséges, más utat kell keresni.

A költséget a flopok számával mérjük.

9.20. Fólia

Háromszögmátrixú egyenletrendszerek

Sémájuk:

9.8. ábra - Alsó háromszögmátrixú egyenletrendszer sémája

Alsó háromszögmátrixú egyenletrendszer sémája

9.9. ábra - Felső háromszögmátrixú egyenletrendszer sémája

Felső háromszögmátrixú egyenletrendszer sémája

Ez utóbbi kirészletezve:

Ha , akkor , tehát a

megoldó algoritmus:

A ciklus magját írhatjuk így is:

Ezen algoritmus költsége a mai számítógépeken nagy méret esetén is alacsony:

flop.

Cél tehát az eredeti egyenletrendszer ekvivalens (megoldásőrző) átalakítása háromszögmátrixúvá, elfogadható költséggel.

9.21. Fólia

A Gauss-módszer

Ekvivalens (megoldásőrző) átalakítás:

 egyik egyenletből kivonjuk egy másik egyenlet valahányszorosát 

-adik lépés:

Ha a -adik egyenlet alattiakból kivonjuk a -adik egyenlet

szorosát, úgy az -edik egyenletből is kiiktatjuk

(elimináljuk) a -adik ismeretlent:

A nyilván kimarad, de a is, hisz a későbbiekben mindegy, mi áll a főátló alatt.

Fenti algoritmus középső sora az ismert jelöléssel:

A Gauss-módszer MATLAB programja

  1  function [x]=Gauss(A,b) 
  2  n=length(b); 
  3  % 
  4  %eliminációs fázis: 
  5  % 
  6  for k=1:n-1 
  7    for i=k+1:n 
  8       g=A(i,k)/A(k,k); 
  9       A(i,k+1:n)=A(i,k+1:n)-g*A(k,k+1:n); 
 10          b(i)=b(i)-g*b(k); 
 11    end % i 
 12  end % k; 
 13  % 
 14  %vége az eliminációs fázisnak 
 15  %visszahelyettesitő fázis: 
 16  [x]=fhmatrixszal(A,b) 
 17  return 
     
  0  % a háromszögmátrixú e.r.-t megoldó program 
  1  function [x]=fhmatrixszal(A,b) 
  2  n=length(A); 
  3  x=zeros(n,1); 
  4  for i=n:-1:1 
  5    x(i)=(b(i)-A(i,i+1:n)*x(i+1:n))/A(i,i); 
  6  end 
  7  return 

A Gauss-módszer műveletigénye: flop.

A Gauss-elimináció determináns számításra is alkalmas.

9.22. Fólia

Főelemkiválasztás

Ha a pivotelem, azaz , akkor szükséges a sorok (és/vagy oszlopok) cseréje. De a numerikus stabilitás érdekében legyen mindenkor lehetőleg nagy abszolút értékű a pivotelem.

Részleges főelemkiválasztás: A -adik lépésben a -adik oszlop () elemei közül kiválasztjuk a maximális abszolút értékűt, majd annak sorát a -adik sorral felcseréljük. A pivotálás után teljesül, hogy

Teljes főelemkiválasztás: A -adik lépésben az () mátrixelemek közül választjuk ki a maximális abszolút értékűt. Ha ennek indexpárja (), akkor a -adik és a -edik sort, valamint a -adik és az -edik oszlopot felcseréljük. A pivotálás után teljesül, hogy

Figyelem: oszlopcsere esetén változócsere is történik.

9.23. Fólia

A pivotelemek növekedési tényezője

A Gauss-módszer eliminációs fázisában egy

ekvivalens egyenletrendszerekből álló sorozatot képezünk, ahol

A végső állapotban

ahol a -adik fő-, vagy pivotelem. A pivotelemek növekedési tényezője:

befolyásolja az eljárás numerikus stabilitását.

Példa. Az egyenletrendszer jól kondicionált (, kicsi), a pontos megoldás: , , bináris duplapontos aritmetikában: , . (Kézi számítással is tanulságos végigszámolni a 16 jegyű decimális aritmetikában.) Főelemkiválasztás nélkül igen nagy és a nagyon rossz , adódik, míg részleges főelemkiválasztással és lehető legjobb , közelítő megoldást kapjuk.

A Gauss-módszer főelemkiválasztással numerikusan nem stabil, míg főelemkiválasztással az.

9.24. Fólia

Gauss-Jordan elimináció részl. főelemkiválasztással

egyenletrendszer:

  1  function [X]=Gauss_Jordan(A,B) 
  2  [n,m]=size(B); r=n+m; 
  3  A=[A,B]; 
  4  % I.1 fazis: foatlo alatti eliminacio 
  5  for k=1:n-1 
  6    for i=k+1:n 
  7       [amax,ind]=max(abs(A(k:n,k))); 
  8          t=ind+k-1; 
  9  %      ind: az A(k:n, k)-beli index, t az A-beli 
 10       if t k 
 11  %         a k-adik es a t-edik sor csereje: 
 12          A([k,t], :)=A([t,k], :); 
 13       end % if 
 14       g=A(i,k)/A(k,k); 
 15       A(i,k+1:r)=A(i,k+1:r)-g*A(k,k+1:r); 
 16    end % i 
 17  end % k; 
 18  % I.2. fazis: foatlo folotti eliminacio 
 19  for k=n:-1:1 
 20  % A kovetkezo i ciklus k=1 eseten ures 
 21    for i=1:k-1 
 22       g=A(i,k)/A(k,k); 
 23       A(i,n+1:r)=A(i,n+1:r)-g*A(k,n+1:r); 
 24    end % i 
 25    X(k,1:m)=A(k,n+1:r)/A(k,k); 
 26    end % k; 
 27  % II. fazis itt nincs: barmely j-re    () 
 28  % az Ax=B(:, j) egyenletrendszer megoldasa: X(:, j) 
 29  % pl. B=I eseten X=inv(A) 
 30  return 

Műveletigény, ha n=m (pl. invertálás): flop.

9.25. Fólia

Az -felbontás

A=LU, ahol L alsó-, U felső háromszögmátrix:

9.10. ábra - Az Az LU -felbontás sémája-felbontás sémája

Az LU -felbontás sémája

Egy nemszinguláris mátrixnak akkor és csak akkor létezik -felbontása, ha összes főminor mátrixa is nemszinguláris, azaz

Gauss-eliminációval:

azaz: oszlopai osztva a főátlóbeli elemmel

()

9.26. Fólia

Egyenletrendszer megoldása -módszerrel

Minden nemszinguláris mátrixhoz létezik olyan permutációmátrix, hogy a mátrixnak van -felbontása.

Az -módszer algoritmusa

1. Határozzuk meg a felbontást

(így: adódik)

2. Oldjuk meg az egyenletrendszert -ra

3. Oldjuk meg az egyenletrendszert -re

Cholesky-felbontás

Ha az mátrix szimmetrikus és pozitív definit, akkor létezik olyan -felbontása, hogy (ez a Cholesky-felbontás).

Műveletigénye: flop.

9.27. Fólia

Iteratív módszerek

Legyen az pontos megoldása , az pedig az -hoz konvergáló sorozat ().

Jacobi-iteráció

Skalár alakban:

Seidel-iteráció

MATLAB program a Seidel-iterációra

  1  function [x, iszam]=Seidel(B, c, x0, q, tol, imax) 
  2    n=length(c); 
  3    p=q/(1-q); 
  4  %   hogy legalabb egy iteraciot biztos vegrehajtson: 
  5    d=inf; 
  6    iszam=0; 
  7    x=x0; 
  8    while (d>tol) & (iszam<imax) 
  9       iszam=iszam+1; 
 10       for i=1:n 
 11          x(i)=B(i,:)*x; 
 12       end % i 
 13       d=p*norm(x-x0, inf); 
 14       x0=x; 
 15    end % while 
 16    return 

Konvergencia-tétel

Ha , akkor mind a Jacobi-, mind a Seidel-iteráció val előállított sorozat tetszőleges választása esetén az egyetlen megoldáshoz konvergál. Érvényes továbbá a következő hibabecslés:

Iteratív módszerek előnyei

(i) a kerekítési hibák nem halmozódnak

(ii) egzakt hibabecslés, ha a feltételek teljesülnek

(iii) adott esetben kisebb műveletigény

            (egy lépés csak flop)

9.28. Fólia

Lineáris egyenletrendszerek hibabecslése

elméleti megoldás: , a közelítő megoldás:

direkt hiba:

reziduális hiba:

inverz hibák: ,

Inverz hiba a jobboldali vektorban

Ha nemszinguláris és , akkor

Inverz hiba az együtthatómátrixban

Ha nemszinguláris, és , akkor

Egyidejű inverz hiba a jobb- és baloldalon

Ha nemszinguláris, és , akkor

A , és a kondíciószám becsléséhez utalunk a szabadon hozzáférhető LINPACK programcsomagra.

Az inverz hibákra több helyen ajánlják a , választást, ahol az abszolútérték elemenként értendő, pedig a gépi epszilon.

A mindig használható.

A inverz hiba becslése Wilkinson tételével

Az egyenletrendszer Gauss-módszerrel lebegőpontos aritmetikában kapott közelítő megoldása kielégíti az

egyenletrendszert, ahol

és a pivot elemek növekedési tényezője, pedig az aritmetika egységnyi kerekítése.

A relatív direkt hiba becslése Wilkinson tételével

A tételből:

Kis kondíciószám esetén feltehetjük, hogy

és ezzel

Utólagos hibabecslés Auchmuty tételével

ahol a (konstansnak nevezett) az -tól és az () hibavektor irányától függ, továbbá .

A kísérletileg becsülhető, tapasztalat szerint többnyire .

9.29. Fólia

Iteratív javítás

Mivel a közelítést határozhatjuk meg, így reményeink szerint is csupán javíthatjuk az közelítést.

Algoritmusa:

Legyen

for

    

    Számítsuk ki az közelítő megoldását

    

end

Az egyenletrendszerek megoldására az -módszer előnyös.

Tapasztalatok szerint sem a , sem az nem csökken monoton, ezért az irodalomban kilépési feltételként ajánlott

nem mindig megbízható.

Helyette néhány lépés megtételét, és a kapott közelítések közül a legkisebb reziduális hibájú választását javasoljuk.

Irodalomjegyzék

[1] Bahvalov, N. Sz.. A gépi matematika numerikus módszerei. Műszaki Könyvkiadó, Bp.. 1977.

[2] Bjezikovics, Ja. Sz.. Közelítő számítások. Tankönyvkiadó. 1952.

[3] Dr. Bálint Elemér. Numerikus és grafikus közelítő módszerek (Műszaki matematikai gyakorlatok, szerk. Dr. Fazekas Ferenc). Tankönyvkiadó, Budapest. 1978.

[4] Égertné Molnár É.-Kálovics F.-Mészáros J.-né. Numerikus Analizis. Miskoci Egyetemi Kiadó. 1992.

[5] Forsythe, G. E., Moler, C. B.. Lineáris algebrai problémák megoldása számítógéppel. Műszaki Könyvkiadó. 1976.

[6] Galántai, A., Jeney, A.. Numerikus módszerek. Miskolci Egyetemi Kiadó. 2005.

[7] Galántai, A., Jeney, A., Mészáros, J.-né. Numerikus módszerek példatár. (szerk. alatt).

[8] Golub, G. H., Van Loan, C. F.. Matrix Computations (second edition). The Johns Hopkins University Press, Baltimore. 1993.

[9] Iványi A. (szerk.). Informatikai algoritmusok 1. TT.. ELTE Eötvös Kiadó. 2004.

[10] P. Henrici. Numerikus analízis. Műszaki Könyvkiadó, Bp.. 1985.

[11] Higham, N. J.. Accuracy and Stability of Numerical Algorithms. SIAM, Philadelphia. 1996.

[13] Jennings, A., McKeown, J. J.. Matrix Computation, (second edition). John Wiley & Sons. 1992.

[14] Móricz F.. Numerikus módszerek az algebrában és analízisben. Polygon. 1997.

[15] Moler, C. B.. Numerical Computing with MATLAB. SIAM, Philadelphia. 2004.

[16] Overton, M. L.. Numerical Computing with IEEE Floating Point Arithmetic. SIAM, Philadelphia. 2001.

[17] Popper Gy., Csizmás F.. Numerikus módszerek mérnökök. Akadémiai Kiadó-Typotex. 1993.

[18] Ralston, A.. Bevezetés a numerikus analízisbe. Műszaki Könyvkiadó. 1969.

[19] Rice, J. E.. Numerical Methods, Software, and Analysis. McGraw-Hill. 1983.

[20] Rice, J. E.. Matrix Computations and Mathematical Software. McGraw-Hill. 1983.

[21] Rózsa P.. Lineáris algebra és alkalmazása. Műszaki Könyvkiadó. 1974.

[22] Stewart, G. W.. Matrix Algorithms, Vol. I: Basic Decompositions. SIAM, Philadelphia. 1998.

[23] Stewart, G. W.. Matrix Algorithms, Vol. II: Eigensystems. SIAM, Philadelphia. 2001.

[24] Stoyan G., Takó G.. Numerikus módszerek 1. ELTE-Typotex. 1993, 1995, 1997.

[25] Stoyan G. (szerk.). Matlab (frissített kiadás). Typotex Kiadó, Budapest. 2005.

[26] Szamarszkij, A. A.. Bevezetés a numerikus módszerek elméletébe. Tankönyvkiadó. 1989.

[27] Szidarovszky Ferenc. Bevezetés a numerikus módszerekbe. Közgazdasági és Jogi Kiadó. 1974.

[28] Ueberhuber, C. W.. Numerical Computation 1-2 (Methods, Software, and Analysis). Springer. 1997.

[29] Watkins, D. S.. Fundamentals of Matrix Computations. John Wiley & Sons. 1991.

[30] Wilkinson, J. H.. Rounding Errors in Algebraic Processes. Dover. 1994.