7. fejezet - Numerikus módszerek

Ábrázoljuk a függvényt, és a grafikon alapján elkülöníthetjük a gyököket tartalmazó intervallumokat. Mindez valós gyökök esetén – főleg a mai számítástechnikai eszközökkel – roppant egyszerű. Ha komplex gyököket keresünk, akkor az f(z)=0 egyenletre alkalmazzuk a következő módszert:

Olyan módszerről, amely egy tetszőleges függvény esetén megadja a gyököket tartalmazó intervallumot, nem tudunk. Általában a következőket tesszük:

Talán ez a legegyszerűbb módszer. Előnye, hogy mindig használható, a gyök tetszőleges pontosságú megközelítésére biztosít lehetőséget. Hátránya a közelítés sebessége.

Legyen adott az f(x)=0 egyenlet, ahol f(x) folytonos függvény, és a ₁ valamint b ₁ , ahol f(a ₁ )‑nek és f(b ₁ )-nek ellentétes az előjele, azaz f(a ₁ )⋅f(b ₁ )<0. Vegyük fel c ₁ -et az a ₁ , b ₁ intervallum felező pontjaként, azaz:

A két félintervallum közül válasszuk azt, ahol teljesül a Bolzano-tétel, azaz

Az algoritmus tetszőlegesen sokáig folytatható, így egyre jobb közelítést nyerhetünk. Ha például 10 lépésben közelítünk, akkor az eredeti intervallum hosszának már csak 2¹⁰-en, azaz 1024-ed része lesz az új intervallum hossza. Másképpen, mivel 1024≈10³, ezért kb. hárommal több pontos tizedesünk lesz, tehát minden újabb pontos tizedes jegy meghatározásához kb. három iterációs lépésre van szükség. Természetesen ennél gyorsabb módszerek is léteznek, de akkor az f(x) függvénynek bizonyos követelményeknek meg kell felelnie.

A példa megoldását táblázatkezelővel Táblázat 7.4 tartalmazza, míg az eredményeket a (Táblázat 7.5) táblázatban tanulmányozhatjuk.

Az algoritmus lépései sokban hasonlítanak az előzőekben leírt intervallum-felezéshez, csak a c _k pontok meghatározása történik másként. Most az a ₁ ,f(a ₁ ) valamint a b ₁ ,f(b ₁ ) pontokon át fektessünk egy egyenest (ami a függvénynek egy „húrja”), és ennek az x tengellyel való metszéspontja legyen c ₁ . Ezt az értéket úgy is tekinthetjük, mint a gyök egy közelítő értékét, bár tudjuk, hogy az érték „hamis”, azaz „false” (innen a módszer elnevezése is). Ha az eredményt nem tarjuk elég pontosnak, akkor a két intervallum közül az előzőekben leírtak szerint választunk, majd egy újabb húr segítségével pontosítunk. A módszer tetszőlegesen sokáig folytatható, tehát a közelítés pontossága határtalanul fokozható, csak kellően türelmesnek kell lennünk. Egy idő után általában már csak az intervallum egyik oldalán lesz változás, ettől a húrok alig változnak, tehát a közelítés „lelassul”. Számos módszer van, amely ezen a problémán próbál segíteni (például minden n-edik lépésnek egy intervallum-felezést iktatunk be).

Nézzük most az előzőekben szövegesen leírt módszert egyenletekkel, a 7.4. ábra jelöléseit használva! Írjuk fel a két ponton átmenő egyenes egyenletét

majd határozzuk meg a metszéspontot, azaz y=0 helyen x=c ₁ értékét:

A következő lépés az új intervallum meghatározása. Eljárhatunk az intervallum-felezésnél ismertetett módon is, de néha alkalmazhatunk feltételek nélküli rekurziót is. Abban az esetben, ha az -ben – tehát nincs inflexiós pont, azaz a görbe homorú, vagy domború – az új intervallum meghatározása egyértelmű. Nézzük például az ábrán látható esetet! Ekkor az egyenes mindig a görbe fölött fog haladni, tehát a tengelyt mindig „előbb” metszi, mint a görbe, azaz az intervallumnak mindig az „ a ” vége fog mozogni. Jelöljük ekkor b ₁ -et c _-1 -gyel, a ₁ -et c ₀ -val:

Mivel az új intervallumok esetén c ₀ =c ₁ , majd c ₁ =c ₂ helyettesítéseket alkalmazunk, ezért a rekurziós összefüggés a következőképpen írható fel:

Kérdés most, hogy mi a teendő akkor, ha a görbe nem ilyen helyzetű. A rekurziós módszer mindig alkalmazható, ha . Ekkor a -1-es indexűnek azt a pontot kell választani, amelynek előjele megegyezik a második derivált előjelével.

A példa megoldása táblázatkezelővel a (Táblázat 7.6) táblázatban látható.

A B3 és a C3 cellák tartalma megegyezik az intervallum-felezésnél leírtakkal. Az eredményeket a Táblázat 7.7 tartalmazza.

Mint látható, sokkal gyorsabban közelítette a gyököt, mint az intervallum-felezés, bár ezt általánosan nem lehet garantálni!

Az eddigi módszerek minden esetben a gyököt közelítették, de néha kissé lassan. A következőkben két olyan módszert nézünk meg, amelyek nem mindig konvergensek, de vagy nagyon könnyen használhatók, vagy nagyon gyorsan adnak eredményt.

Az első módszer az iteráció. Az egyenlet sokszor f(x)=0 alakban adott, és a gyököt az [a,b] intervallumban keressük. Ilyen esetben a következőképpen kell átalakítanunk az egyenletet:

Most képezzük a következő sorozatot:

Az így kapott x₀, x₁, …, x _n számsorozat (ha konvergál) az f(x) egyenlet gyökéhez tart, hiszen ha x_∞ a sorozat határértéke, akkor:

A következőkben nézzük meg grafikusan, milyen esetek fordulhatnak elő!

Az előző ábrákból is leolvasható, hogy konvergenciára csak akkor számíthatunk, ha g(x) meredeksége kisebb, mint 45°, azaz abban az intervallumban, ahol a gyököt keressük. Az egyszerű iterációt előszeretettel használják, mert rendkívül könnyen alkalmazható akár egy táblázatkezelővel is.

A mintapéldánk (7.14. egyenlet) esetén az iteráció nem működne, mert a derivált értéke:

A gondot az okozza, hogy a függvényünk túl „meredek”. Ha egy alkalmas konstanssal megszorozzuk, akkor a zérus helye nem változik, de a meredeksége igen, hiszen:

Nézzük, milyen feltétel van f-re:

Legyen , ahol igaz, hogy , azaz M az első derivált maximális, vagy – negatív esetben – minimális értékével egyezik meg az adott intervallumon. Most, ha , akkor a konvergencia feltétele:

ami mindig teljesül, hiszen M nagyobb, vagy egyenlő, mint a derivált, ezért a tört értéke egynél kisebb! Írjuk fel még az így kapott iterációs alakot:

Nézzük először az eredeti példa megoldását képletekkel (Táblázat 7.8), majd pedig értékekkel (Táblázat 7.9)!

Ahogy azt vártuk, az x _i sorozat divergens. Mivel az függvény maximuma akkor lesz, ha a szinusz értéke -1, így az M=2. Ebből α-ra -0,5 adódik. Az átalakított egyenletet tartalmazza a Táblázat 7.10, valamint a megoldást értékekkel szemlélteti a Táblázat 7.11.

Mint látható, a konvergencia viszonylag gyors.

Legyen most adott az f(x)=0 egyenlet, és tudjuk, hogy az [a,b] intervallumban van gyöke! Legyen , és helyettesítsük a függvényt az x ₀ -nál húzott érintőjével. Az érintő és az x tengely metszéspontját (x ₁ ) tekintsük a gyök egy közelítésének! Az eljárás folytatásaként tegyünk mindent ugyan így, csak most az x ₁ értékkel! Írjuk fel az egyenleteket és az összefüggéseket a 7.11. ábra alapján! Az egyenes egyenlete:

vagy másképpen:

Az x tengellyel való metszéspontban y=0, így x kifejezhető:

Az eljárás folytatásaként írhatjuk, hogy:

Vegyük észre, hogy ez az összefüggés és a 7.29 egyenlet nagyon hasonlít egymásra, az ott szereplő M tulajdonképpen az itteni nevezőben lévő derivált maximális értéke! Gyanítható ezek után, hogy az érintő módszer konvergenciájával nem lesz gond. Ez általában így is van. Ha elég közel vesszük fel a kezdeti x₀ pontot a gyökhöz, akkor majdnem mindig konvergens a közelítés. A gond az, hogy általában fogalmunk sincs, hogy hol a gyök, örülünk, ha egy intervallumra sikerül leszűkíteni a keresést! Alkalmazása során általában meg szoktak adni egy maximális iteráció számot is.

Miért használják mégis? Például azért, mert komplex gyökök meghatározására is alkalmas, vagy azért, mert ha konvergens, akkor majdnem mindig „másodrendűen” gyors, azaz:

ahol x a gyök, C pedig egy állandó érték. Ez azt jelenti, hogy lépésenként általában megduplázódik a helyes jegyek száma a gyökközelítő értékében. Nézzük meg most a megszokott példánkat a (Táblázat 7.12 és Táblázat 7.13) táblázatokban!

A (Táblázat 7.13) táblázatban szereplő számsor mindennél jobban illusztrálja a konvergencia sebességét! Milyen problémák léphetnek fel mégis, mikor lehet gond? Például, ha az érintő mentén kilépünk az [a,b] intervallumból (7.12. ábra), vagy szerencsétlenül választjuk meg a kezdőértéket: , és . Ekkor , és így . A számítást tovább folytatva: . Látható, hogy körbe-körbe járunk (7.13. ábra)!

A konvergenciával kapcsolatban azonban nemcsak rosszat mondhatunk. Azt várnánk, hogy nem működik az eljárás, ha a gyök például egy minimumhely is egyben, hiszen ott a derivált értéke zérus. A deriválttal pedig osztunk! Nem szabad elfelejtenünk, hogy a gyököt csak közelítjük, de sosem érjük el, így a derivált sem lesz nulla a közelítés közben. Az [a,b] intervallumban lehet szélsőértékhely, de akkor könnyen úgy járhatunk, hogy kilépünk az intervallumból, hasonlóan, mint az a (7.13. ábra) ábrán is látható. Arra mindenképpen ügyeljünk, hogy a közelítést ne indítsuk minimum-, vagy maximumhelyről! Az előzőek illusztrálására nézzük meg, hogy az függvény gyökét hogyan közelíti a módszer (Táblázat 7.14)!

Mint látható, a konvergencia itt is megvan, csak sokkal lassabb. Bizonyítható, hogy többszörös gyökök esetén a gyököt közelítése csak lineáris.

A módszer egy utolsó érdekessége, hogy segítségével könnyedén készíthetünk gyökvonó algoritmust. Legyen az egyenletünk az . Ennek zérus helye nyilvánvalóan: . Ekkor a közelítő formula:

Speciálisan, ha k=2, akkor négyzetgyökvonásra a következő összefüggést kapjuk:

Az egyik leggyakoribb tevékenység, amit egy polinommal el kell végeznünk, hogy meghatározzuk egy adott helyen a helyettesítési értékét. A Horner-séma ennek a számításnak a műveletek számára optimalizált eljárása, azaz itt kell a legkevesebbet számolnunk. Gyakorlatilag csak egy átrendezésről van szó.

