Numerikus gyökvonás

Kollár Zsolt
4.évfolyamos villamosmérnök szakos hallgató
Neptun kód: CTQWYF

1.  Bevezetés

A gyökvonás alapvető műveletként van manapság kezelve .  Ezt a műveletet a mérnöki gyakorlatban sűrűn alkalmazzák:nem csak a Pithagorsz tétel esetén,hanem sok más fizikai számításkor is.
Az f : x->(Sqrt[x])függvény az f:x->(x^2) függvény inverze, de csak az xϵ[0, ∞] intervallumon van értelmezve.

\(Plot[{x\^2}, {x, 0, 2}, \ Ticks -> {{0, 1, 2}, {0, 1}}, 
      PlotStyle -> {{Hue[0.99]}}, 
      PlotLabel -> 
        StyleForm[\*""Az \!\(x\^2\) függvény"", 
          FontSize -> 12]]\)
Plot[{Sqrt[x], x}, {x, 0, 2}, 
        Ticks -> {{0, 1, 2}, {0, 1}}, 
        PlotStyle -> {{Hue[0.19]}, {Hue[0.99]}}, 
        PlotLabel -> 
          StyleForm["Az sqrt(x) és az y=x \
függvény", FontSize -> 12]]

⁃  Graphics  ⁃

Null (⁃  Graphics  ⁃)
Null

Jól látható, hogy a (Sqrt[x]) függvény a [0, 
          1] intervallumon mindig nagyobb értéket vesz fel, mint az a szám, amiből gyököt vonunk .  Az [1, ∞] intervallumon pedig a gyök mindig kisebb,mint maga a szám .  Negatív számokra pedig nincsenen értelmezve a függvény.
Azonban az "sqrt" függvény numerikus kiszámítása korántsem ilyen egyszerű és triviális .  Fontos szempont függvényértékek numerikus kiszámítása esetén,hogy visszavezessük egy olyan kifejezésre,amely csak alapműveleteket tartalmaz (összeadás,  kivonás,  
          szorzás,  osztás), mivel ezek már a számítógép számára könnyen és gyorsan elvégezhetőek.

A továbbiakban a numerikus gyökvonás intervallumát leszűkítjük a [0, 1] intervallumra .  Mivel a DSP-k 
(Digital Signal Processor) törtszámok ábrázolása és aritmetikája miatt  ez jól megvalósítható .  Egynél nagyobb számok esetén egy előskálázást alkalmazunk,ami beskálázza az adott számot a [0, 1]intervallumba .  A gyökvonás után pedig a skálázási faktor
gyökével visszaskálázzuk a kapott értéket.

A skálázási faktor gyökét kiszámolni már triviális? Csak néhány, rögzített faktort használunk?

\(Ebben\ a\ dolgozatban\ két\ módszer\ lesz\ bővebben\ ismeretetve\ és\ \(tárgyalva : \ 
          A\ Newton\ módszer\ és\ a\ gyökfüggvény\ sorba\ fejtése . 

A\ Newton\ módszer\ egy\ rekurzív\ algoritmus\ szerint\ közelíti\ a\ keresett\ szám\ gyökét\)\), \(ahol : 
		x\_\(k + 1\) = \(1\/2\) \((x\_k + 
              a\/x\_k)\)\), 
\(A\ képletben\ szereplő\ \(változók : \ 
        x\_k\ az\ "a"\ szám\ feltételezett\ gyöke\), \ 
    x\_\(k + 1\)\ a\ következők\ közelítés . 
Megfigyelhető, \ 
    hogy\ ha\ x\_k = \@\(\(a\)\(\ \)\)\ választással\ élünk\ , \ 
    vagyis\ x\_k\ pontosan\ az\ "a"\ szám\ gyöke, \ 
akkor\ a\ kifejezés\ \ \@a\ értéket\ veszi\ fel . \ 
        Az\ algoritmus\ szempontjából\ fontos\ \(kérdés : \ 
        Hogyan\ válasszuk\ meg\ x\_0\ értéket?\ 
            A\ Newton\ módszer\)\)