Legyen az egyenletünk

alakú.

Ekkor átrendezhetjük a következő alakra:

Példaként nézzük azt az egyenletet, amely a (7.14. ábra) ábrán is szerepel:

Valós x esetén a következő függvény valósítja meg az eljárást:

Function BehelyValos(x AsDouble, n AsInteger, a() AsDouble) AsDouble
' a() az együttható vektor
' n a polinom fokszáma
Dim i AsInteger, y AsDouble
y = a(n) ' a legmagasabb fokú tag együtthatója
For i = n - 1 To 0 Step -1
If a(i) <> 0 Then         ' ellentetes elojellel megegyeznek
  If Abs(1 + y * x / a(i)) < 0.0000000000001 Then
    y = 0 ' a numerikus hibák csökkentésére
Else
    y = y * x + a(i)
End If
Else
  y = y * x ' mert a(i)=0
End If
Next i
BehelyValos = y
End Function

Komplex esetben eljárást kellett írni, mert két visszaadott érték van. A program listája a következő:

Sub BehelyKomplex0(xr AsDouble, xi AsDouble, n AsInteger, a() AsDouble,_
ByRef yr AsDouble, ByRef yi AsDouble)
' a() az együttható vektor
' n a polinom fokszáma
' xr x reális része; xi x imaginárius része
' y az eredmény
' yr y reális része; yi y imaginárius része
Dim i AsInteger, y AsDouble
yr = a(n)
yi = 0
For i = n - 1 To 0 Step -1
If a(i) <> 0 Then
If Abs(1 + (yr * xr - yi * xi) / a(i)) < 0.0000000000001 Then
    y = 0
Else
    y = yr * xr - yi * xi + a(i)
End If
Else
  y = yr * xr - yi * xi
End If
If yi * xr <> 0 Then
If Abs(1 + yr * xi / yi / xr) < 0.0000000000001 Then
    yi = 0
Else
    yi = yr * xi + yi * xr
End If
Else
  yi = yr * xi + yi * xr
End If
yr = y
Next i
End Sub

Az eljárás teljes mértékben a megszokott kézi módszert követi, ami tulajdonképpen az általános iskolákban tanított osztás általánosítása. Ismétlésként álljon itt egy példa:

Az osztás eredménye, és a maradék (18) is kiadódott. Az ezt megvalósító programkód pedig a következő:

Sub PolDiv(n As Integer, m As Integer, _
ByRef a() As Double, ByRef b() As Double, ByRef c() As Double)
' a/b=c, ahol a, b és c a polinomok együtthatóit tartalmazó vektorok
' a 0. fokú tag a vektor 0 elemében van
' ha van maradék, akkor a-ban lesz
Dim i As Integer, j As Integer
ReDim c(n - m)
For i = n - m To 0 Step -1
c(i) = a(i + m) / b(m)
For j = 0 To m
If j = 0 Then
a(i + m - j) = 0
Else
If a(i + m - j) <> 0 Then
If Abs(1 - c(i) * b(m - j) / a(i + m - j)) < 0.0000000001 Then
a(i + m - j) = 0 ' a különbség közel nulla
Else
a(i + m - j) = a(i + m - j) - c(i) * b(m - j)
End If
Else
a(i + m - j) = -c(i) * b(m - j) ' mert a(i+m-j)=0
End If
  End If
Next j
Next i
End Sub

Az egyszerűen, szinte ránézésre meghatározható gyökökről van szó ebben a fejezetben. Ilyen például, ha a 7.37. egyenletben a₀ értéke zérus. Ekkor biztosan van olyan gyök, amely nulla, és a polinom fokszáma csökkenthető. A későbbi eljárásoknál fontos lesz, hogy az egyenlet nulladfokú tagja nem lehet zérus, mert osztunk vele!

Általában ritkán fordul elő, hogy egy polinomnak az 1 gyöke, de könnyen ellenőrizhető, és lehet, hogy a későbbiekben a közelítések során nem tudnánk pontosan előállítani. Az 1 akkor gyök, ha a polinom együtthatóinak összege zérus, hiszen behelyettesítéskor minden x hatvány helyére 1 kerül! Ha az 1 gyök, akkor az egyenletet osszuk el gyöktényezővel!

Hasonló a helyzet a -1-el is, de ekkor a páros kitevőjű tagok együtthatóinak összege fog megegyezni a páratlan kitevőjű tagok együtthatóinak összegével, hiszen a párosaknál a behelyettesítés eredményeként az x hatvány 1 lesz, a páratlanoknál pedig -1. Ha a -1 gyök, akkor az egyenletet osszuk el gyöktényezővel!

Természetesen ne feledkezzünk meg arról, hogy mindegyik primitív gyök többszörös gyök is lehet, tehát az eljárást addig kell folytatni, amíg mindet meg nem kapjuk. Ha szerencsénk van, a fokszám csökkenés miatt nem is kell közelítő módszert alkalmaznunk.

A pontosság miatt (mivel a közelítő eljárások csak „közelítenek”) hacsak lehet, szeretjük az egyenlet megoldását valamilyen képletből megkapni. Az első- és a másodfokú esetet mindenki ismeri, de a harmadfokút már nem. Ennek megoldását – ugyan nem ő találta ki, csak ő publikálta – Gerolamo Cardanoról nevezték el. Van azonban valami, ami nem alaptalanul viseli az ő nevét: a kardántengely.

Ennyi technika történet után térjünk vissza a matematikához! Az általános harmadfokú egyenlethez nem létezik megoldó képlet, csak olyan változatához, amelyben nem szerepel másodfokú tag:

Az általános egyenletből azonban ez az alak mindig előállítható a teljes négyzethez hasonló eljárással, csak itt természetesen köbről van szó. Ekkor az átszámítás a következő módon végezhető el:

Az (7.41) egyenlet megoldása:

ahol

Mint látható, ezt a számítást elvégezni nem olyan egyszerű, mint a másodfokú egyenlet gyökeit meghatározni. Érdekeség talán, hogy abban az esetben, ha a megoldás három valós gyökből áll, már u és v meghatározása során is komplex számokat kapunk! Sokszor nem is bajlódnak vele, inkább közelítő módszert alkalmaznak. Mi nem így tettünk, itt a programkód:

Sub harmadfoku(s As Integer, ByRef a() As Double)
' s annak a sornak a száma, ahová a gyököt ki kell írni
' a az együttható vektor
Dim p As Double, q As Double, u As Double, v As Double, diszkr As Double
Dim uabsz As Double, fi As Double
Ifa(3) = 0 Then
Call masodfoku(s, a)
Else ' az egyenlet átalakítása z^3+p*z+q=0 -ra
kesz = True
p = (3 * a(1) - a(2) * a(2)) / 3
q = a(0) - a(2) * a(2) * a(2) / 27 - a(2) / 3 * p
If p = 0 Then ' z^3+q=0
  p = Sgn(-q) * (Abs(q)) ^ (1 / 3)
' gyökök kiírása az eredménylapra
Range("F" + Szam(s)).Value = p - a(2) / 3
s = s + 1
Range("F" + Szam(s)).Value = -p / 2 - a(2) / 3
Range("G" + Szam(s)).Value = 3 ^ 0.5 * p / 2
s = s + 1
Range("F" + Szam(s)).Value = -p / 2 - a(2) / 3
Range("G" + Szam(s)).Value = -3 ^ 0.5 * p / 2
Exit Sub
End If
If q = 0 Then ' z^3+p*z=0 => z=0 => x=-B/3/A, de A=1
Range("F" + Szam(s)).Value = -a(2) / 3
a(0) = a(1) - 2 * a(2) * a(2) / 9
a(1) = 2 / 3 * a(2)
a(2) = 1
Call masodfoku(s + 1, a)
Exit Sub
End If
diszkr = q * q / 4 + p * p * p / 27
If diszkr < 0Then ' a valós gyök meghatározása
fi = Sqr(Abs(q * q / 4 + p * p * p / 27))
  uabsz = (q * q / 4 + fi * fi) ^ (1 / 6)
fi = (1 - Sgn(-q)) / 2 * 3.1415926536 + Atn(fi / (-q / 2))
  u = 2 * uabsz * Cos(fi / 3) - a(2) / 3 / a(3)
Else ' diszkr > 0, az egyetlen valós gyök meghatározása
fi = -q / 2 + Sqr(diszkr)
  u = Sgn(fi) * Abs(fi) ^ (1 / 3)
  v = -p / 3 / u
  u = u + v - a(2) / 3
End If
v = Round(u) ' egész megoldások pontosítása
If Abs(((u + a(2)) * u + a(1)) * u + a(0)) > _
Abs(((v + a(2)) * v + a(1)) * v + a(0)) Then u = v
Range("F" + Szam(s)).Value = u
a(0) = a(1) + u * (u + a(2))
a(1) = u + a(2)
a(2) = 1
Call masodfoku(s + 1, a)
End If 'harmadfokú-e
End Sub

Több oka is van annak, hogy szeretnénk előállítani egy olyan polinomot, amelynek az eredeti polinom minden gyöke a megoldása, egyszeres multiplicitással. Az egyik ok, hogy ez nyilvánvalóan fokszámcsökkenéssel jár, és lehet, hogy nem is lesz szükségünk közelítő (azaz pontatlanabb) módszer alkalmazására. A másik ok, hogy – a későbbiekben ismertetésre kerülő – a valós gyökök számát meghatározó módszer csak ilyen esetben működik.

Legyen adott az f(x) polinom, melynek főegyütthatója (a _n ) egy^[3], és amelynek s darab különböző gyöke (x _i ) van, rendre α _i multiplicitással. Ekkor a polinom gyöktényezős alakja a következő:

A polinom deriváltja előállítható a következő alakban:

ahol egyszeres gyökök esetén α _i -1=0, így az adott gyök nem szerepel a szorzatban. Legyen g(x) (7.46) és (7.47) legnagyobb közös osztója, azaz:

Látható, hogy lesz a keresett csak egyszeres gyököket tartalmazó polinom. A g(x) polinomot az ún. Euklidesz-féle algoritmussal lehet előállítani. Állítsuk elő az f(x) és deriváltja hányadosát (q₀(x)) és az osztás maradékát (r₁(x))! Ekkor felírható a következő egyenlőség:

Folytassuk, és készítsünk egy sorozatot mindaddig, amíg a maradék nullává nem válik:

Bizonyítható, hogy r _k (x) lesz a legnagyobb közös osztó, azaz g(x). Most már csak p(x) meghatározása van hátra:

Az algoritmust megvalósító programrészlet a következő:

……………
'a az eredeti polinom együtthatói
Call Derivalt1(n, a, bb)'bb a derivált együtthatót tartalmazza
Call AegyenloB(n, aa, a) ' az eredeti tömb aa-ba
i = n ' a polinom fokszáma
j = n - 1 ' a derivált fokszáma
Do
Call PolDiv(i, j, aa, bb, cc)
Call PolFokszam(i, aa) ' a maradék fokszáma
' aa-t es bb-t fel kell cserélni, valamint a fokszámaikat is
Call AegyenloB(i, cc, aa)
Call AegyenloB(j, aa, bb)
Call AegyenloB(i, bb, cc)
k = i
i = j
j = k
LoopUntil (k = 0) And (cc(0) = 0) ' osztás, amíg a maradék 0nem lesz
Call AegyenloB(n, bb, a)      ' az eredeti polinom bb-be
Call PolDiv(n, i, bb, aa, cc) ' osztva a legnagyobb közös osztóval
k = n - i
Call PolFokszam(k, cc)
' cc-ben csak egyszeres gyökök vannak, fokszáma=k
……………