előnye DSP esetén , hogy a skálázás kettes számrendszerben, vagyis a 2 - vel való szorzás, könnyen elvégezhető: egyszerű eltolás.
(analógia :  
              mint a 10 - es számrendszerben a 10 - 
                zel való szorzás) .  
          
          Problémát okoz az osztás művelete,  bár szerencsére ez 
visszavezethető a kivonás műveletére.

Ha muszáj vmit mondani, akkor a [0,1]-beliszám gyökét ezzel közelíteném: (1+x)/2.

\(A\ gyökfüggvény\ sorba\ fejtése\ a\ [0, 1]\ \(intervallumon : 
	\@x\) = 
          1 - \(1\/2\) \((1 - 
                    x)\) - \(\(1*1\)\/\(2*4\)\) \((1 - x)\)\^2 - \(\(1*1*3\)\/\(2*4*6\)\) \((1 - x)\)\^3 -  ... \)
Ezen módszer implementálása esetén fontos kérdés,hogy hány tagot veszünk figyelembe .  DSP-n történő implementálásának előnye, hogy a processzorok támogatják a polinomok kiszámítását!

Series[(Sqrt[x]), {x, 1, 9}] // No
rmal

A Horner szabály alkalmazásával egy polinom gyorsan kiszámítható a minimális számú szorzás és összeadás felhasználásával:
Pl .  
          f (x) = d + cx + (bx^2) + 
          (ax^3) = ((((ax) + b) x) + c) 
            x + d.

<< Algebra`Horner`

Horner[d + c x + b (x^2) + a (x^3), x]

2.  A Newton - 
        módszer konvergálása a [0, 
            1] intervallumon

A fent említett módszer segítségével keressük meg az a=0.25 gyökét,(Subscript[x, 0])=1  kezdőértéket választva és vizsgáljuk meg,  
              hogy milyen gyorsan konvergál a tényleges gyökhöz (jelen esetben a 0.5-höz).

ClearAll[f]; f[a_, x_] := (x + a/x)/2.0

NestList[ff, x0, 4]