A következőkben igazolás nélkül közlünk egy olyan eljárást, amellyel intervallum határozható meg a gyökök abszolút értékére. Legyen adott a 7.37. alakú egyenlet. Tegyük fel, hogy a _n és a₀ nem zérus. Határozzuk meg a következő értékeket:

Bizonyítható, hogy az egyenlet minden gyökére igaz a következő egyenlőtlenség:

Az algoritmust megvalósító eljárás a következő:

Sub GyokInterv(n As Integer, ByRef a() As Double, _
ByRef al As Double, ByRef fel As Double)
Dim ak1 As Double, ak2 As Double
Dim i As Integer
ak1 = Abs(a(0))
ak2 = Abs(a(1))
For i = 1 To n - 1
If Abs(a(i)) > ak1 Then ak1 = Abs(a(i))
If Abs(a(i + 1)) > ak2 Then ak2 = Abs(a(i + 1))
Next i
fel = 1 + ak1 / Abs(a(n))
al = 1 / (1 + ak2 / Abs(a(0)))
End Sub

Az algoritmusnak a _n =1 főegyütthatójú, csak egyszeres gyököket tartalmazó polinomra van szüksége. Ilyet mi már elő tudunk állítani.

Az eljárás két részből áll. Az első részben egy polinomokból álló sorozatot (Sturm-lánc) kell meghatározni. Az algoritmus emlékeztet az Euklideszi algoritmusra, mert itt is polinomosztással kell dolgoznunk, csak a kiinduló két polinom lesz más, és az osztási maradék mínusz egyszeresére lesz szükségünk. Nézzük ezek után a számítás leírását, ahol p(x) a (7.51) szerint előállított k-ad fokú polinom:

Az így előállított p₀(x),…, p _k (x) polinomsorra lesz szükségünk. Ha egy [a,b] intervallumban szeretnénk megtudni a valós gyökök számát, akkor állítsuk elő a p₀(a),…,p _k (a) sorozatot, hagyjuk ki belőle az esetleges nullákat, majd számoljuk meg, hogy p₀(a)-tól p _k (a)-ig hány előjelváltás volt. Az így kapott értéket jelöljük v _a -val. Hasonlóan járjunk el b-vel is, az eredmény legyen v _b . Bizonyítható, hogy az [a,b]-ben lévő valós gyökök száma megegyezik a v _a -v _b különbséggel.

A következő két eljárás az algoritmus egy lehetséges megvalósítására mutat példát. Vegyük észre, hogy az előző pontban bemutatott intervallumkereső eljárás két intervallumot – egyet a pozitív számokra, egyet a negatívokra – adott. Ez a magyarázata, hogy az eljárások két darabszámot állítanak elő:

Function ElojelValtas(n As Integer, a() As Double) As Integer
Dim s As String, i As Integer
' Megszámolja, hogy az a-ban levő sorozatban hány előjelváltás van.
' A nulla nem számít.
s = ""
' a számokat az előjelükkel helyettesíti, a nullát kihagyja
For i = 0 To n
Select Case Sgn(a(i))
Case 1: s = s + "+"
Case -1: s = s + "-"
End Select
Next i
Do While InStr(s, "++") > 0 'az ismétlődések kiszűrése
i = InStr(s, "++")
s = DeleteStr(s, i, 1)
Loop
Do While InStr(s, "--")> 0 ' az ismétlődések kiszűrése
i = InStr(s, "--")
s = DeleteStr(s, i, 1)
Loop ' csak +-+- szerű elemek maradnak a sorozatban
ElojelValtas = Len(s) - 1
EndFunction

SubValosDb(n As Integer, ByRef a() As Double, _
ByRef negdb As Integer, ByRef pozdb As Integer)
Dim aa() As Double, bb() As Double, cc() As Double
Dim al As Double, fel As Double
Dim v1() As Double, v2() As Double, v3() As Double, v4() As Double
Dim i As Integer, j As Integer, k As Integer, L As Integer
Dim vege As Boolean
' Sturm-féle gyökszámlálás
' A két intervallum szélein vizsgál.
' A tömbökben a Sturm-sorozat helyettesítési értékei vannak.
ReDim v1(n)
ReDim v2(n)
ReDim v3(n)
ReDim v4(n)
Call GyokInterv(n, a, al, fel)
v1(0) = BehelyValos(fel, n, a)
v2(0) = BehelyValos(al, n, a)
v3(0) = BehelyValos(-fel, n, a)
v4(0) = BehelyValos(-al, n, a)
Call AegyenloB(n, aa, a)
i = n
L = 2
Call Derivalt1(n, a, bb)
j = n - 1
Call PolFokszam(j, bb)
v1(1) = BehelyValos(fel, j, bb) ' a derivált helyettesítési értékei
v2(1) = BehelyValos(al, j, bb)
v3(1) = BehelyValos(-fel, j, bb)
v4(1) = BehelyValos(-al, j, bb)
vege = False
Do
Call PolDiv(i, j, aa, bb, cc)
Call PolFokszam(i, aa)
Call NegalA(i, aa)
v1(L) = BehelyValos(fel, i, aa)
v2(L) = BehelyValos(al, i, aa)
v3(L) = BehelyValos(-fel, i, aa)
v4(L) = BehelyValos(-al, i, aa)
L = L + 1
If i = 0 Then vege = True
Call AegyenloB(i, cc, aa)
Call AegyenloB(j, aa, bb)
Call AegyenloB(i, bb, cc)
k = i
i = j
j = k
Loop Until vege
L = L - 1
pozdb = ElojelValtas(L, v2) - ElojelValtas(L, v1)
negdb = ElojelValtas(L, v3) - ElojelValtas(L, v4)
End Sub