NestList[f[0.25, #] &, 1, 6]

NestList[(# + 0.25/#)/2 &, 1, 6]

x0 = {1};
a = 0.25;
For[i = 0, i < 6, i++; 
        x0 = Join[x0, {(Extract[x0, i] + (a/Extract)[x0, i])/2}]; 
        Print[i, ". ", "x=", Extract[N[x0], i]]]
ListPlot[Table[Extract[x0, k], {k, 1, 7}], 
        PlotStyle -> PointSize[0.02], PlotRange -> {0, 1}]

Jól látható, hogy az algoritmus exponenciálisan közelíti a szám gyökét .  Jelen esetben már az 5. lépésben elérte a tényleges gyököt .  Hogyan viselkedik az algoritmus akkor, ha a kezdőérték a tényleg gyöknél kisebbnek választjuk?
\(A\ következő\ példában\ továbbra\ is\ a\ a\)=0.25 gyökét keressük, de mostkezdőértéknek a kiválasztott szám gyökénél egy kisebb számot 
pl .  (Subscript[x, 0]) = 
          0.25 - t választunk.

x0 = { .25};
a = 0.25;
For[i = 0, i < 6, i++; 
        x0 = Join[x0, {(Extract[x0, i] + (a/Extract)[x0, i])/2}]; 
        Print[i, ". ", "x=", Extract[N[x0], i]]]
ListPlot[Table[Extract[x0, k], {k, 1, 7}], 
        PlotStyle -> PointSize[0.02], PlotRange -> {0, 1}]

ListPlot[NestList[(# + 0.25/#)/2 &, 0.25, 6], 
      PlotStyle -> {PointSize[0.02], Red}]

Jól látható, \(hogy\ az\ első\ rekurzíós\ lépés\ után\ a\ módszer\ "felugrik"\ a\ tényleges\ gyök\ fölé\), \(és\ \ onnan\ konvergál\ exponenciálisan\ a\ tényleges\ \ gyökhöz\);jelen esetben ugyanazokkal a részeredményekkel, \(mintha\ x\_0 - t\ \ 1 - 
          nek\ választottuk\ volna . \ \ A\ Newton\ módszernél\ általánosságban\ mondhatjuk\), \(hogy\ a\ rekurzíó\ azzal\ gyorsítható\), ha  (Subscript[x, 0])-t (Sqrt[a])-\(hoz\ \ minél\ közelebb\ eső\ értéknek\ választjuk . \ 
              Persze\ kérdés\),hogy hogyan?

3.  
        A gyökfüggvény sorbafejtésénak konvergálása

A sorbafejtés implementálásának fő kérdése, hogy hány tagot érdemes figyelembe venni .  Ebben a fejezetben ezt fogjuk vizsgálni .  Első közelítésben nézzük egy egyszerű esetet: vegyük figyelembe a sorbafejtés első 5 tagját. 

Vizsgáljunk 5 esetet: a számok amiből gyököt akarunk vonni legyenek: 0.1, 0.25, 0.5, 0.81,1.

ai={1/2,  1/8,  
              3/48,  (3*5)/(2*4*6*8),  (3*5*7)/(48*8*10)};

x = 0.1;
xn = {1};
For[i = 0, i < 5, i++; 
        xn = Join[
            xn, {(Extract[xn, i] - Extract[ai, i]*(1 - x)^i)}]]
v1 = Table[{i, Extract[xn, i]}, {i, 1, 5}]

x = 0.25;
xn = {1};
For[i = 0, i < 5, i++; 
        xn = Join[
            xn, {(Extract[xn, i] - Extract[ai, i]*(1 - x)^i)}]]
v2 = Table[{i, Extract[xn, i]}, {i, 1, 5}]

x = 0.5;
xn = {1};
For[i = 0, i < 5, i++; 
        xn = Join[
            xn, {(Extract[xn, i] - Extract[ai, i]*(1 - x)^i)}]]
v3 = Table[{i, Extract[xn, i]}, {i, 1, 5}]

x = 0.81;
xn = {1};
For[i = 0, i < 5, i++; 
        xn = Join[
            xn, {(Extract[xn, i] - Extract[ai, i]*(1 - x)^i)}]]
v4 = Table[{i, Extract[xn, i]}, {i, 1, 5}]

x = 1;
xn = {1};
For[i = 0, i < 5, i++; 
        xn = Join[
            xn, {(Extract[xn, i] - Extract[ai, i]*(1 - x)^i)}]]
v5 = Table[{i, Extract[xn, i]}, {i, 1, 5}]

v = 
        Join[v1,  v2, v3, v4, v5];
ListPlot[v, PlotStyle -> PointSize[0.02], 
        PlotRange -> {0, 1.5}]
Null

\(A\ konvergencia\ az\ 1\ körüli\ számok\ esetén\ a\ legjobb . \ A\ polinom\ 0 - 
              hoz\ közeli\ értékekre\ ad\ rossz\ közelítést\), mivel ilyenkor a polinom magasabb hatványú tagjai is még számottevőek, vagyis ezeket is figyelembe kéne venni .  Minél közelebb van a szám az  1-hez,annál jobb a közelítés .  x=1 gyöke már az első lépésben jó, az x=0.81 gyöke (0.9) már a második lépésben nagyon jó közelítés .  
              Megfigyelhető, hogy az x=0.1 gyöke (0.31622) még az 5. tag figyelembevételével is eltér a második 
tizedesjegyben . 
Ezen polinomszámítási eljárás javítható, ha a 0 és 0.5 között számokat felskálázzuk a 0.5 és 1 közöttitartományba. Természetesen a kapott értéket majd utólag vissza kell skálázni a skálázási faktor gyökével.

Normal[Series[(Sqrt[x]), {x, 1, 6}]] /. x ->  .1

Thread[x -> {0.1, 0.25, 0.5, 0.81, 1.0}]

(Normal[Series[(Sqrt[x]), {x, 1, 6}]] /. # &) /@ 
      Thread[x -> {0.1, 0.25, 0.5, 0.81, 1.0}]

4.  A két módszer megvalósítása processzorral

4.1. DSP Aritmetikai egysége, számábrázolás

A DSP (Digital Signal Processzor) számábrázolása bináris .  Típustól, gyártmánytól függően  véges helyiértéken képes csak ábrázolni .  M  anapság  8,16 ,\(\((32)\) 
              bites\ processzorokat\ \ alkalamaznak . \ 
                Ez\ azt\ jelenti\ hogy\ \(a : \ 

\ \ \ 8\ bites\)\) = 2^8=0-tól 256-ig képes a számokat kezelni.
	16 bites = 2^16=0-tól 65536-ig képes a számokat kezelni.

Az adott DSP-től függen természetesen képes előjeles (signed)számokat is kezelni,ebben az esetben a legmagasabb helyiértékű bitet 
\(előjelnek\ van\ fentartva, \ így : 

\ \ \ \ 8\ bites\  = \ \(\(+\(\(/\)\(-\ \ 2^7\)\)\) = \(-128\) - tól\ 128 - 
            ig\ képes\ a\ számokat\ kezelni\)\)
       16 bites  =  + /- 2^15 = -32768 - 
              tól 32768 - 
              ig képes a számokat kezelni

    \(Vagyis\ a\ számábrázolási\ tartomány\ nem\ csökkent\ csupán\ el\ lett\ tolva\ a\ tartomány\ felével\ negatív\ irányba . 
\ 
              Természetesen\ néhány\ processzor\ képes\ törtszámüzemmódban\ \ is\ számolni\), \(ilyenkor\ a\ számábrázolási\ tartomány\ a\ [\(-1\), 1]\), 
amit feloszt  annyi részre,\(ahány\ számot\ képes\ kezelni . \ 
              Az\ ilyen\ típusú\ processzorokra\ próbáljuk\ a\ gyökvonást\ optimalizálni . 
\(Példul : 
                8\ bites\ processzor\ \(esetén : \ \ 
0100\ 0000\_2\ Jelöli\ a\  + 0.5\_10 - t\)\)\), vagy 11  1  0 (Subscript[0000, 2]) Jelöli a -(Subscript[0.75, 10])-t.  
Ez azt jelenti  tulajdonképpen az előjelbit után van a bináris pont .  A következőkben azt fogjuk vizsgálni, \(hogy\ az\ előzőekben\ ismeretet\ Newton\ módszer\ illetve\ a\ gyökvonás\ sorbafejtése\ hogyan\ működik\),\(ha\ a\ számábrázolásunk\ csak\ véges\ egész\), \(illetve\ melyik\ módszer\ implementálható\ \(\(egyszerűbben\)\(.\)\)\)

4.2. A Newton módszer viselkedése és implementálása

Az osztás művelete processzorok esetén vissza van vezetve a kivonás műveletére .   Az algoritmus programozásának  2 féle változatának vázlatos lépései:

I. Változat:
 1. Lépés: kiszámoljuk (a/2) értékét
\(2. \ \(Lépés : \ 
          kiszámoljuk\ a\/\(2  x\_k\)\ értékét\)\)
3.  Lépés :  
            kiszámoljuk (a/(2  (Subscript[x, k]))) + (1/2) 
                (Subscript[x, k]) kifejezést     

II. Változat:
1. Lépés: kiszámoljuk (a/(Subscript[x, k])) értékét
2.  Lépés :  
          kiszámoljuk (a/(Subscript[x, k])) + 
            (Subscript[x, k]) értékét
3.  Lépés :  
            kiszámoljuk 
               (1/2) ((a/(Subscript[x, k])) + (Subscript[x, k])) kifejezést

\(Digitális\ világban\ a\ kettő\ hatványaival\ történő\ osztás\ \((szorzás)\)\ művelete\ könnyen\ és\ gyorsan\ elvégezhető, \ hiszen\ csak\ jobbra\ történő\ \ eltolást\ jelent . \ \ 
\((Analógia : \ 
                mintha\ a\ tizes\ számrendszerben\ 10 - 
                  el\ \(\(osztanánk\)\(.\)\))\)\ Mivel\ nincs\ törtszámábrázolásunk\ ezért\ a\ szimuláció\ során\ , \ 
        osztás\ műveletet\ után\) 
\(mindig\ lefele\ kerekítünk, \ 
      azaz\ a\ Floor[]\ függvényt\ használjuk . 
          Ezzel\ próbáljuk\ szimumálni\ a\ véges\ számábrázolás\ okozta\ \(problémát!\)\)

(Subscript[x, 0]) jó megválasztása nagyon fontos, mivel a DSP-ben osztás művelte történik, amely az eredmény törtrészét elveszítjük .  Ezért érdemes 
(Subscript[x, 0])= 1-t választani .  Ezen választás hátránya, hogy az algoritmus kiszámításakor, a második esetben tulcsordulás léphet fel a 2. lépésben, \(hiszen\ az\ összeg\ túllépheti\ az\ 1 - et . \ 
              A\ véges\ \ számábrázolást\ ebben\ az\ esetben\ nem\ törtszámokra\),hanem egész számokra fogjuk bemutatni, vagyis most
\(eltérünk\ az\ előzőekben\ előírt\ feltételtől\), \(hogy\ csak\ a\ [0, 
              1]\ intervallumon\ számolunk . 
Első\ lépésben\ teszteljük\ mindkét\ módszer\ egy\ ismert\ pl . \ a\)=25 értékkel.

I. Módszer eredményei

x0 = {1};
a = 25;
For[i = 0, i < 9, i++; 
        x0 = Join[
            x0, {(Floor[Extract[x0, i]/2] + 
                  Floor[Floor[a/2]/Extract[x0, i]])}]; 
        Print[i, ". ", "x=", N[Extract[x0, i]]]]
ListPlot[Table[Extract[x0, k], {k, 1, 8}], 
        PlotStyle -> PointSize[0.02], PlotRange -> {0, 13}]

Jól látható, hogy az algoritmus 4 és 5 között oszcillál. Ezen annak a következménye, hogy a 2-vel való osztásnál páratlan számok esetén 0.5-öt tévedünk. Ezért előbb érdemes
összeadni csak az után elosztani kettővel az összeget.
    
    II. Módszer eredményei

x0 = {1};
a = 25;
For[i = 0, i < 8, i++; 
      x0 = Join[
          x0, {(Floor[(Extract[x0, i] + Floor[a/Extract[x0, i]])/
                  2])}]; 
      Print[i, ". ", "x=", N[Extract[x0, i]]]]
ListPlot[Table[Extract[x0, k], {k, 1, 8}], 
      PlotStyle -> PointSize[0.02], 
      PlotRange -> {0, 15}]

A második módszer esetén az algoritmus a helyes értékhez konvergál. Azonban az  algoritmus ezen változatában is oszcilláció léphet fel, ha pl. a=24-et választunk.

x0 = {1};
a = 24;
For[i = 0, i < 8, i++; 
        x0 = Join[
            x0, {(Floor[(Extract[x0, i] + Floor[a/Extract[x0, i]])/
                    2])}]; 
        Print[i, ". ", "x=", N[Extract[x0, i]]]]
ListPlot[Table[Extract[x0, k], {k, 1, 8}], 
        PlotStyle -> PointSize[0.02], PlotRange -> {0, 15}]

Vizsgáljuk meg végül a  hibafüggvényt. Kiszámoljuk 1 és 100 közötti egész számok gyökeit az  algoritmus segítségével és minden esetben kivonjuk ezt  tényleges gyökből. 
 Jól látható, hogy a hiba a legkisebb helyiértéken belül marad, jelen esetben 0 és 1 között, valamint az is jól látható, hogy a négyzetszámoknál a hiba értéke 0.

dif = {}
x0 = 1;
For[i = 1, i < 101, i++, 
      x0 = 1; 
For[j = 0, j < 12, j++; 
        x0 = Floor[(Floor[i/x0] + x0)/2]]; 
      dif = Join[dif, {Sqrt[i] - x0}]]
ListPlot[Table[Extract[dif, k], {k, 1, 100}], 
      PlotRange -> { 0, 1}, 
      PlotLabel -> 
        StyleForm["A hibafüggvény", 
          FontSize -> 14]]

4.3. A gyökfüggvény sorba fejtésének viselkedés és implementálása

A gyökfüggvény sorbafejtésének előnye, hogy sokkal gyorsabban kiértékelhető, mint a Newton-féle rekurzív módszer. A Newton módszer leállási feltétele nem egyértelmű, hiszen,
mint már korábban láttuk, oszcilláció léphet fel. Ha (Subscript[x, 0])-t rosszul választjuk meg és távol van a kérdéses szám valós gyökétől, akkor Newton módszer csak sok lépésben fog konvergálni.
A korábban ismeretett gyökfüggvény első 3,4 vagy 5 tagját használva, illetve a korábban említett előskálázást alkalmazva már igen jó közelítést tudunk mondani az "a" szám gyökére. 
f(x)=a+bx+c(x^2) + (dx^3) alakban próbáljuk a kiértékelni, akkor 3 összeadást és 6 szorzást kell elvégezni. Ellenben a Horner alakban: f(x)=((dx+c) x+b)x+a, ami megegyezik az előző alakkal, már csak 3 szorzást és 3 összeadást kell elvégezni, ezzel az algoritmus sebessége gyorsítható. Néhány processzor rendelkezik kifejezetten polinomok kiértékelését támogató utasításokkal is.

5.  Összegzés

Mindkét módszernek megvan a maga előnye .   Jól látható, hogy a Newton módszer kiszámítása sok időt igényel, a gyökfüggvény sorbafejtés gyorsan kiértékelhető ám kevésbé pontos.  Igazán hatékony és pontos algoritmust a kettő előnyeinek ötvözésével tudunk elkészíteni: Először  a gyökfüggvény sorbefejtését értékeljük ki az első 3-4-5 tagra,majd a kapott értéket behelyettesítve  (Subscript[x, 0])-ba, a Newton módszert alkalmazva néhány (2-3) rekurziós lépésre,már igen pontos értéket kaphatunk a kérdéses szám gyökére .  
              Problémát okoz a 0 gyökének kiszámítása, mivel skálázni nem lehet, ezért külön esetként kell kezelni:Ebben az esetben nem fut le az algoritmus,hanem rögtön megoldásként visszaadja a helyes értéket, vagyis a nullát.

6.  Felhasznált irodalom

1,Texas Instruments :  
            TMS320C54x DSP  -  CPU and Peripherals  -  
              Reference Set Volume 2, 2001.
2,Texas Instruments :  
            TMS320C54x DSP  -  Mnemonic Instruction Set  -  
              Reference Set Volume 1, 2001.
3,Dr .  Arató Péter :  
            Logikai rendszerek tervezése, Műegyetemi kiadó, 1992.
4,Hack Frigyes :  
            Függvénytáblázatok, Matematikai összefüggések, Tankönyvkiadó, Budapest, 1969.
5,Szili  László  és  Tóth János :  Matematika és Mathematica, ELTE Eötvös Kiadó, Budapest,1996.
6,Werner, Burkhard: First Steps in Mathematica, Springer Verlag London Limited, 1994.
7,Bronstejn,Szemengyajev: Matematikai zsebkönyv, Műszaki kiadó, Budapest,1980.