Közelítő módszernek az érintőmódszert választottuk. Mint láttuk, az eljárás általában konvergens, de ha rossz helyről indulunk, akkor körbe-körbe járhat. A másik probléma, hogy kiléphet az intervallumból kevésbé zavaró, mert akkor is talál gyököt, csak másikat. Mivel mi minden gyököt szeretnénk megkeresni, és a sorrend nem számít, ezzel a „hibával” nem foglalkozunk. A körbe-körbe járás csak speciális esetben fordul elő, ezért a kezdeti értékét véletlenszám-generátor felhasználásával úgy állítottuk elő, hogy még a feltétel is teljesüljön. Ha ezer lépés után sem elég pontos a közelítés, akkor a program átvált az intervallumfelezésre.

A valós esetet megvalósító programrészlet a következő:

Call GyokInterv(k, cc, al, fel)
al = -al
fel = -fel
' első derivált cv, nv
Call Derivalt1(k, cc, cv)
nv = k - 1
Call PolFokszam(nv, cv)
' második derivált cvv, nvv
Call Derivalt1(nv, cv, cvv)
nvv = nv - 1
Call PolFokszam(nvv, cvv)
' megfelelő kezdeti érték keresése ->sign(f)=sign(f")
j = 0
Do
If negdb > 0 Then
x = (fel - al) * Rnd + al
Else
x = (-al + fel) * Rnd - fel
End If
j = j + 1
If j = 25 Then
al = al / 1.1
fel = 1.1 * fel
j = 0
End If
Loop UntilSgn(BehelyValos(x, k, cc)) = Sgn(BehelyValos(x, nvv, cvv))
' megoldás érintővel
x0 = x
j = 0
Do
x = x0 - BehelyValos(x0, k, cc) / BehelyValos(x0, nv, cv)
If (Abs(1-x/x0)<epsz) And (Abs(BehelyValos(x, k, cc)) < epsz) ThenExit Do
x0 = x
j = j + 1
If j = 1000 Then
x = IntervFelezes(al, fel, k, cc)
Exit Do
End If
Loop
Range("F" + Szam(gysz)).Value = x
If x < 0Then
negdb = negdb - 1
Else
pozdb = pozdb - 1
End If
gysz = gysz + 1
bb(1) = 1
bb(0) = -x
Call PolDiv(k, 1, cc, bb, aa) ' osztás a gyöktényezővel
k = k - 1
Call PolFokszam(k, aa)
Call AegyenloB(k, cc, aa)

Hasonlóan jártunk el komplex esetben is, de ott az első közelítő érték valós részét választottuk véletlenszerűen, a képzetes mindig egy lett:

Call GyokInterv(k, cc, al, fel)
' elsőderivált cv, nv
Call Derivalt1(k, cc, cv)
nv = k - 1
Call PolFokszam(nv, cv)
' második derivált cvv, nvv
Call Derivalt1(nv, cv, cvv)
nvv = nv - 1
Call PolFokszam(nvv, cvv)
' megfelelő kezdeti érték keresése ->sign(f)=sign(f')
Do
xr = (fel - al) * Rnd + al
Loop Until Sgn(BehelyValos(xr, k, cc)) = Sgn(BehelyValos(xr, nvv, cvv))
' megoldás érintővel
xr0 = xr
xi = 1
xi0 = 1
Do
Call BehelyKomplex0(xr0, xi0, k, cc, yr, yi)
Call BehelyKomplex0(xr0, xi0, nv, cv, yr0, yi0)
nev = yr0 * yr0 + yi0 * yi0
xr = xr - (yr * yr0 + yi * yi0) / nev
xi = xi - (yr0 * yi - yi0 * yr) / nev
If (Abs(1 - Sqr((xi * xi + xr * xr) / (xi0 * xi0 + xr0 * xr0))) < epsz) _
And (Sqr(BehelyKomplex(xr, xi, k, cc)) < epsz) Then ExitDo
xr0 = xr
xi0 = xi
Loop
If Abs(BehelyValos(xr, k, cc)) < epsz Then
Range("F" + Szam(gysz)).Value = xr
gysz = gysz + 1
n = n - 1
bb(1) = 1
bb(0) = -xr
Call PolDiv(k, 1, cc, bb, aa)
k = k - 1
Else
Range("F" + Szam(gysz)).Value = xr
Range("G" + Szam(gysz)).Value = xi
Range("F" + Szam(gysz + 1)).Value = xr
Range("G" + Szam(gysz + 1)).Value = -xi
gysz = gysz + 2
n = n - 2
bb(2) = 1
bb(1) = -2 * xr
bb(0) = xr * xr + xi * xi
Call PolDiv(k, 2, cc, bb, aa)
k = k - 2
End If
Call PolFokszam(k, aa)
Call AegyenloB(k, cc, aa)

Mint látható, mindkét programrészben szerepel a gyöktényezővel való osztás. Valós esetben ez magától értetődő, de vajon mi a helyzet a komplex gyökök esetén? Ha megtaláltunk egy komplex gyököt, akkor igazából kettőt találtunk, mert a konjugáltja is gyök. Legyen például a gyök. Ekkor behelyettesítve egy 7.37. egyenletbe, majd a valós és képzetes részeket szétválogatva:

Az egyenlőség nyilván csak akkor állhat fenn, ha mind a valós, mind pedig a képzetes rész zérus. Mivel gyököt helyettesítettünk be, ez nyilván így is van. Abban az esetben, ha a konjugáltat írjuk az egyenletbe, akkor csak annyi változik, hogy a képzetes rész előjele ellentétes lesz:

Azonban az előzőekben leírtak miatt ez is kielégíti az egyenletet, tehát a konjugált is gyök. Ezek után lehetőségünk van egyszerre két gyöktényezővel osztani. Ennél is jobb azonban a helyzet, mert:

Azaz egy valós együtthatójú másodfokú polinommal oszthatunk.

A klasszikus matematikai módszer nagyon jól használható a mátrixokkal kapcsolatos tételek bizonyítására, de gyakorlatilag használhatatlan számítógépes algoritmusként. Egy 12x12-es mátrix invertálása esetén a tárigény kevéssel nagyobb, mint double (8 bájt helyigényű) számábrázolás esetén. Több különböző numerikus módszer is ismert a mátrix invertálásra, melyek közül most egy olyat tárgyalunk, amely végeredményben nagyon hasonlít a Gauss-eliminációnál megismert algoritmusra.

Az eljárás alapgondolatát nézzük meg egy konkrét példán! Legyen az invertálandó mátrixunk a következő:

Szorozzuk meg -t a következő transzformációs mátrixszal:

Az eredmény:

Tulajdonképpen az első sor -4-szeresét hozzáadtuk a második sorhoz. A transzformációs mátrix valójában egy módosított egységmátrix, amely a Gauss-eliminációnál már ismertetett szorzót tartalmazza abban a pozícióban, ahová nullát szeretnénk kapni. Hasonlóan kaphatunk egy következő mátrixot is és így tovább, mindaddig, amíg a szorzat eredménye egységmátrix nem lesz:

A szorzat legyen egyenlő -vel. Azt látjuk, hogy , azaz az egységmátrix. Ez azonban csak úgy lehet, ha . Ha kezdetben felvettünk volna egy az -val azonos dimenziójú egységmátrixot, és sorba megszoroztuk volna a transzformációs mátrixokkal, akkor megkaptuk volna az inverz mátrixot. A szorzások azonban – mint láttuk – tulajdonképpen sorműveletek voltak, éppen úgy, mint a Gauss-eliminációnál. Összefoglalva, az algoritmus főbb lépései az alábbiakban láthatók:

Az algoritmust megvalósító Basic nyelvű függvények a következők:

Function Csere(n As Integer, i As Integer, ByRef a() As Double, _
ByRefinv() As Double) As Integer
' sorokat cserél mindkét mátrixban
' n a mátrix dimenziója
' az i. sor alatt kell keresni nem nulla értékeket
' a ebből lesz egység mátrix
' inv ez volt az egységmátrix
Dim j As Integer, k As Integer, cs As Double
Csere = 1 ' nem sikerült
For j = i + 1 To n
If a(j, i) <> 0 Then Exit For
Next j
If a(j, i) = 0 Then Exit Function
For k = 1 To n
cs = a(i, k)
a(i, k) = a(j, k)
a(j, k) = cs
cs = inv(i, k)
inv(i, k) = inv(j, k)
inv(j, k) = cs
Next k
Csere = 0 ' sikerült
End Function

Function MInverz(n As Integer,a() As Double,ByRef inv() As Double) As Integer
' a mátrix dimenziója
' a az invertálandó mátrix
' inv itt lesz az inverz mátrix
Dim i As Integer, j As Integer, k As Integer
Dim sz As Double
ReDim inv(n, n)
For i = 1 To n ' kezdetben egységmátrix
inv(i, i) = 1
Next i
For i = 1 To n - 1 ' alsó háromszög kinullázása
For j = i + 1 To n
If a(i, i) = 0 Then err = Csere(n, i, a, inv)
If err <> 0 Then Exit Function
sz = -a(j, i) / a(i, i)
For k = 1 To n
  a(j, k) = a(j, k) + sz * a(i, k)
  inv(j, k) = inv(j, k) + sz * inv(i, k)
Next k
Next j
Next i
If a(n, n) = 0 Then Exit Function
For i = n To 2 Step -1 ' felső háromszög kinullázása
For j = i - 1 To 1 Step -1
sz = -a(j, i) / a(i, i)
For k = 1 To n
  a(j, k) = a(j, k) + sz * a(i, k)
  inv(j, k) = inv(j, k) + sz * inv(i, k)
Next k
Next j
Next i
For i = 1 To n ' ettől lenne az egységmátrix
sz = a(i, i)
For j = 1 To n
inv(i, j) = inv(i, j) / sz
Next j
Next i
MInverz = 0
EndFunction

Abban az esetben, ha a mátrix rosszul kondicionált, ugyanazok a lehetőségek állnak rendelkezésünkre, mint a lineáris egyenletrendszer megoldásánál. Az eljárás közben – ha szükségünk van rá – előállíthatjuk a mátrix determinánsának az értékét is. Ezt akkor tehetjük meg, amikor az eredeti mátrixunk főátlója alatt már minden nulla. Ekkor a főátló elemeit összeszorozva megkapjuk a determináns értékét.

A gyakorlatban sokszor előfordul, hogy egy feladat megoldása során felírunk minden egyenletet, amit csak tudunk. Ekkor könnyen előfordulhat, hogy néhány egyenlet összefügg a többivel, tehát azokból előállítható. Szükségünk van tehát egy olyan módszerre, amellyel a „rengeteg” egyenletből egy olyan egyenletrendszert kapunk, amelyben már csak független egyenletek szerepelnek.

Mit is jelent az, hogy egy egyenlet független a többitől? Azt, hogy nem lehet előállítani a többiek lineáris kombinációjaként. Vegyük észre a hasonlóságot a vektorok függetlenségével! Legyenek az egyenletek együtthatói egy n elemű vektor koordinátái, és jelöljük őket v ₁, v ₂,… v _k -val. Ha most például n=k és a v _i vektorokkal minden n dimenziós vektor előállítható, akkor v _i -k függetlenek és teljes bázist alkotnak. Ha azonban van olyan v _j a vektorok között, amely előállítható a következőképpen:

akkor a vektorok összefüggőek, és nem alkotnak bázist. A bázisfaktorizáció célja az, hogy a v ₁, v ₂,… v _k vektorokból olyan w ₁, w ₂,… w _l vektorokat állítson elő, amelyek függetlenek, így bázist alkotnak (általánosan ).

Nézzünk egy konkrét példát! Legyen adva a háromdimenziós térben három vektor: . Ha ezek egy síkba esnek, akkor nem feszítik ki a háromdimenziós teret, azaz nem lehet tetszőleges 3-dimenziós vektort előállítani a segítségükkel. Annak eldöntése, hogy függetlenek-e ránézésre nem mindig egyszerű. Erre szolgál a bázisfaktorizáció.

Az eljárást az előző példán mutatjuk be. Hozzuk létre a sorvektorokként kezelt v ₁, v ₂, v ₃ vektorokból az mátrixot! Végezzük el a felső háromszög mátrixszá alakítást a Gauss‑eliminációnál ismertetett módon:

Látható, hogy a harmadik vektor zérussá vált, így a bázisvektorok:

Ez a módszer általánosságban is így működik. Eredményként egy trapéz (szerencsés eseten háromszög) alakú elrendezéshez jutunk, amely már megkönnyíti az egyenletrendszer megoldását (7.78).

Természetesen előfordulhat, hogy kevesebb az egyenletünk, mint az ismeretlenek száma. Ekkor vagy fel tudunk még írni kellő számú független egyenletet, vagy marad néhány szabad paraméterünk. Az irodalomban erre a feladatra is – nevezetesen kevesebb az egyenletünk, mint az ismeretlen – sok módszer található, például az ún. projektor mátrixokra építő Purcell-módszer [8.] .

Nemlineáris egyenletrendszereket alapvetően kétféleképpen szokás megadni:

Bármelyiket is ismerjük, a másik egyszerűen előállítható a következő helyettesítéssel:

Az iterációs sorozatot (7.84) alakból állíthatjuk elő a szokásos módon:

feltéve, hogy az x _i vektorok benne maradnak ɸ értelmezési tartományában. A konvergencia feltétele, hogy:

Ha egy kicsit figyelmesebben megnézzük ezt a feltételt, akkor észrevehetünk egy differenciahányados szerű kifejezést is benne. Abban az esetben, ha és folytonosak, az értelmezési tartományon a feltétel így írható fel:

ahol

az ún. Jacobi mátrix. Ez már nagymértékben hasonlít az egyváltozós iterációnál megismert konvergencia feltételre. A feltétel eléggé szigorú, azaz sokszor nem elégíti ki az egyenletrendszerünk. Az egyváltozós esetben bemutatott módszer azonban – természetesen átalakítva – itt is alkalmazható. Induljunk ki a (7.83) alakból! Ekkor a részletek mellőzésével, ha x ₀ az első közelítés:

és

Bár kicsi az esélye, de ha nem állítható elő, mert a mátrix nem invertálható, akkor próbáljunk másik kezdeti értéket választani.

Példaként nézzünk egy feladatot!. Legyen adott a (7.15. ábra) ábrán látható mechanizmus. Szeretnénk megtudni, mi lesz a C csukló pozíciója, amikor B csukló vízszintesen az x=-1, y=0 pontban van. Az adatok: egység.

A megoldáshoz két kör metszéspontját kell meghatároznunk, mivel ekkor a C pont egy (-1,0) középpontú, 5 egység sugarú, és egy D középpontú, 2 egység sugarú kör közös pontjában lesz. Az egyenletrendszer tehát a következő:

Ennek az egyenletrendszernek viszonylag könnyű előállítani a pontos megoldását, éppen ezért lesz kiváló arra, hogy a módszereket ellenőrizni tudjuk vele. A levezetés mellőzésével a megoldás:

Nézzük most a klasszikus iterációt! Ekkor az egyenletrendszer alakja:

Kezdeti értéknek válasszuk az x=1 és y=2 kiindulási adatokat. Táblázatkezelőbe behelyettesítve a (Táblázat 7.16) táblázatban látható adatokat kapjuk.

Próbálkozzunk a javított változattal is! Ebben az esetben az

alakú egyenletrendszerre van szükségünk. A derivált mátrix pedig az Jacobi mátrix, mely most a következő:

Behelyettesítve a kezdőértéket és invertálva a mátrixot:

Most már felírhatjuk az iterációs összefüggést:

Táblázatkezelő alkalmazásával a (Táblázat 7.17) táblázatban látható eredményeket kapjuk.

Az iteráció viszonylag gyors, mert tíz lépés után leállíthattuk. Meg kell azonban jegyezni, hogy nem mindegy, hogy honnan indítjuk a számítást, mert lehet divergens is a sorozat.

Az érintőmódszer általánosításához az előző fejezetben bemutatott, javított iterációs módszert kell csak aktualizálni. Az mátrix meghatározása során mindig az aktuális értékeket helyettesítsük be, azaz:

Nézzük most meg ennek hatását az előző feladat esetében! A számításokat a Táblázat 7.18 tartalmazza.

A konvergencia itt – hasonlóan az egyváltozós változathoz – gyors. Annak vizsgálata azonban, hogy az adott esetben a számítások eredményre fognak-e vezetni, azaz az adott paraméterekkel az iteráció konvergens-e, eléggé bonyolult, időigényes, éppen ezért nem is szokták elvégezni. Ehelyett azt figyelik, hogy nem tartanak-e az értékek (elég, ha csak egyikük) végtelenhez.

Mint a 7.98. összefüggésből látható, minden iterációs lépésnél el kell végeznünk egy mátrix invertálását. Ennél valamivel kevesebb számítást igényel, ha az egyenletet megszorozzuk a Jacobi mátrixszal, mert akkor csak egy lineáris egyenletrendszer kell megoldanunk minden lépés során: