Sekos riba

3. Monotoninės sekos

1. Monotoninių sekų apibrėžimas.

Apibrėžimas. Seka $\{x_{n}\}$ vadinama nemažėjančia (nedidėjančia), jei kiekvienas tos sekos elementas ne mažesnis (ne didesnis) už pirmesnįjį elementą, t. y., jei su visais numeriais n teisinga nelygybė

$x_{n}\leq x_{n+1}\;\;$ ($x_{n}\geq x_{n+1}$).

Nemažėjančios ir nedidėjančios sekos vadinamos monotoninėmis sekomis. Jei monotoninės sekos $\{x_{n}\}$ elementai su bet kuriuo numeriu n tenkina nelygybę $x_{n}<x_{n+1}\;\;$ ($x_{n}>x_{n+1}$), tai seka $\{x_{n}\}$ vadinama didėjančia (mažėjančia). Didėjančios ir mažėjančios sekos dar vadinamos griežtai monotoninėmis sekomis.

Monotoninės sekos yra aprėžtos arba iš viršaus, arba iš apačios; būtent, nedidėjančios sekos aprėžtos iš viršaus, o nemažėjančios aprėžtos iš apačios pirmaisiais savo elementais. Todėl nedidėjanti seka yra aprėžta iš abiejų pusių, kai ji aprėžta iš apačios, o nemažėjanti seka yra aprėžta iš abiejų pusių, kai ji aprėžta iš viršaus.

Pateiksime tris monotoninių sekų pavyzdžius.

1. Seka $1,\;1,\;{\frac {1}{2}},\;{\frac {1}{2}},\;...,\;{\frac {1}{n}},\;{\frac {1}{n}},\;...\;$ nedidėjanti. Ji aprėžta iš viršaus savo pirmuoju elementu, lygiu vienetui, o iš apačios - skaičiumi 0.

2. Seka $1,\;1,\;2,\;2,\;...,\;n,\;n,\;...\;$ nemažėjanti. Ji aprėžta iš apačios savo pirmuoju elementu, o iš viršaus neaprėžta.

3. Seka ${\frac {1}{2}},\;{\frac {2}{3}},\;{\frac {3}{4}},\;...,\;{\frac {n}{n+1}},\;...\;$ didėjanti. Ji aprėžta iš abiejų pusių: iš apačios pirmuoju elementu ${\frac {1}{2}},$ iš viršaus, pavyzdžiui, skaičiumi 1.

2. Monotoninės sekos konvergavimo požymis.

Įrodysime pagrindinę teoremą.

3.15 teorema. Jei nemažėjanti (nedidėjanti) seka $\{x_{n}\}$ aprėžta iš viršaus (iš apačios), tai ji konverguoja.

Pagal praeitą skirsnį seka $\{x_{n}\},$ atitinkanti 3.15 teoremos sąlygas, yra aprėžta. Todėl 3.15 teoremą galima trumpai formuluoti taip: jei monotoninė seka $\{x_{n}\}$ yra aprėžta iš abiejų pusių, tai ji konverguoja.

Įrodymas. Kadangi seka $\{x_{n}\}$ yra aprėžta, tai jos elementų aibė turi tikslųjį viršutinį rėžį ${\overline {x}}$ ir tikslųjį apatinį rėžį ${\underline {x}}.$ Įrodysime štai ką: jei $\{x_{n}\}$ – nemažėjanti seka, tai jos riba yra minėtasis tikslusis viršutinis rėžis; jei $\{x_{n}\}$ – nedidėjanti seka, tai jos riba yra minėtasis tikslusis apatinis rėžis ${\underline {x}}.$ Išnagrinėsime tik nemažėjančią seką, nes samprotavimai apie nedidėjančią seką būtų analogiški.

Kadangi ${\overline {x}}$ yra sekos $\{x_{n}\}$ elementų aibės tikslusis viršutinis rėžis, tai kiekvieną $\varepsilon >0$ atitinka toks elementas $x_{N},$ kad $x_{N}>{\overline {x}}-\varepsilon$ ir $x_{N}\leq {\overline {x}}$ (bet koks elementas $x_{n}$ ne didesnis už tikslųjį viršutinį rėžį ${\overline {x}},$ t. y. $x_{n}\leq {\overline {x}}$). Iš šitų nelygybių gauname nelygybes $0\leq {\overline {x}}-x_{N}<\varepsilon .$ Kadangi $\{x_{n}\}$ – nemažėjanti seka, tai nelygybės $x_{N}\leq x_{n}\leq {\overline {x}}$ yra teisingos, kai $n\geq N.$ Iš to išplaukia, kad $0\leq {\overline {x}}-x_{n}\leq {\overline {x}}-x_{N},\;$ kai $n\geq N.$ Anksčiau pabrėžėme, kad ${\overline {x}}-x_{N}<\varepsilon ,$ todėl, kai $n\geq N,$ teisingos nelygybės $0\leq {\overline {x}}-x_{n}<\varepsilon ,$ iš kurių gauname nelygybę $|x_{n}-{\overline {x}}|<\varepsilon .$ Vadinasi, įsitikinome, kad ${\overline {x}}$ yra sekos $\{x_{n}\}$ riba. Teorema įrodyta.

1 pastaba. Monotoninės sekos aprėžtumas yra būtinas ir pakankamas jos konvergavimo požymis.

Iš tikrųjų, jei monotoninė seka yra aprėžta, tai pagal 3.15 teoremą ji konverguoja; jei monotoninė seka konverguoja, tai, remiantis 3.8 teorema, ji yra aprėžta.

2 pastaba. Konverguojanti seka gali ir nebūti monotoninė. Pavyzdžiui, seka $\{x_{n}\},$ apibrėžta formule $x_{n}={\frac {(-1)^{n}}{n}},$ konverguoja, o jos riba yra skaičius 0. Tačiau ši seka nėra monotoninė, nes po teigiamo elemento eina neigiamas, o po neigiamo – teigiamas.

3 pastaba. Jei $\{x_{n}\}$ yra nemažėjanti ir aprėžta seka, o ${\overline {x}}$ – tos sekos riba, tai visi jos elementai tenkina nelygybę $x_{n}\leq {\overline {x}}.$ Nedidėjančios ir aprėžtos sekos $\{x_{n}\},$ konverguojančios į ${\underline {x}},$ elementai tenkina nelygybę ${\underline {x}}\leq x_{n}.$ Kad šie teiginiai teisingi, įsitikinome, įrodinėdami 3.15 teoremą.

3.15 teoremos išvada. Sakykime, duota begalinė segmentų sistema $[a_{1};\;b_{1}],\;[a_{2};\;b_{2}],\;...,\;[a_{n};\;b_{n}],\;...\;$ Jei kiekvienas segmentas yra pirmesniajame*, o skirtumas $b_{n}-a_{n}$ (jį vadinsime segmento $[a_{n};\;b_{n}]$ ilgiu) artėja prie nulio, kai $n\to \infty$ (segmentų sistemą, turinčią šią savybę, vadinsime susitraukiančia), tai egzistuoja vienintelis taškas c, priklausantis visiems tos sistemos segmentams.

Įrodymas. Pirmiausia pastebėsime, kad taškas c, priklausantis visiems segmentams, gali būti tik vienas. Iš tikrųjų, jei būtų dar vienas taškas d, priklausantis visiems segmentams, tai segmentas ** [c; d] priklausytų visiems segmentams $[a_{n};\;b_{n}].$ Bet tokiu atveju su bet kokiu numeriu n būtų teisingos nelygybės $b_{n}-a_{n}\geq d-c>0;$ tai neįmanoma, nes $b_{n}-a_{n}\to 0,$ kai $n\to \infty .$ Dabar įsitikinsime, kad taškas c, priklausantis visiems segmentams $[a_{n};\;b_{n}],$ egzistuoja. Kadangi segmentų sistema yra susitraukianti, tai kairiųjų galų seka $\{a_{n}\}$ yra nemažėjanti, o dešiniųjų galų seka $\{b_{n}\},$ – nedidėjanti. Kadangi abi tos sekos yra aprėžtos (visi sekų $\{a_{n}\}$ ir $\{b_{n}\}$ elementai priklauso segmentui $[a_{1};\;b_{1}]$), tai pagal 3.15 teoremą jos konverguoja. Iš to, kad skirtumas $b_{n}-a_{n}$ nyksta, išplaukia, jog minimosios sekos turi bendrą ribą. Tą ribą žymėkime raide c. Remiantis 3 pastaba, su bet kokiu numeriu n teisingos šios nelygybės: $a_{n}\leq c\leq b_{n},$ t. y. taškas c priklauso visiems segmentams $[a_{n};\;b_{n}].$

_______________

* Tai reiškia, kad $a_{n-1}\leq a_{n}<b_{n}\leq b_{n-1}.$

** Siekdami konkretumo, tariame, kad d > c.

3. Keli konverguojančių monotoninių sekų pavyzdžiai.

Išnagrinėsime kelias sekas, kurių ribas apskaičiuosime, remdamiesi 3.15 teorema apie monotoninės sekos ribą. Be to, šiame skirsnyje susipažinsime su vienu bendru sekos ribos ieškojimo metodu, kai seka apibrėžiama rekurentine formule (Rekurentinė formulė (lot. recurrens – grįžtas) – formulė, pagal kurią $(n+1)$-ąjį sekos elementą galima išreikšti jos pirmųjų n elementų reikšmėmis.).

1 pavyzdys. Tirsime seką $\{x_{n}\},$ kurios elementas $x_{n}$ lygus

$x_{n}={\sqrt {a+{\sqrt {a+{\sqrt {a+...+{\sqrt {a}}}}}}}},\;\;a>0$

(n radikalų). Tą pačią seką galima, savaime aišku, nusakyti šitokia rekurentine formule:

$x_{1}={\sqrt {a}},\quad x_{n}={\sqrt {a+x_{n-1}}}.$

Norėdami įsitikinti, kad sekos $\{x_{n}\}$ riba egzistuoja, įrodysime, jog ta seka yra didėjanti ir aprėžta. Kad ta seka didėja, matyti tiesiog. Įsitikinsime, kad seka $\{x_{n}\}$ aprėžta iš viršaus skaičiumi A, kai A – didedsnysis iš dviejų skaičių a ir 2. Jei $x_{n}\leq a,$ tai teiginys jau įrodytas. Jei $x_{n}>a,$ tai, dešinėje nelygybės $x_{n}^{2}=a+x_{n-1}\leq a+x_{n}\;$ pusėje skaičių a pakeitę didesniu už jį skaičiumi $x_{n},$ gauname $x_{n}^{2}<2x_{n},$ o iš čia $x_{n}<2.$ Taigi įrodėme, kad seka $\{x_{n}\}$ aprėžta iš viršaus. Pagal 3.15 teoremą ta seka turi ribą. Pažymėkime tą ribą raide c. Savaime aišku, $c>0.$ Iš rekurentinės formulės turime lygybę

$x_{n}^{2}=a+x_{n-1},$

iš kurios matyti, kad sekos $\{x_{n}^{2}\}$ ir $\{a+x_{n-1}\}$ sutampa. Todėl jos turi bendrą ribą. Kadangi pirmosios sekos riba lygi $c^{2},$ o antrosios $a+c$ (kai n varo į begalybę skaičiams $x_{n}$ ir $x_{n-1},$ tai $x_{n}$ praktiškai nesiskiria nuo $x_{n-1}$ ir abu yra labai arti c reikšmės (abu labai panašūs į c)), tai $c^{2}=a+c.$ Iš čia, turėdami mintyje, kad $c>0,$ randame c:

$c^{2}=a+c,$

$c^{2}-c-a=0,$

$c_{1,2}={\frac {-(-1)\pm {\sqrt {(-1)^{2}-4\cdot (-a)}}}{2}}={\frac {1\pm {\sqrt {1+4a}}}{2}};$

kadangi $c>0,$ tai

$c={\frac {1+{\sqrt {1+4a}}}{2}}.$

2 pavyzdys. Tirsime seką $\{x_{n}\},$ kuria dažniausiai naudojamasi, apskaičiuojant teigiamo skaičiaus a kvadratinę šankį elektronine skaičiavimo mašina. Šita seka apibrėžiama tokia rekurentine formule:

$x_{n+1}={\frac {1}{2}}{\Big (}x_{n}+{\frac {a}{x_{n}}}{\Big )},\quad n=1,\;2,\;3,\;...$

Pirmuoju elementu $x_{1}$ čia galima laikyti bet kokį teigiamą skaičių.

Įrodysime, kad ta seka konverguoja, o jos riba yra skaičius ${\sqrt {a}}.$ Pirmiausia įsitikinsime, kad seka $\{x_{n}\}$ turi ribą. Tam užtenka įrodyti, kad seka $\{x_{n}\}$ yra aprėžta iš apačios ir kad, pradedant antruoju elementu, ji nedidėja. Iš pradžių įsitikinsime, kad seka $\{x_{n}\}$ aprėžta iš apačios. Sąlygoje pasakyta, jog $x_{1}>0,$ o iš rekurentinės formulės, kai $n=1,$ gauname $x_{2}>0.$ Taip samprotaudami toliau, įsitikiname, kad visi $x_{n}$ yra teigiami.

Dabar įrodysime, kad visi $x_{n},$ kai $n\geq 2,$ tenkina nelygybę $x_{n}\geq {\sqrt {a}}.$ Parašę rekurentinę formulę šitaip: $x_{n+1}={\frac {\sqrt {a}}{2}}{\Big (}{\frac {x_{n}}{\sqrt {a}}}+{\frac {\sqrt {a}}{x_{n}}}{\Big )},$ remsimės beveik savaime aiškia nelygybe $t+{\frac {1}{t}}\geq 2$*, kuri teisinga, kai $t>0$ (imame $t={\frac {x_{n}}{\sqrt {a}}}$). Gausime

$x_{n+1}={\frac {\sqrt {a}}{2}}{\Big (}{\frac {x_{n}}{\sqrt {a}}}+{\frac {\sqrt {a}}{x_{n}}}{\Big )}\geq {\frac {\sqrt {a}}{2}}\cdot 2={\sqrt {a}},$

kai $n\geq 1,$ t. y. $x_{n}\geq {\sqrt {a}},$ pradedant numeriu $n=2.$

Pagaliau įsitikinsime, kad seka $\{x_{n}\},$ kai $n\geq 2,$ nedidėja. Iš rekurentinės formulės matome, kad ${\frac {x_{n+1}}{x_{n}}}={\frac {1}{2}}{\Big (}1+{\frac {a}{x_{n}^{2}}}{\Big )},$ o iš čia, turėdami mintyje, kad $x_{n}\geq {\sqrt {a}},$ gauname ${\frac {x_{n+1}}{x_{n}}}\leq 1,$ arba $x_{n}\geq x_{n+1}$ (kai $n\geq 2$).

Kadangi seka $\{x_{n}\},$ kai $n\geq 2,$ nedidėja ir yra aprėžta iš apačios skaičiumi ${\sqrt {a}}$ (nes $x_{n}\geq {\sqrt {a}}$), tai ji turi ribą, ne mažesnę kaip ${\sqrt {a}}$ (žr. 3.15 ir 3.13 teoremą). Tą ribą pažymėkime raide c ir atsižvelgę į tai, kad $\lim _{n\to \infty }x_{n+1}=c,$ o $\lim _{n\to \infty }{\frac {1}{2}}{\Big (}x_{n}+{\frac {a}{x_{n}}}{\Big )}={\frac {1}{2}}{\Big (}c+{\frac {a}{c}}{\Big )},$ gauname lygybę

$c={\frac {1}{2}}{\Big (}c+{\frac {a}{c}}{\Big )}$**,

$c^{2}={\frac {1}{2}}{\Big (}c^{2}+a{\Big )},$

$2c^{2}=c^{2}+a,$

$2c^{2}-c^{2}=a,$

$c^{2}=a.$

Vadinasi, $c={\sqrt {a}}.$

1 pastaba. Spręsdami pateiktuosius pavyzdžius, naudojomės plačiai taikomu metodu sekos ribai skaičiuoti. Metodo esmė: iš pradžių įsitikinama, kad sekos riba egzistuoja, o paskui iš lygties, kuri gaunama iš rekurentinės formulės, vietoje $x_{n}$ ir $x_{n+1}$ įrašius ieškomąją sekos $\{x_{n}\}$ ribos reikšmę c, apskaičiuojama tos ribos skaitinė reikšmė.

2 pastaba. Rekurentinės formulės dažnai naudojamos skaičiavimo matematikoje, nes jas taikant daug kartų, atliekamos to paties tipo skaičiavimo operacijos, o tai labai patogu skaičiuojant elektroninėmis skaičiavimo mašinomis.

Išnagrinėtoji rekurentinė formulė aprašo, kaip įsitikinome, ${\sqrt {a}}$ skaičiavimo algoritmą (įrodėme, kad $\lim _{n\to \infty }x_{n}={\sqrt {a}}$).

Vėliau bus tiriamas sekos $\{x_{n}\}$ konvergavimo į ${\sqrt {a}}$ greitis. Įrodoma, kad, tinkamai pasirinkus pirmąjį artinį $x_{1},$ kai $a>1,$ jau ketvirtasis artinys $x_{4}$ nuo skaičiaus ${\sqrt {a}}$ skiriasi mažiau kaip $10^{-10}.$

3 pavyzdys. Įsitikinsime, kad sekos $\{c_{n}\},$ kai $c_{n}={\frac {x^{n+1}}{(n+1)!}},$ o x – bet koks fiksuotas skaičius, riba lygi nuliui. Kai n – pakankamai didelis natūralusis skaičius, ${\frac {|x|}{n+1}}<1.$ Todėl, pradedant kuriuo nors numeriu N, turi būti $|c_{n+1}|<|c_{n}|,$ nes $|c_{n+1}|={\frac {|x|^{n}}{n!}}\cdot {\frac {|x|}{n+1}}=|c_{n}|\cdot {\frac {|x|}{n+1}}.$

Vadinasi, pradedant numeriu N, seka $\{|c_{n}|\}$ mažėja. Kadangi, be to, ji aprėžta iš apačios (pavyzdžiui nuliu), tai pagal 3.15 teoremą seka $\{|c_{n}|\}$ konverguoja. Tarkime, kad c yra tos sekos riba. Iš lygybės $|c_{n+1}|=|c_{n}|\cdot {\frac {|x|}{n+1}}$ išplaukia, kad $c=0,$ nes sekos $\{|c_{n+1}|\}$ riba lygi c, o sekos $\{{\frac {|x|}{n+1}}\}$ – nuliui.

_____________________

* Tą nelygybę įrodinėjant, užtenka pastebėti, kad, kai $t>0,$ ji ekvivalenti nelygybei $t^{2}-2t+1\geq 0.$

** Ta lygybė gaunama iš rekurentinės formulės $x_{n+1}={\frac {1}{2}}{\Big (}x_{n}+{\frac {a}{x_{n}}}{\Big )}.$

Sekos, aproksimuojančios a šaknyje, konvergavimo greitis

Šio skyriaus 3 paragrafo 3 skirsnyje įrodėme, kad sekos $\{x_{n}\}$ apibrėžiamos rekurentine formule

$x_{n+1}={\frac {1}{2}}{\Big (}x_{n}+{\frac {a}{x_{n}}}{\Big )},\quad n=1,\;2,\;3,\;...,\quad (3.10)$

kai $a>0,$ o $x_{1}$ – bet koks teigiamas skaičius, riba lygi ${\sqrt {a}}.$ Skaičiaus ${\sqrt {a}}$ artiniu galima laikyti bet kurį tos sekos narį $x_{n+1}.$ Tokiu atveju, be abejo, reikia išsiaiškinti, koks turi būti įteracijų* skaičius n, kad ${\sqrt {a}}$ artinys būtų nurodyto tikslumo.

______________

* Iteracija (lot. iteratio – kartojimas) – kokios nors matematinės operacijos pakartojimas. Šiuo atveju viena iteracija yra $x_{n+1}$ apskaičiavimas, kai žinomas $x_{n},$ pagal rekurentinę (3.10) formulę.

______________

Nagrinėsime seką $x_{n},$ apibrėžtą rekurentine (3.10) formule. Tos sekos elementą $x_{n}$ vadinsime skaičiaus $\gamma ={\sqrt {a}}\;$ n-uoju artiniu. Skaičių

$\varepsilon _{n}={\frac {x_{n}-\gamma }{\gamma }}\quad (3.11)$

vadinsime n-ojo artinio santykine paklaida.

Įrodysime teiginį apie santykinės paklaidos $\varepsilon _{n+1}$ įvertinimą, remiantis pirmojo artinio santykine paklaida $\varepsilon _{1}.$

Jei $x_{1}$ taip parinktas, kad $|\varepsilon _{1}|<{\frac {1}{2}},$ tai su bet kokiu $n\geq 1$ teisinga šitokia nelygybė:

$0<\varepsilon _{n+1}<\varepsilon _{1}^{2^{n}}.\quad (3.12)$

[ $0\leq \varepsilon _{n+1}\leq \varepsilon _{1}^{2^{n}},$ kai $\varepsilon _{1}$ gali būti lygus nuliui.]

Įrodymas. Iš (3.11) formulės

$x_{n}=\gamma (1+\varepsilon _{n}).\quad (3.13)$

Remdamiesi (3.10) ir (3.13) formulėmis bei lygybe ${\frac {a}{\gamma }}=\gamma ,$ gauname

$x_{n+1}={\frac {1}{2}}{\Big (}x_{n}+{\frac {a}{x_{n}}}{\Big )}={\frac {1}{2}}\left[\gamma (1+\varepsilon )+{\frac {a}{\gamma (1+\varepsilon _{n})}}\right]={\frac {1}{2}}\left[{\frac {a(1+\varepsilon _{n})^{2}+a}{\gamma (1+\varepsilon _{n})}}\right]=$

$={\frac {1}{2}}\left[{\frac {a(1+2\varepsilon _{n}+\varepsilon _{n}^{2})+a}{\gamma (1+\varepsilon _{n})}}\right]={\frac {a}{\gamma }}\left[{\frac {(1+2\varepsilon _{n}+\varepsilon _{n}^{2})+1}{2(1+\varepsilon _{n})}}\right]=\gamma \left[{\frac {2+2\varepsilon _{n}+\varepsilon _{n}^{2}}{2(1+\varepsilon _{n})}}\right]=\gamma \left[1+{\frac {\varepsilon _{n}^{2}}{2(1+\varepsilon _{n})}}\right].$

Kadangi $x_{n+1}=\gamma (1+\varepsilon _{n+1}),$ tai savaime aišku, kad

$\gamma (1+\varepsilon _{n+1})=\gamma \left[1+{\frac {\varepsilon _{n}^{2}}{2(1+\varepsilon _{n})}}\right],$

$\varepsilon _{n+1}={\frac {1}{2(1+\varepsilon _{n})}}\cdot \varepsilon _{n}^{2}.\quad (3.14)$

Pagal sąlygą $|\varepsilon _{1}|<{\frac {1}{2}}.$ Iš to išplaukia nelygybės $0<{\frac {1}{2(1+\varepsilon _{n})}}<1.$ Bet tada iš (3.14) lygybės, kai $n=1,$ išplaukia nelygybė $\varepsilon _{2}>0.$ Panašiai naudodamiesi (3.14) lygybe, kai $n=1,\;2,\;...,$ įsitikiname, kad visi $\varepsilon _{n+1}$ yra neneigiami ($n\geq 1$).

Iš (3.14) lygybės, nelygybių $0<{\frac {1}{2(1+\varepsilon _{n})}}<1$ ir iš to, kad visi $\varepsilon _{n},$ kai $n>1,$ yra neneigiami, išplaukia nelygybė $\varepsilon _{n+1}<\varepsilon _{n}^{2},$ kai $n\geq 1.$ Iš čia gauname dešiniąją (3.12) nelygybę. Teiginys įrodytas.

[$|\varepsilon _{1}|<{\frac {1}{2}},\;$ $\varepsilon _{2}<\varepsilon _{1}^{2},\;$ $\varepsilon _{3}<\varepsilon _{2}^{2}<(\varepsilon _{1}^{2})^{2}=\varepsilon _{1}^{4},\;$ $\varepsilon _{4}<\varepsilon _{3}^{2}<(\varepsilon _{1}^{4})^{2}=\varepsilon _{1}^{8}=\varepsilon _{1}^{2^{3}}$ ir taip toliau.]

Iš (3.12) nelygybių matyti, kad ${\sqrt {a}}$ artinio po n iteracijų santykinė paklaida įvertinama, remiantis pirmojo artinio $x_{1}$ santykine paklaida $\varepsilon _{1}$ ir iteracijų skaičiumi n. Vėliau įsitikinsime, jog tuo atveju, kai $a>1$*, pirmąjį artinį galima pasirinkti taip, kad santykinės paklaidos $\varepsilon _{1}$ modulis nebūtų didesnis už 0.05. Savaime aišku, taip pasirinkus $x_{1},$ santykinė paklaida tenkins įrodytojo teiginio sąlygas. Tuomet bus aišku ir koks turi būti iteracijų skaičius n, kad ${\sqrt {a}}$ artinio santykinė paklaida nebūtų didesnė už duotąjį skaičių $\varepsilon$: tą skaičių n galima rasti iš nelygybės* (ta nelygybė gaunama tiesiog iš (3.12) nelygybių)

$(0.05)^{2^{n}}<\varepsilon .\quad (3.15)$

Taigi tarkime, kad $a>1.$ Skaičių a išreikšime šitaip:

$a=2^{2k+i}M;\quad (3.16)$

čia k – sveikasis neneigiamas skaičius, skaičius i lygus arba nuliui, arba vienetui, o skaičius M tenkina sąlygą

$1\leq M<2.\quad (3.17)$

Pabrėšime, kad skaičius a išreiškiamas (3.16) pavidalu vieninteliu budu.

Pirmąjį artinį $x_{1}$ pasirinkime šitaip:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}.\quad (3.18)$

Įsitikinsime, kad tuo atveju, kai M – bet koks skaičius, tenkinantis (3.17) sąlygas, pirmojo ${\sqrt {a}}$ artinio $x_{1},$ rasto pagal (3.18) formulę, santikinės paklaidos $\varepsilon _{1}$ modulis ne didesnis kaip 0.05.

Įrodinėdami šį teiginį, remsimės tikslia santykinės paklaidos išraiška $\varepsilon _{1}={\frac {x_{1}-\gamma }{\gamma }}.$ Iš (3.16) lygybės gauname $\gamma ={\sqrt {a}}=2^{k}{\sqrt {2^{i}M}};$ todėl iš $\varepsilon _{1}$ išraiškos ir iš (3.18) formulės

$\varepsilon _{1}={\frac {x_{1}-\gamma }{\gamma }}={\frac {2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}-2^{k}{\sqrt {2^{i}M}}}{2^{k}{\sqrt {2^{i}M}}}}={\frac {{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}-{\sqrt {2^{i}M}}}{\sqrt {2^{i}M}}}.\quad (3.19)$

3.4 pav.

Alternatyvus 3.4 pav. Oy ašimi grafikas ištemptas tik maždaug 2 kartus daugiau nei Ox ašimi.

Kadangi skaičius i lygus arba nuliui, arba vienetui, o $M\geq 1,$ tai ${\sqrt {2^{i}M}}\geq 1.$ Todėl iš (3.19) lygybės išplaukia nelygybė

$|\varepsilon _{1}|\leq |{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}-{\sqrt {2^{i}M}}|.\quad (3.20)$

Pažymėkime ${\sqrt {2^{i}M}}$ raide X. Kadangi $1\leq M<2,$ o i lygus arba nuliui, arba vienetui, tai visos galimos X reikšmės priklauso intervalui [1; 2):

$1\leq X<2.\quad (3.21)$

Naudodami naujai įvestą žymėjimą, (3.20) nelygybę perrašysime šitaip:

$|\varepsilon _{1}|\leq |{\frac {1}{3}}X^{2}-X+{\frac {17}{24}}|.\quad (3.22)$

Pagal (3.22) nelygybę maksimalioji $|\varepsilon _{1}|$ reikšmė ne didesnė už maksimaliąją reiškinio $|{\frac {1}{3}}X^{2}-X+{\frac {17}{24}}|$ reikšmę, kai X tenkina (3.21) sąlygas. Kad būtų aiškiau, išnagrinėkime funkcijos $f(X)={\frac {1}{3}}X^{2}-X+{\frac {17}{24}}$ grafiką. Iš elementariosios matematikos kurso žinome, kad šios funkcijos grafikas yra parabolė (3.4 pav.)**, kurios viršūnė atitinka $X={\frac {3}{2}}.$

[ $f'(X)={\Big (}{\frac {1}{3}}X^{2}-X+{\frac {17}{24}}{\Big )}'={\frac {2}{3}}X-1,$

${\frac {2}{3}}X-1=0,\;\;\;{\frac {2}{3}}X=1,\;\;\;X={\frac {3}{2}}.$ ]

Kadangi $f(1)=f(2)={\frac {1}{24}},\;$ o $\;f{\Big (}{\frac {3}{2}}{\Big )}=-{\frac {1}{24}},\;$ tai X reikšmės, tenkinančios (3.21) sąlygas, atitinka $f(X)$ reikšmėms, esančioms tarp $-{\frac {1}{24}}$ ir ${\frac {1}{24}}.$ Kitaip sakant,

$|f(X)|=|{\frac {1}{3}}X^{2}-X+{\frac {17}{24}}|\leq {\frac {1}{24}}.$

Iš paskutinės ir (3.22) nelygybių gauname reikiamą $\varepsilon _{1}$ įvertinimą:

$|\varepsilon _{1}|\leq {\frac {1}{24}}<0.05.$

Pastaba. Jei nurodyta, kad santykinė paklaida $\varepsilon$ lygi $10^{-10},$ tai, norint tokiu tikslumu apskaičiuoti bet kokio skaičiaus $a>1$ kvadratinę šaknį, reikia, suradus $x_{1}$ pagal (3.18) formulę, atlikti tik tris iteracijas ($n=3$), nes

$(0.05)^{2^{3}}=0.05^{8}=0.0000000000390625=3.90625\cdot 10^{-11}<10^{-10}.$

[ $(1/24)^{8}=(0.041(6))^{8}\approx 9.0846890375538\cdot 10^{-12}<10^{-11}<10^{-10}.$ Pasirodo, santykinė paklaida $\varepsilon _{4}$ visada bus mažesnė už $10^{-11}$ po 3 iteracijų ($n=3$). ]

_______________

* Jei $a<1,$ tai $a={\frac {1}{b}}$ ir $b>1.$ Tada ${\sqrt {a}}={\frac {1}{\sqrt {b}}}.$

** 3.4 paveiksle ašies Oy mastelis 20 kartų didesnis už ašies Ox mastelį.

Pavyzdžiai

• Apskaičiuosime ${\sqrt {a}}={\sqrt {7}}.$ Iš (3.16) formulės $a=2^{2k+i}M.$ Matome, kad $2^{2k+i}$ turi būti lygus 4, t. y. $2^{2k+i}=4.$ Tada aišku, kad $k=1,\;\;i=0.$ Toliau randame M:

$7=2^{2k+i}M,$

$7=4M,$

$M={\frac {7}{4}}.$

Toliau randame $x_{1}$ pagal (3.18) formulę:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{1}{\Big (}{\frac {1}{3}}\cdot 2^{0}\cdot {\frac {7}{4}}+{\frac {17}{24}}{\Big )}=$

$=2{\Big (}{\frac {1}{3}}\cdot {\frac {7}{4}}+{\frac {17}{24}}{\Big )}={\frac {7}{6}}+{\frac {17}{12}}={\frac {14+17}{12}}={\frac {31}{12}}=$

= 2.58(3).

Skaičių $x_{2}$ randame pagal (3.10) rekurentinę formulę:

$x_{1+1}={\frac {1}{2}}{\Big (}x_{1}+{\frac {a}{x_{1}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {31}{12}}+{\frac {7}{\frac {31}{12}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {31}{12}}+{\frac {84}{31}}{\Big )}={\frac {1}{2}}\cdot {\frac {31^{2}+84\cdot 12}{12\cdot 31}}={\frac {1}{2}}\cdot {\frac {961+1008}{372}}={\frac {1969}{744}}=$

= 2.646505376344086021505376344086.

Toliau rasime $x_{3}$:

$x_{3}={\frac {1}{2}}{\Big (}x_{2}+{\frac {a}{x_{2}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {1969}{744}}+{\frac {7}{\frac {1969}{744}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {1969}{744}}+{\frac {5208}{1969}}{\Big )}={\frac {1}{2}}\cdot {\frac {1969^{2}+5208\cdot 744}{744\cdot 1969}}={\frac {7751713}{2929872}}=$

= 2.6457514184920023809913880196814.

[Atėmę tikslią ${\sqrt {7}}$ reikšmę iš $x_{3}$ reikšmės gauname:

2.6457514184920023809913880196814 - 2.6457513110645905905016157536393 = 1.0742741179048977226604214545141e-7.]

Toliau rasime $x_{4}$:

$x_{4}=x_{3+1}={\frac {1}{2}}{\Big (}x_{3}+{\frac {a}{x_{3}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {7751713}{2929872}}+{\frac {7}{\frac {7751713}{2929872}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {7751713}{2929872}}+{\frac {20509104}{7751713}}{\Big )}=$

$={\frac {1}{2}}\cdot {\frac {7751713^{2}+20509104\cdot 2929872}{2929872\cdot 7751713}}={\frac {1}{2}}\cdot {\frac {60,089,054,434,369+60,089,049,554,688}{2929872\cdot 7751713}}={\frac {120,178,103,989,057}{45,423,053,741,472}}=$

= 120178103989057/45423053741472 = 2.6457513110645927714791488926518.

Atėmę tikslią ${\sqrt {7}}$ reikšmę iš gautos $x_{4}$ reikšmės, gauname:

2.6457513110645927714791488926518 - 2.6457513110645905905016157536393 = 2.1809775331390125063383535564728e-15.

Tai yra gavome paklaidą $\varepsilon \approx 2.180977533139\cdot 10^{-15}<10^{-10}.$

Santykinė paklaida yra tokia:

$\varepsilon _{4}={\frac {x_{4}-\gamma }{\gamma }}=$

(2.6457513110645927714791488926518 - 2.6457513110645905905016157536393)/2.6457513110645905905016157536393 =

= 8.2433202395785127055094193690687e-16.

Taigi $\varepsilon _{4}\approx 8.2433202395785\cdot 10^{-16}<10^{-10}.$

• Apskaičiuosime ${\sqrt {a}}={\sqrt {39}}.$ Iš (3.16) formulės $a=2^{2k+i}M.$ Pagal (3.17) $1\leq M<2.$ Todėl $2^{2k+i}=32.$ Tada $2k+i=5.$ Iš čia gauname, kad $k=2,\;\;i=1.$

Toliau randame M:

$a=2^{2k+i}M,$

$39=2^{2\cdot 2+1}M,$

$39=2^{5}M,$

$M={\frac {39}{32}}.$

Toliau randame $x_{1}$ pagal (3.18) formulę:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{2}{\Big (}{\frac {1}{3}}\cdot 2^{1}\cdot {\frac {39}{32}}+{\frac {17}{24}}{\Big )}=4{\Big (}{\frac {2}{3}}\cdot {\frac {39}{32}}+{\frac {17}{24}}{\Big )}={\frac {312}{96}}+{\frac {17}{6}}={\frac {312+17\cdot 16}{96}}={\frac {584}{96}}={\frac {73}{12}}=$

= 6.08(3).

[ ${\sqrt {39}}=$ 6.2449979983983982058468931209398. ]

Skaičių $x_{2}$ randame pagal (3.10) rekurentinę formulę:

$x_{2}={\frac {1}{2}}{\Big (}x_{1}+{\frac {a}{x_{1}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {73}{12}}+{\frac {39}{\frac {73}{12}}}{\Big )}=$

= 0.5*(6.0833333333333333333333333333333 + 6.4109589041095890410958904109589) =

= 0.5 * 12.494292237442922374429223744292 = 6.2471461187214611872146118721461.

Toliau rasime $x_{3}$:

$x_{3}={\frac {1}{2}}{\Big (}x_{2}+{\frac {a}{x_{2}}}{\Big )}=$

= 0.5*(6.2471461187214611872146118721461 + 39/6.2471461187214611872146118721461) =

= 6.2449983677207123204232036062421.

[Atėmę tikslią ${\sqrt {39}}$ reikšmę iš $x_{3}$ reikšmės, gauname:

6.2449983677207123204232036062421 - 6.2449979983983982058468931209398 = 3.693223141145763104853022918526e-7.]

Toliau rasime $x_{4}$:

$x_{4}=x_{3+1}={\frac {1}{2}}{\Big (}x_{3}+{\frac {a}{x_{3}}}{\Big )}=$

= 0.5*(6.2449983677207123204232036062421 + 39/6.2449983677207123204232036062421) =

= 6.2449979983984091265040071378015.

Atėmę tikslią ${\sqrt {39}}$ reikšmę iš $x_{4}$ reikšmės, gauname:

6.2449979983984091265040071378015 - 6.2449979983983982058468931209398 = 1.0920657114016861678930991073622e-14.

Taigi, gavome paklaidą $1.0920657114\cdot 10^{-14}.$ O santykinė paklaida $\varepsilon _{4}$ lygi

$\varepsilon _{4}={\frac {x_{4}-\gamma }{\gamma }}=$

= (6.2449979983984091265040071378015 - 6.2449979983983982058468931209398)/6.2449979983983982058468931209398 =

= 1.7487046620059110065547770263338e-15.

T. y. $\varepsilon _{4}\approx 1.748704662\cdot 10^{-15}<10^{-10}.$

• Apskaičiuosime ${\sqrt {a}}={\sqrt {2}}.$ Iš (3.16) formulės $a=2^{2k+i}M.$ Pagal (3.17) $1\leq M<2.$ Todėl $2^{2k+i}=2.$ Tada $2k+i=1.$ Iš čia gauname, kad $k=0,\;\;i=1.$

Kadangi $2^{2k+i}=2,$ tai $M=1.$

Toliau randame $x_{1}$ pagal (3.18) formulę:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{0}{\Big (}{\frac {1}{3}}\cdot 2^{1}\cdot 1+{\frac {17}{24}}{\Big )}={\frac {2}{3}}+{\frac {17}{24}}={\frac {2\cdot 8+17}{24}}={\frac {33}{24}}=1.375.$

Skaičių $x_{2}$ randame pagal (3.10) rekurentinę formulę:

$x_{2}={\frac {1}{2}}{\Big (}x_{1}+{\frac {a}{x_{1}}}{\Big )}={\frac {1}{2}}{\Big (}1.375+{\frac {2}{1.375}}{\Big )}=$

= 0.5*(1.375 + 1.4545454545454545454545454545455) = 0.5 * 2.8295454545454545454545454545455 =

= 1.4147727272727272727272727272727.

Toliau rasime $x_{3}$:

$x_{3}={\frac {1}{2}}{\Big (}x_{2}+{\frac {a}{x_{2}}}{\Big )}=$

= 0.5*(1.4147727272727272727272727272727 + 2/1.4147727272727272727272727272727) =

= 1.4142136728733114275282949981745.

[1.4142136728733114275282949981745 - 1.4142135623730950488016887242097 = 0.0000001105002163787266062739648 $\approx 1.105002163787266\cdot 10^{-7}.$]

Toliau rasime $x_{4}$:

$x_{4}={\frac {1}{2}}{\Big (}x_{3}+{\frac {a}{x_{3}}}{\Big )}=$

= 0.5*(1.4142136728733114275282949981745 + 2/1.4142136728733114275282949981745) =

= 1.4142135623730993657935457387426.

Atėmę tikslią ${\sqrt {2}}$ reikšmę iš $x_{4}$ reikšmės, gauname:

1.4142135623730993657935457387426 - 1.4142135623730950488016887242097 = 4.316991857014532913283559342618e-15.

Santykinė paklaida yra tokia:

$\varepsilon _{4}={\frac {x_{4}-\gamma }{\gamma }}=$

=(1.4142135623730993657935457387426 - 1.4142135623730950488016887242097)/1.4142135623730950488016887242097 =

= 3.0525742164220827729222598916484e-15.

Gavome, kad $\varepsilon _{4}\approx 3.052574216422\cdot 10^{-15}<10^{-10}.$

Šiaip, skaičiuojant su kompiuteriu, santikinę paklaidą $\varepsilon _{4}$ (arba $\varepsilon _{n+1}$) reikia padauginti iš 10, nes procesoriaus x87 FPU be eksponentės dirba su skaičiais nuo 1.00000000000000001 iki 9.99999999999999999. Tada nebus prarastas tikslumas, jei įvertinsime, kad FPU paklaida yra $\delta <10\varepsilon _{n+1}.$

Šaknies traukimas iš didelių skaičių

Bet kokį skaičių a ($a>1$) galima užrašyti taip:

$a=p\cdot 10^{j}\cdot 10^{v}.\quad (1)$

Čia $1<p<10;\;$ skaičius v gali būti nulis arba teigiamas lyginis skaičius (v gali būti 0, 2, 4, 6, 8, 10 ir t. t.); skaičius j lygus arba nuliui arba vienetui.

Jeigu ištrauksime kvadratinę šaknį iš a, tai gausime:

${\sqrt {a}}={\sqrt {p}}\cdot {\sqrt {10^{j}}}\cdot 10^{v \over 2}.\quad (2)$

Jegiu v nelygus nuliui, tai $10^{v/2}$ bus sveikasis skaičius su dvigubai mažiau nulių nei pas $10^{v}.$ Jeigu $j=1,$ tai tereikia apskaičiuoti ${\sqrt {10}},$ kurį bus galima visada taikyti (2) formulėje. Vadinasi, tereikia mokėti greitai ištraukti šaknį iš skaičiaus p.

1) Jeigu $1<p<2,\;$ tai pagal (3.16) formulę

$p=2^{2k+i}M;$

iš čia $k=0,\;\;i=0,\;\;M=p.$

Tada pagal (3.18) formulę randame $x_{1}$:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{0}{\Big (}{\frac {1}{3}}\cdot 2^{0}\cdot p+{\frac {17}{24}}{\Big )}={\frac {1}{3}}p+{\frac {17}{24}}.$

Ir toliau, darant iteracijas, apskaičiuojamas $x_{2},\;x_{3},\;x_{4},\;...$ pagal (3.10) formulę:

$x_{n+1}={\frac {1}{2}}{\Big (}x_{n}+{\frac {p}{x_{n}}}{\Big )},\quad n=1,\;2,\;3,\;...$

2) Jeigu $2\leq p<4,\;$ tai pagal (3.16) formulę

$p=2^{2k+i}M;$

iš čia $k=0,\;\;i=1.$ Tada $p=2M,\;\;M={\frac {p}{2}}.$

Pagal (3.18) formulę randame $x_{1}$:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{0}{\Big (}{\frac {1}{3}}\cdot 2^{1}\cdot {\frac {p}{2}}+{\frac {17}{24}}{\Big )}={\frac {1}{3}}p+{\frac {17}{24}}.$

Matome, kad šiuo atveju artinio $x_{1}$ formulė tokia pati kaip ir 1) atveju.

Toliau skaičiuojama kaip ir pirmu atveju.

3) Jeigu $4\leq p<8,\;$ tai pagal (3.16) formulę

$p=2^{2k+i}M;$

iš čia $k=1,\;\;i=0.$ Tada $p=4M,\;\;M={\frac {p}{4}}.$

Pagal (3.18) formulę randame $x_{1}$:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{1}{\Big (}{\frac {1}{3}}\cdot 2^{0}\cdot {\frac {p}{4}}+{\frac {17}{24}}{\Big )}={\frac {1}{6}}p+{\frac {17}{12}}.$

Toliau skaičiuojama kaip ir pirmu atveju.

4) Jeigu $8\leq p<10,\;$ tai pagal (3.16) formulę

$p=2^{2k+i}M;$

iš čia $k=1,\;\;i=1.$ Tada $p=8M,\;\;M={\frac {p}{8}}.$

Pagal (3.18) formulę randame $x_{1}$:

$x_{1}=2^{k}{\Big (}{\frac {1}{3}}\cdot 2^{i}M+{\frac {17}{24}}{\Big )}=2^{1}{\Big (}{\frac {1}{3}}\cdot 2^{1}\cdot {\frac {p}{8}}+{\frac {17}{24}}{\Big )}={\frac {1}{6}}p+{\frac {17}{12}}.$

Matome, kad šiuo atveju artinio $x_{1}$ formulė tokia pati kaip ir 3) atveju.

Toliau skaičiuojama kaip ir pirmu atveju.

Apibendrinant, skaičiaus ${\sqrt {p}}\,$ ($1<p<10$) apskaičiavimui, taikoma $x_{1}$ reikšmė

$x_{1}={\frac {1}{3}}p+{\frac {17}{24}},\quad (3)$

kai $1<p<4;$

ir taikoma $x_{1}$ reikšmė

$x_{1}={\frac {1}{6}}p+{\frac {17}{12}},\quad (4)$

kai $4\leq p<10.$

• Apskaičiuosime ${\sqrt {a}}={\sqrt {1.5\cdot 10^{8}}}.$ Čia $p=1.5,\;\;j=0,\;\;v=8.$ Tada pagal (1) formulę

$a=p\cdot 10^{j}\cdot 10^{v}=1.5\cdot 10^{0}\cdot 10^{8}.$

Pagal (2) formulę

${\sqrt {a}}={\sqrt {p}}\cdot {\sqrt {10^{j}}}\cdot 10^{v \over 2}={\sqrt {p}}\cdot 10^{8 \over 2}={\sqrt {p}}\cdot 10^{4}.$

Toliau randame $x_{1}$ skaičiui $p=1.5$:

$x_{1}={\frac {1}{3}}p+{\frac {17}{24}}={\frac {1}{3}}\cdot {\frac {3}{2}}+{\frac {17}{24}}={\frac {1}{2}}+{\frac {17}{24}}={\frac {12+17}{24}}={\frac {29}{24}}=$

= 1.2083333333333333333333333333333.

[ ${\sqrt {1.5}}=$ 1.2247448713915890490986420373529. ]

[1.2247448713915890490986420373529 - 1.2083333333333333333333333333333 = 0.0164115380582557157653087040196; tai yra ne santykinė paklaida, o reali paklaida.]

Toliau pagal (3.10) formulę apskaičiuojame $x_{2}$:

$x_{2}=x_{1+1}={\frac {1}{2}}{\Big (}x_{1}+{\frac {p}{x_{1}}}{\Big )}={\frac {1}{2}}{\Big (}{\frac {29}{24}}+{\frac {1.5}{\frac {29}{24}}}{\Big )}=$

= 0.5*(1.2083333333333333333333333333333 + 1.5/1.2083333333333333333333333333333) = 1.2248563218390804597701149425287.

Toliau pagal (3.10) formulę apskaičiuojame $x_{3}$:

$x_{3}=x_{2+1}={\frac {1}{2}}{\Big (}x_{2}+{\frac {p}{x_{2}}}{\Big )}=$

= 0.5*(1.2248563218390804597701149425287 + 1.5/1.2248563218390804597701149425287) = 1.2247448764620622240199548319682.

[1.2247448764620622240199548319682 - 1.2247448713915890490986420373529 = 0.0000000050704731749213127946153 = $5.07047\cdot 10^{-9}.$]

Toliau pagal (3.10) formulę apskaičiuojame $x_{4}$:

$x_{4}=x_{3+1}={\frac {1}{2}}{\Big (}x_{3}+{\frac {p}{x_{3}}}{\Big )}=$

= 0.5*(1.2247448764620622240199548319682 + 1.5/1.2247448764620622240199548319682) = 1.2247448713915890595945823395751.

Iš ką tik gautos $x_{4}$ reikšmės atėmę ${\sqrt {1.5}}$ reikšmę, gauname:

1.2247448713915890595945823395751 - 1.2247448713915890490986420373529 = 0.0000000000000000104959403022222 = $1.049594\cdot 10^{-17}.$

Santykinė paklaida yra tokia:

$\varepsilon _{4}={\frac {x_{4}-\gamma }{\gamma }}=$

(1.2247448713915890595945823395751 - 1.2247448713915890490986420373529)/1.2247448713915890490986420373529 =

= 8.5698993703859497564530032133562e-18.

Taigi $\varepsilon _{4}\approx 8.56989937\cdot 10^{-18}<10^{-10}.$

Iš anksčiau apskaičiuoto

${\sqrt {a}}={\sqrt {p}}\cdot 10^{4}=x_{4}\cdot 10^{4}=$

= 1.2247448713915890595945823395751 * 10000 = 12247.448713915890595945823395751.

Atėmę tikslią ${\sqrt {a}}$ iš 12247.448713915890595945823395751, gauname paklaidą:

12247.448713915890595945823395751 - 12247.448713915890490986420373529 = 1.0495940302222154304017026259951e-13.

O santykinė paklaida skaičiui ${\sqrt {a}}$ yra

(12247.448713915890595945823395751 - 12247.448713915890490986420373529)/12247.448713915890490986420373529 =

= 8.5698993703859497564530032133562e-18.

Tokia pati, kaip skaičiui ${\sqrt {p}}.$

• Apskaičiuosime ${\sqrt {a}}={\sqrt {9.6\cdot 10^{11}}}.$ Čia $p=9.6,\;\;j=1,\;\;v=10.$ Tada pagal (1) formulę

$a=p\cdot 10^{j}\cdot 10^{v}=9.6\cdot 10^{1}\cdot 10^{10}.$

Pagal (2) formulę

${\sqrt {a}}={\sqrt {p}}\cdot {\sqrt {10^{j}}}\cdot 10^{v \over 2}={\sqrt {p}}\cdot {\sqrt {10}}\cdot 10^{10 \over 2}={\sqrt {p}}\cdot {\sqrt {10}}\cdot 10^{5}.$

Šaknį iš 10 lengva apskaičiuoti (nes (4) formulė tinka skaičiams $4\leq a<16$).

${\sqrt {10}}=$ 3.1622776601683793319988935444327.

Toliau randame $x_{1}$ skaičiui $p=9.6$ pagal (4) formulę:

$x_{1}={\frac {1}{6}}p+{\frac {17}{12}}={\frac {1}{6}}\cdot 9.6+{\frac {17}{12}}=$

= 0.16666666666666666666666666666667*9.6 + 1.4166666666666666666666666666667 = 3.0166666666666666666666666666667.

[ ${\sqrt {9.6}}=$ 3.0983866769659335081434123198259. ]

[3.0983866769659335081434123198259 - 3.0166666666666666666666666666667 = 0.0817200102992668414767456531592; tai yra ne santykinė paklaida, o reali paklaida.]

Toliau pagal (3.10) formulę apskaičiuojame $x_{2}$:

$x_{2}=x_{1+1}={\frac {1}{2}}{\Big (}x_{1}+{\frac {p}{x_{1}}}{\Big )}=$

= 0.5*(3.0166666666666666666666666666667 + 9.6/3.0166666666666666666666666666667) = 3.0994935543278084714548802946593.

Toliau pagal (3.10) formulę apskaičiuojame $x_{3}$:

$x_{3}=x_{2+1}={\frac {1}{2}}{\Big (}x_{2}+{\frac {p}{x_{2}}}{\Big )}=$

= 0.5*(3.0994935543278084714548802946593 + 9.6/3.0994935543278084714548802946593) = 3.0983868746074953885623653564026.

[3.0983868746074953885623653564026 - 3.0983866769659335081434123198259 = 0.0000001976415618804189530365767 =

$\approx 1.9764156188\cdot 10^{-7}.$]

Toliau pagal (3.10) formulę apskaičiuojame $x_{4}$:

$x_{4}=x_{3+1}={\frac {1}{2}}{\Big (}x_{3}+{\frac {p}{x_{3}}}{\Big )}=$

= 0.5*(3.0983868746074953885623653564026 + 9.6/3.0983868746074953885623653564026) = 3.0983866769659398117763289613096.

Iš ką tik gautos $x_{4}$ reikšmės atėmę ${\sqrt {9.6}}$ reikšmę, gauname:

3.0983866769659398117763289613096 - 3.0983866769659335081434123198259 = 0.0000000000000063036329166414837 =

$\approx 6.3036329166414837\cdot 10^{-15}.$

Santykinė paklaida yra tokia:

$\varepsilon _{4}={\frac {x_{4}-\gamma }{\gamma }}=$

(3.0983866769659398117763289613096 - 3.0983866769659335081434123198259)/3.0983866769659335081434123198259 =

= 2.0344887755631127106428103904763e-15.

Taigi $\varepsilon _{4}\approx 2.03448877556\cdot 10^{-15}<10^{-10}.$

Iš ankstesnių skaičiavimų

${\sqrt {a}}={\sqrt {p}}\cdot {\sqrt {10}}\cdot 10^{5}=x_{4}\cdot {\sqrt {10}}\cdot 10^{5}=$

= 3.0983866769659398117763289613096 * 3.1622776601683793319988935444327 * 100000 = 979795.89711327323266266864962311.

Tiksli ${\sqrt {a}}$ reikšmė yra tokia:

${\sqrt {a}}={\sqrt {9.6\cdot 10^{11}}}=$ 979795.89711327123927891362988236.

Atėmę tikslią ${\sqrt {a}}$ reikšmę iš 979795.89711327323266266864962311, gauname paklaidą:

979795.89711327323266266864962311 - 979795.89711327123927891362988236 = 0.00000000199338375501974075 =

$\approx 1.993383755\cdot 10^{-9}.$

Santykinė paklaida skaičiui ${\sqrt {a}}$ yra

(979795.89711327323266266864962311 - 979795.89711327123927891362988236)/979795.89711327123927891362988236 =

= 2.0344887755631126970242116177837e-15,

t. y. tokia pati (su 16 pirmų tokių pačių skaitmenų) kaip skaičiui ${\sqrt {p}}.$

[Čia kažkas negerai su tais 16 skaitmenų, nes iš viso ciferblate yra 32 skaitmenys... Jeigu įdėti (padaryti copy-paste) į kalkuliatorių tokį reiškinį:

(3.0983866769659398117763289613096 - 3.0983866769659335081434123198259)/3.0983866769659335081434123198259 -(979795.89711327323266266864962311 - 979795.89711327123927891362988236)/979795.89711327123927891362988236 =

tai gaunamas toks rezultatas:

= 1.3618598772692596082816173047622e-32.

Kažkokie dvigubi skaičiavimo standartai su "Windows 10" kaluliatorium... Gal po kokių nors dalybos iteracijų labai prasiranda tikslumas, nes ne visi skaitmenys Winows 10 kalkuliatoriuje rodomi (jų yra virš 40, gal apie 47 skaitmenys iš viso, iš kurių rodomi tik 32 skaitmenys). Nors gal ir ne tai, nes ten tik atimtis atliekama... Gal koks bug'as?

Bendrai, viskas teisingai, nes

2.0344887755631127106428103904763e-15 - 2.0344887755631126970242116177837e-15 = 1.36185987726926e-32.

Tai rodo, kad skaičiavimo tikslumas yra 32 skaitmenys, nes jų šituose skaičiavimuose tiek ir yra.]

Jeigu $0<a<1,$ tai pažymėję $b={\frac {1}{a}},$ galime, jei b didelis skaičius, apskaičiuoti ${\sqrt {b}}$ pagal pateiktą aukščiau būdą. Apskaičiavę ${\sqrt {b}},$ galime rasti ${\sqrt {a}}$ taip:

${\sqrt {a}}={\frac {1}{\sqrt {b}}}.$

Tačiau toks ${\sqrt {a}}\;$ ($0<a<1$) skaičiavimo būdas nėra ekonomiškas, nes reikia atlikti dvi [papildomas] dalybos operacijas.

Pateiksime ekonomiškesnį, greitesnį ${\sqrt {a}}$ skaičiavimo būdą. Jei $0<a<1,$ tai a galima užrašyti taip:

$a=p\cdot 10^{j}\cdot 10^{v}.\quad (5)$

Čia $1<p<10;\;$ skaičius v gali būti nulis arba neigiamas lyginis skaičius (v gali būti 0, -2, -4, -6, ...); skaičius j lygus arba nuliui arba $-1.$

Jeigu ištrauksime kvadratinę šaknį iš a, tai gausime:

${\sqrt {a}}={\sqrt {p}}\cdot {\sqrt {10^{j}}}\cdot 10^{v \over 2}.\quad (6)$

Jeigu $j=-1,$ tai tereikia apskaičiuoti $({\sqrt {10}})^{-1},$ kurį bus galima visada taikyti (6) formulėje. Vadinasi, tereikia mokėti greitai ištraukti šaknį iš skaičiaus p. O toliau jau lengva pagal (6) formulę surasti kam lygus ${\sqrt {a}}.$ Kaip rasti ${\sqrt {p}}$ pateikta aukščiau.

Kiek pamenu, skaičiau, kad jeigu norima padauginti skaičių dvejetainėje sistemoje iš 2, tai reikia visus bitus pastumti per vieną bito poziciją į kairę. O jeigu norima skaičių dvejetainėje sistemoje padalinti iš 2, tai reikia visus bitus pastumti viena bito pozicija į dešinę. Tokiu budu, jeigu skaičius j nelygus nuliui, tai reikia skaičių a padauginti iš 10, jei 0<a<1 (arba padauginti iš 0.1, jei a>1). Tada skaičius v bus lyginis skaičius. Toliau norint padalinti v iš 2, tereikia rodiklio v visus bitus pastumti per vieną bitą į dešinę. Taip mes gausime v/2. Ir nereikia procesoriaus FPU daryti sudetingų eksponentės dalybos operacijų (skaičiuojant šaknį iš kokio nors skaičiaus).

Lieka tik paslaptis kodėl procesorių FPU šaknies traukimo operacijai sunaudoja tiek pat arba tik apie 1.5 karto (priklausomai nuo procesoriaus) daugiau ciklų negu dalybos operacijai. Net jeigu artinį $x_{1}$ kažkokiu budu galima gauti su santykine paklaida apie $10^{-6},$ tai vis tiek reikės dviejų iteracijų ir to pasekoje reikės dviejų dalybos operacijų (kiekvienoje iteracijoje yra viena dalybos operacija). Kad užtektų apskaičiuoti ${\sqrt {a}}$ iš vienos iteracijos, artinio $x_{1}$ santykinė paklaida $\varepsilon _{1}$ turi būti apie $10^{-9}$ (kad gauti ${\sqrt {a}}$ su santykine paklaida apie $10^{-18}$). Jeigu kažkokiu budu galima pagal skaičių p pažiūrėt į kažkokią artinių $x_{1}$ lentelę, tai ta lentelė turės apie milijardą artinių ir užims apie $80\;({\text{bit}})\cdot 10^{9}=$ 10 GB, t. y. apie 10 gigabaitų. Į procesorių akivaizdu to neįstatysi, tuo labiau, kad, pavyzdžiui, Pentium II/III naudojo 17...38 ciklus dalybos operacijai ir 27...50 ciklų šaknies traukimo operacijai (tik apie 1.5 karto daugiau ciklų šakniai nei dalybai). Tais laikais ir 1 MB (megabaito) kokios nors atminties (pvz., cache atminties) įdėjimas į procesorių buvo techniškai neįmanomas ar sunkiai pasiekiamas dalykas.

Kadangi artinio $x_{1}$ skaičiavimui [pateiktam aukščiau] nereikia daryti dalybos operacijos, tai galbūt yra, tarkim, apie 100 ar 1000 artinio $x_{1}$ skaičiavimo grupių skaičiui ${\sqrt {p}}.$ Kaip matėme, tik dviejų grupių (1<p<4 ir $4\leq p<10$) užtenka ${\sqrt {p}}$ skaičiavimui, kad gauti pirmo artinio $x_{1}$ santykinę paklaidą $\varepsilon _{1}<1/24<0.05.$ Jeigu sukurti apie 10 tokių artinio $x_{1}$ grupių (1<p<2, 2<p<3, 3<p<4, ..., 9<p<10) arba apie 100 tokių artinio $x_{1}$ grupių (1<p<1.1, 1.1<p<1.2, 1.2<p<1.3, ..., 9.9<p<10), tai peršasi mintis, kad 10-iai grupių tikslumas gali būti apie $10^{-5},$ o šimtui tokių grupių santykinės paklaidos $\varepsilon _{1}$ tikslumas gali būti apie $10^{-10}.$ Tik tai gali užimti darbo suradimui artinio $x_{1}$ formulės kiekvienai grupei (tokia formulė padidinanti tikslumą siauro/trumpo intervalo grupei turėtų egzistuot matematiškai mąstant ir lažinantis...). Tik gal tikslumo padidėjimas gali būti ne toks žymus, kai grupių padaugėja apie 10 kartų (arba 100 kartų), nežinau...

Kitas klausimas yra, kaip kiekvieną skaičių p nusiusti į jį atitinkančią vieną iš 10 (ar vieną iš 100) grupę (kad paskui gauti artinį $x_{1}$). Programuotojas sakytų, kad reikia vidutiniškai praifint 5 kartus, jei yra 10 grupių (arba ~50 kartų, jei yra 100 grupių). Pavyzdžiui, jei yra skaičius 3.5346, tai pirmą kartą paklaust ar skaičius p [jei yra 10 grupių] yra mažiau už 2, antrą kartą paklaust ar skaičius p yra mažiau už 3 ir trečią kartą palyginus 3.5346 su 4 nustatoma, kad skaičius p patenka į trečią grupę, nes yra mažiau už 4. Bet su ifinimais šaknies skaičiavimai bus dar ilgesni. Todėl su tranzistoriais galima sukurt painią schemą, kuri, jeigu skaičiaus p pirmas skaitmuo yra 3, siunčia skaičių p į trečią grupę, kurioje skaičiuoja artinį $x_{1}$ skaičiams nuo 3 iki 4 ir gauna numanomai minėtą tikslumą $10^{-5}.$ O jeigu skaičiaus p pirmas skaitmuo yra 4, tai siunčiasi per tą loginę schemą toks skaičius į ketvirtą grupę (į kurią įeina skaičiai nuo 4 iki 5) ir t. t. Kas ten žino, gal taip ir randamas [skaičiaus p] $x_{1}$ artinys (apie $10^{-10}$ tikslumu) naudojant apie 100 tokių grupių. Ir tada užtenka tik vienos iteracijos (t. y. užtenka vienos dalybos operacijos [kuri yra iteracijoje]), kad gauti apie $10^{-18}$ tikslumą.

Pastebėsime, kad jeigu artinys $x_{1}={\sqrt {a}},$ tai pagal (3.10) formulę

$x_{2}={\frac {1}{2}}{\Big (}x_{1}+{\frac {a}{x_{1}}}{\Big )}={\frac {1}{2}}{\Big (}{\sqrt {a}}+{\frac {a}{\sqrt {a}}}{\Big )}={\frac {1}{2}}{\Big (}{\sqrt {a}}+{\sqrt {a}}{\Big )}={\sqrt {a}}.$

Čia aiškinama kaip dalybos operacija vyksta procesoriuje:

https://www.quora.com/How-is-integer-division-implemented-in-Computer-Hardware

https://en.wikipedia.org/wiki/Division_algorithm#SRT_division

Pavyzdžiui čia:

https://electronics.stackexchange.com/questions/22410/how-does-division-occur-in-our-computers

skaitant, gali būti neaišku kas yra "reciprocal", "denominator", "numerator". Jie yra atitinkamai 1/x, vardiklis, skaitiklis.

Šis tekstas:

"For floating point numbers, however, it may be quicker to first compute the reciprocal of the denominator and then multiply that times your numerator."

verčiamas per Google translate taip:

"Tačiau slankiojo kablelio skaičiams gali būti greičiau pirmiausia apskaičiuoti vardiklio atvirkštinį skaičių ir tada padauginti jį iš skaitiklio."

Ten pat duotas pavyzdys, jeigu norima apskaičiuoti 1/3.5. Iš pradžių spėjama, koks turi būti g, kad g*3.5 būtų lygus 1. Spėjimas, kad g=0.3. Tada g*3.5=0.3*3.5=1.05. Tada toliau naudojama pagerinta formulė g'=g(2-gd). Gauname 0.3*(2-1.05)=0.285. Jau gana arti 1/3.5=0.28571428571428571428571428571429. Toliau kartojamas procesas: 0.285*(2-0.285*3.5)=0.2857125. Dar toliau: 0.2857125*(2-0.2857125*3.5)=0.285714285703125. Paklaida yra:

0.28571428571428571428571428571429 - 0.285714285703125 = 0.00000000001116071428571428571429 $\approx 1.11607142857\cdot 10^{-11}.$

Čia

https://www.reddit.com/r/explainlikeimfive/comments/yc5qkp/eli5_how_does_a_computer_do_division_efficiently/

dar papildomai rašo apie šitą būdą, kuris ką tik buvo pateiktas skaičiuot 1/D, o suradus kam tai lygu, galima paskui sužinot kam lygu N/D, padauginus N iš apskaičiuoto 1/D (kai D=p, t. y. kai 1<D<10, o jei D>10, tai skaičiaus D eksponentės rodiklis atimamas iš skaičiaus N eksponentės rodiklio...).

Šito metodo išvedimas aprašytas angliškoje Vikipedijoje: https://en.wikipedia.org/wiki/Division_algorithm#Newton–Raphson_division

Jame naudojamas Niutono metodas (liesinių metodas), kuris taikomas polinomo realiųjų šaknų apytikslėms reikšmėms rasti.

Greitas šaknies iš skaičiaus a traukimo būdas

Parodysime praktiškai greičiausią kvadratinės šaknies traukimo būdą. Šitas būdas paremtas Niutono metodu (liestinių metodu). Todėl pirma trumpai apžvelgsime Niutono metodą.

Niutono metodas. Naudojant Niutono metodą galima rasti polinomo (arba kokios nors funkcjos) f(x) šaknį (sprendinį) bet kokiu norimu tikslumu. Nagrinėkime atkarpą [a; b], kurioje yra viena šaknis c (taške c funkcija f(x) kerta ašį Ox). Toliau, priklausomai nuo to ar funkcija f(x) atkarpoje [a; b] įgaubta ar išgaubta bei priklausomai nuo to ar funkcija f(x) didėja ar mažėja atkarpoje [a; b], pasirenkamas tinkamas atkarpos vienas iš galų (a arba b). Pažymėkime tą galą $a_{0}$ ($a_{0}=a$ arba $a_{0}=b$). Toliau taške $a_{0}$ surandame funkcijos f(x) liestinės lygtį, kuri atrodo taip:

$y=f(a_{0})+f'(a_{0})(x-a_{0}).$

Arba taip:

$y-f(a_{0})=f'(a_{0})(x-a_{0}).\quad (1)$

Toliau randamas liestinės y(x) susikirtimo taškas su Ox ašimi. Tą tašką pažymėkime (x; y)=(d; f(d))=(d; 0). Įstatę šio taško koordinates į (1) lygtį gausime:

$0-f(a_{0})=f'(a_{0})(d-a_{0}),$

$-{\frac {f(a_{0})}{f'(a_{0})}}=d-a_{0},$

$d=a_{0}-{\frac {f(a_{0})}{f'(a_{0})}}.\quad (2)$

Dabar turime siauresnį intervalą, kuriame ant Ox ašies yra šaknis c (per kurią pereina funkcjos f(x) grafikas). Ir dabar reikšmė $d=a_{1}$ yra arčiau šaknies c. Pakartoję šį procesą rasime tašką $a_{2},$ esantį dar arčiau prie c taško (ant Ox ašies):

$a_{2}=a_{1}-{\frac {f(a_{1})}{f'(a_{1})}}.$

Taip kartodami (2) formulę galime labai tiksliai rasti funkcijos f(x) susikirtimo su Ox ašimi tašką c. Šios formulės paklaida mažėja kvadratu po kiekvienos iteracijos, kai $a_{n}$ konverguoja į tašką c (po kiekvienos iteracijos, tikslumas didėja daug greičiau, nei jeigu spėliotume šaknį c, dalindami intervalą [a; b] kaskart per pusę ir tikrindami ar skirtinguose intervalų pusėse funkcija f(x) turi priešingus ženklus). Išsamiai apie Niutono metodą galima paskaityti aukštosios algebros arba matematinės analizės vadovėliuose.

Čia pasakojama apie ${\sqrt {a}}$ skaičiavimo metodus: https://www.quora.com/How-do-computers-calculate-square-roots

O čia

https://en.wikipedia.org/wiki/Fast_inverse_square_root#Newton's_method

aprašoma, kaip apskaičiuoti $1/{\sqrt {a}}.$ Kai apytiksliai apskaičiuojama $1/{\sqrt {a}},$ reikia gautą rezultatą padauginti iš a. Tada ir gausime ${\sqrt {a}}.$ Be to $1/{\sqrt {a}}$ rezultatas yra net naudingesnis, nes nereikia dalinti ${\sqrt {a}}$ iš a, kad gauti $1/{\sqrt {a}}.$

Pereikime prie ${\frac {1}{\sqrt {a}}}$ skaičiavimo algoritmo, kuris vadinasi Niutono metodu [šaknies skaičiavimui]. Pažymėkime $y={\frac {1}{\sqrt {a}}},$ o skaičių a pažymėkime per x. Tada

$f(y)={\frac {1}{y^{2}}}-x=0.$

Jeigu $y_{n}$ yra apytikslė skaičiaus y reikšmė, tai $y_{n+1}$ yra dar tikslesnė skaičiaus y reikšmė (pagal liestinių metodą):

$y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}.\quad (3)$

Surandame f(y) išvestinę:

$f'(y)=({\frac {1}{y^{2}}}-x)'={\frac {-2}{y^{3}}}.$

$f'(y_{n})={\frac {-2}{y_{n}^{3}}}.$

Toliau įstatome $f'(y_{n})$ ir $f(y_{n})$ į (3) formulę:

$y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}=y_{n}-{\frac {{\frac {1}{y_{n}^{2}}}-x}{\frac {-2}{y_{n}^{3}}}}=y_{n}+{\frac {{\frac {y_{n}^{3}}{y_{n}^{2}}}-xy_{n}^{3}}{2}}=$

$=y_{n}+{\frac {y_{n}-xy_{n}^{3}}{2}}=y_{n}+{\frac {y_{n}}{2}}-{\frac {xy_{n}^{3}}{2}}={\frac {3y_{n}}{2}}-{\frac {xy_{n}^{3}}{2}}=$

$=y_{n}\left({\frac {3}{2}}-{\frac {x}{2}}y_{n}^{2}\right)={\frac {y_{n}\left(3-xy_{n}^{2}\right)}{2}}.$

Gavome rekurentinę formulę $1/{\sqrt {a}}$ reikšmei skaičiuoti:

$y_{n+1}=y_{n}\left({\frac {3}{2}}-{\frac {x}{2}}y_{n}^{2}\right)={\frac {y_{n}\left(3-xy_{n}^{2}\right)}{2}}.\quad (4)$

Pavyzdžiai

• Apskaičiuosime ${\frac {1}{\sqrt {x}}}={\frac {1}{\sqrt {2.3}}}$ ir rasime kam lygus ${\sqrt {2.3}}={\sqrt {a}}={\frac {x}{\sqrt {x}}}={\frac {2.3}{\sqrt {2.3}}}.$

Pirmiausia turime žinoti $y={\frac {1}{\sqrt {2.3}}}$ apytikslią reikšmę. Tą reikšmę procesoriaus FPU gali gauti iš tam tikros lentelės, kuri galbūt yra procesoriaus L1 Data cache'e arba L2 cache. Tiksli ${\sqrt {2.3}}$ reikšmė yra tokia:

${\sqrt {2.3}}=$ 1.5165750888103101108513650872564.

O $y={\frac {1}{\sqrt {2.3}}}=$ 0.65938047339578700471798482054627.

Mes paimsime, kad $y_{1}=0.6.$ Tada pagal formulę

$y_{n+1}=y_{n}\left({\frac {3}{2}}-{\frac {x}{2}}y_{n}^{2}\right)={\frac {y_{n}\left(3-xy_{n}^{2}\right)}{2}}$

rasime

$y_{2}={\frac {y_{1}\left(3-xy_{1}^{2}\right)}{2}}={\frac {0.6\left(3-2.3\cdot 0.6^{2}\right)}{2}}={\frac {0.6\left(3-0.828\right)}{2}}={\frac {0.6\cdot 2.172}{2}}={\frac {1.3032}{2}}=0.6516.$

Toliau pagal (4) formulę rasime $y_{3}$:

$y_{3}={\frac {y_{2}\left(3-xy_{2}^{2}\right)}{2}}={\frac {0.6516\left(3-2.3\cdot 0.6516^{2}\right)}{2}}={\frac {0.6516\left(3-2.3\cdot 0.42458256\right)}{2}}={\frac {0.6\cdot 2.023460112}{2}}={\frac {1.3184866089792}{2}}=0.6592433044896.$

[0.65938047339578700471798482054627 - 0.6592433044896 = 0.00013716890618700471798482054627 $\approx 1.37168906187\cdot 10^{-4}.$]

Toliau pagal (4) formulę rasime $y_{4}$:

$y_{4}={\frac {y_{3}\left(3-xy_{3}^{2}\right)}{2}}=$

= 0.6592433044896*(3 - 2.3*0.6592433044896^2)/2 = 0.65938043059651203399501910136126.

[0.65938047339578700471798482054627 - 0.65938043059651203399501910136126 =

= 0.00000004279927497072296571918501 $\approx 4.279927497\cdot 10^{-8}.$]

Toliau pagal (4) formulę rasime $y_{5}$:

$y_{5}={\frac {y_{4}\left(3-xy_{4}^{2}\right)}{2}}=$

= 0.65938043059651203399501910136126*(3 - 2.3*0.65938043059651203399501910136126^2)/2 = 0.65938047339578283767489142997287.

Taigi, $y_{5}$ paklaida yra tokia:

0.65938047339578700471798482054627 - 0.65938047339578283767489142997287 =

= 0.0000000000000041670430933905734 $=4.1670430933905734\cdot 10^{-15}.$

Po 4 iteracijų gavome paklaidą mažesne nei $10^{-14}.$

Dabar galime rasti a šaknyje:

${\sqrt {2.3}}={\sqrt {a}}={\sqrt {x}}=2.3\cdot y_{5}=$ 2.3*0.65938047339578283767489142997287 = 1.5165750888103005266522502889376.

[1.5165750888103101108513650872564 - 1.5165750888103005266522502889376 = 0.0000000000000095841991147983188 =

$=9.5841991147983188\cdot 10^{-15}.$]

Pažiūrėkime kaip konverguoja $y_{n+1}$ į ${\frac {1}{\sqrt {2.3}}}$ iš dešinės pusės, kai $y_{1}$ parenkamas taip, kad $y_{1}>{\frac {1}{\sqrt {2.3}}}.$

Taigi parinkime $y_{1}=0.66.$ Tada

$y_{2}={\frac {y_{1}\left(3-xy_{1}^{2}\right)}{2}}={\frac {0.66\left(3-2.3\cdot 0.66^{2}\right)}{2}}={\frac {0.66\left(3-2.3\cdot 0.4356\right)}{2}}={\frac {0.66\left(3-1.00188\right)}{2}}={\frac {0.66\cdot 1.99812}{2}}={\frac {1.3187592}{2}}=0.6593796.$

[0.65938047339578700471798482054627 - 0.6593796 = 0.00000087339578700471798482054627 $\approx 8.73395787\cdot 10^{-7}.$]

Toliau pagal (4) formulę rasime $y_{3}$:

$y_{3}={\frac {y_{2}\left(3-xy_{2}^{2}\right)}{2}}={\frac {0.6593796\left(3-2.3\cdot 0.6593796^{2}\right)}{2}}={\frac {0.6593796\left(3-0.999997350861168\right)}{2}}=$

$={\frac {0.6593796\cdot 2.000002649138832}{2}}=$ 1.3187609467881033886272/2 = 0.6593804733940516943136.

[0.65938047339578700471798482054627 - 0.6593804733940516943136 = 0.00000000000173531040438482054627 $\approx 1.7353104\cdot 10^{-12}.$]

Toliau pagal (4) formulę rasime $y_{4}$:

$y_{4}={\frac {y_{3}\left(3-xy_{3}^{2}\right)}{2}}=$

0.6593804733940516943136*(3 - 2.3*0.6593804733940516943136^2)/2 = 0.65938047339578700471797797024742.

Paklaida lygi

0.65938047339578700471798482054627 - 0.65938047339578700471797797024742 =

= 0.00000000000000000000000685029885 $=6.85029885\cdot 10^{-24}.$

Po 3 iteracijų gavome paklaidą mažesne nei $10^{-23}.$

Šį kartą, parenkant $y_{1}$ iš dešinės, $y_{1}$ reikšmė buvo daug tikslesnė negu anksčiau. Bet teoriškai Niutono metodas (liestinių metodas) bent jau polinomams gali greitai konverguoti tik iš vienos šaknies y pusės. Čia gi greitas konvergavimas į ${\frac {1}{\sqrt {2.3}}}$ vyksta tiek iš kairės pusės, tiek iš dešinės pusės, kas atrodo kaip prieštaravimas liestinių metodo teorijai. Gal greito konvergavimo tik iš vienos pusės taisyklė taikoma tik polinomams ir kai kurioms funkcjoms, o ši funkcija $f(y)={\frac {1}{y^{2}}}-x$ yra išimtis.

Dabar galime rasti a šaknyje:

${\sqrt {2.3}}={\sqrt {a}}={\sqrt {x}}=x\cdot y_{4}=2.3\cdot y_{4}=$ 2.3*0.65938047339578700471797797024742 = 1.5165750888103101108513493315691.

[1.5165750888103101108513650872564 - 1.5165750888103101108513493315691 = 0.0000000000000000000000157556873 =

$=1.57556873\cdot 10^{-23}.$]

• Apskaičiuosime ${\frac {1}{\sqrt {x}}}={\frac {1}{\sqrt {5.8}}}$ ir rasime kam lygus ${\sqrt {5.8}}={\sqrt {x}}.$

Tiksli ${\sqrt {5.8}}$ reikšmė yra tokia:

${\sqrt {5.8}}=$ 2.4083189157584590960256482060757.

O $y={\frac {1}{\sqrt {5.8}}}=$ 0.41522739926869984414235313897857.

Mes paimsime, kad $y_{1}=0.4.$ Tada pagal (4) formulę

$y_{2}={\frac {y_{1}\left(3-xy_{1}^{2}\right)}{2}}={\frac {0.4\left(3-5.8\cdot 0.4^{2}\right)}{2}}={\frac {0.4\left(3-0.928\right)}{2}}={\frac {0.4\cdot 2.072}{2}}={\frac {0.8288}{2}}=0.4144.$

[0.41522739926869984414235313897857 - 0.4144 = 0.00082739926869984414235313897857 $\approx 8.27399268699844\cdot 10^{-4}.$]

Toliau pagal (4) formulę rasime $y_{3}$:

$y_{3}={\frac {y_{2}\left(3-xy_{2}^{2}\right)}{2}}=$

= 0.4144*(3 - 5.8*0.4144^2)/2 = 0.4152249278464.

[0.41522739926869984414235313897857 - 0.4152249278464 = 0.00000247142229984414235313897857 $\approx 2.471422299844\cdot 10^{-6}.$]

Toliau pagal (4) formulę rasime $y_{4}$:

$y_{4}={\frac {y_{3}\left(3-xy_{3}^{2}\right)}{2}}=$

= 0.4152249278464*(3 - 5.8*0.4152249278464^2)/2 = 0.41522739924663512944560082817282.

[0.41522739926869984414235313897857 - 0.41522739924663512944560082817282 = 0.00000000002206471469675231080575 =

$\approx 2.206471469675231\cdot 10^{-11}.$]

Toliau pagal (4) formulę rasime $y_{5}$:

$y_{5}={\frac {y_{4}\left(3-xy_{4}^{2}\right)}{2}}=$

= 0.41522739924663512944560082817282*(3 - 5.8*0.41522739924663512944560082817282^2)/2 = 0.41522739926869984414059439797726.

[0.41522739926869984414235313897857 - 0.41522739926869984414059439797726 = 0.00000000000000000000175874100131 =

$=1.75874100131\cdot 10^{-21}.$]

Tada ${\sqrt {5.8}}={\sqrt {x}}=5.8\cdot y_{5}=$ 5.8*0.41522739926869984414059439797726 = 2.4083189157584590960154475082681.

[2.4083189157584590960256482060757 - 2.4083189157584590960154475082681 = 0.0000000000000000000102006978076 =

$=1.02006978076\cdot 10^{-20}.$]

Po 4 iteracijų gavome paklaidą mažesne nei $10^{-19}$ (beveik $10^{-20}$).

Dabar pažūrėsime kaip konverguos $y_{n}$ į $1/{\sqrt {5.8}}$ iš dešinės pusės (kai $y_{1}>1/{\sqrt {5.8}}$). Parenkame $y_{1}=0.42.$ Tada pagal (4) formulę

$y_{2}={\frac {y_{1}\left(3-xy_{1}^{2}\right)}{2}}={\frac {0.42\left(3-5.8\cdot 0.42^{2}\right)}{2}}={\frac {0.42\left(3-1.02312\right)}{2}}={\frac {0.42\cdot 1.97688}{2}}=0.4151448.$

[0.41522739926869984414235313897857 - 0.4151448 = 0.00008259926869984414235313897857 $\approx 8.2599268699844\cdot 10^{-5}.$]

Toliau pagal (4) formulę rasime $y_{3}$:

$y_{3}={\frac {y_{2}\left(3-xy_{2}^{2}\right)}{2}}=$

= 0.4151448*(3 - 5.8*0.4151448^2)/2 = 0.4152273746236976005632.

[0.41522739926869984414235313897857 - 0.4152273746236976005632 = 0.00000002464500224357915313897857 $\approx 2.4645002243579\cdot 10^{-8}.$]

Toliau pagal (4) formulę rasime $y_{4}$:

$y_{4}={\frac {y_{3}\left(3-xy_{3}^{2}\right)}{2}}=$

= 0.4152273746236976005632*(3 - 5.8*0.4152273746236976005632^2)/2 = 0.41522739926869765000924208032344.

[0.41522739926869984414235313897857 - 0.41522739926869765000924208032344 = 0.00000000000000219413311105865513 =

$\approx 2.194133111\cdot 10^{-15}.$]

Toliau pagal (4) formulę rasime $y_{5}$:

$y_{5}={\frac {y_{4}\left(3-xy_{4}^{2}\right)}{2}}=$

= 0.41522739926869765000924208032344*(3 - 5.8*0.41522739926869765000924208032344^2)/2 = 0.41522739926869984414235313896118.

[0.41522739926869984414235313897857 - 0.41522739926869984414235313896118 = 0.00000000000000000000000000001739 =

$=1.739\cdot 10^{-29}.$]

Tada ${\sqrt {5.8}}={\sqrt {x}}=5.8\cdot y_{5}=$ 5.8*0.41522739926869984414235313896118 = 2.4083189157584590960256482059748.

[2.4083189157584590960256482060757 - 2.4083189157584590960256482059748 = 0.0000000000000000000000000001009 =

$=1.009\cdot 10^{-28}.$]

Po 4 iteracijų gavome ${\sqrt {5.8}}$ paklaidą mažesne nei $10^{-27}$ (beveik $10^{-28}$).

[2.4083189157584590960256482059748^2 = 5.7999999999999999999999999995139.]

Pastebėsime, kad abiejuose pavyzdžiuose (skaičiuojant ${\sqrt {2.3}}$ ir ${\sqrt {5.8}}$) skaičiuojant $y_{n}$ iš dešinės (kai $y_{1}>1/{\sqrt {2.3}}$ arba $y_{1}>1/{\sqrt {5.8}}$), jau po pirmos iteracijos $y_{n}<1/{\sqrt {x}}$ (t. y. jau $y_{2}<1/{\sqrt {2.3}}$ arba $y_{2}<1/{\sqrt {5.8}}$). Vadinasi atlikus pirmą iteraciją $y_{n}$ skaičiavimai toliau vyksta iš kairės (pagal Niutono metodą).

Štolco teorema

• Štolco teorema

Trisdešimt antro laipsnio šaknies traukimas taikant Niutono metodą

Skaičiuojant ln(z), kai 0.7<z<1.4, reikalingas greitas šešiolikto ar trisdešimt antro laipsnio šaknies traukimo algoritmas (apie tokius natūrinio logaritmo skaičiavimus parašyta čia https://lt.wikibooks.org/wiki/Matematika/Teiloro_eilutė#Natūrinio_logaritmo_liekamojo_nario_įvertinimas ).

Parodysime kaip apskaičiuoti ${\frac {1}{\sqrt[{32}]{z}}}={\frac {1}{z^{1/32}}}$ taikant Niutono metodą. Pažymėkime $y={\frac {1}{z^{1/32}}}.$ Tada turime funkciją

$f(y)={\frac {1}{y^{32}}}-z=0.$

Pagal liestinių metodą

$y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}.\quad (3)$

Surandame f(y) išvestinę:

$f'(y)=({\frac {1}{y^{32}}}-z)'={\frac {-32}{y^{33}}}.$

$f'(y_{n})={\frac {-32}{y_{n}^{33}}}.$

Toliau įstatome $f'(y_{n})$ ir $f(y_{n})$ į (3) formulę:

$y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}=y_{n}-{\frac {{\frac {1}{y_{n}^{32}}}-z}{\frac {-32}{y_{n}^{33}}}}=y_{n}+{\frac {{\frac {y_{n}^{33}}{y_{n}^{32}}}-zy_{n}^{33}}{32}}=$

$=y_{n}+{\frac {y_{n}-zy_{n}^{33}}{32}}=y_{n}+{\frac {y_{n}}{32}}-{\frac {zy_{n}^{33}}{32}}={\frac {33y_{n}}{32}}-{\frac {zy_{n}^{33}}{32}}=$

$=y_{n}\left({\frac {33}{32}}-{\frac {z}{32}}y_{n}^{32}\right)={\frac {y_{n}\left(33-zy_{n}^{32}\right)}{32}}.$

Gavome rekurentinę formulę ${\frac {1}{\sqrt[{32}]{z}}}$ reikšmei skaičiuoti:

$y_{n+1}={\frac {y_{n}\left(33-zy_{n}^{32}\right)}{32}}.\quad (5)$

Pavyzdžiai

• Apskaičiuosime ${\frac {1}{\sqrt[{32}]{z}}}={\frac {1}{\sqrt[{32}]{1.3}}}$ ir rasime kam lygus ${\sqrt[{32}]{1.3}}={\sqrt[{32}]{z}}.$ Prieš tai dar pastebėsime, kad norint apskaičiuoti ${\sqrt[{32}]{z}}$ reikia 5 kartus ištraukti kvadratinę šaknį iš z. Tiksli ${\sqrt[{32}]{1.3}}$ reikšmė yra tokia:

${\sqrt[{32}]{1.3}}=$ 1.0082325861537202021533989768813.

O $y={\frac {1}{\sqrt[{32}]{1.3}}}=$ 0.99183463590963025992221676151568.

Mes paimsime, kad $y_{1}=0.99.$ Tada pagal (5) formulę

$y_{2}={\frac {y_{1}\left(33-zy_{1}^{32}\right)}{32}}={\frac {0.99\left(33-1.3\cdot 0.99^{32}\right)}{32}}=$

= 0.99*(33 - 1.3*0.72498033595785364231768779197961)/32 = 0.99177969711319507382303549411632.

[0.99183463590963025992221676151568 - 0.99177969711319507382303549411632 = 0.00005493879643518609918126739936 $\approx 5.49\cdot 10^{-5}.$]

Toliau pagal (5) formulę rasime $y_{3}$:

$y_{3}={\frac {y_{2}\left(33-zy_{2}^{32}\right)}{32}}=$

= 0.99177969711319507382303549411632*(33 - 1.3*0.99177969711319507382303549411632^32)/32 = 0.99183458572688580323137637227142.

[0.99183463590963025992221676151568-0.99183458572688580323137637227142=0.00000005018274445669084038924426 $\approx 5.018\cdot 10^{-8}.$]

Toliau pagal (5) formulę rasime $y_{4}$:

$y_{4}={\frac {y_{3}\left(33-zy_{3}^{32}\right)}{32}}=$

= 0.99183458572688580323137637227142*(33 - 1.3*0.99183458572688580323137637227142^32)/32 = 0.99183463590958836578366677638205.

(0.99183458572688580323137637227142*(33 - 1.3*0.76922952379474351345468696687987)/32 = 0.99183463590958836578366677638205.)

[0.99183463590963025992221676151568 - 0.99183463590958836578366677638205 = 0.00000000000004189413854998513363 =

$\approx 4.1894138549985\cdot 10^{-14}.$]

Toliau pagal (5) formulę rasime $y_{5}$:

$y_{5}={\frac {y_{5}\left(33-zy_{5}^{32}\right)}{32}}=$

= 0.99183463590958836578366677638205*(33 - 1.3*0.99183463590958836578366677638205^32)/32 = 0.9918346359096302599222167323178.

Atėmę, ką tik gautą $y_{5}$ reikšmę iš tikslios $y={\frac {1}{\sqrt[{32}]{1.3}}}$ reikšmės, gauname:

0.99183463590963025992221676151568 - 0.9918346359096302599222167323178 = 0.00000000000000000000000002919788 =

$\approx 2.919788\cdot 10^{-26}.$

Gavome po keturių iteracijų paklaidą mažesnę už $10^{-25}.$

Padalijus vienetą iš $y_{5}$ reikšmės gauname

${\sqrt[{32}]{1.3}}={\sqrt[{32}]{z}}={\frac {1}{y_{5}}}=$

= 1/0.9918346359096302599222167323178 = 1.0082325861537202021533990065619.

[1.0082325861537202021533989768813 - 1.0082325861537202021533990065619 = -0.0000000000000000000000000296806 =

$\approx -2.96806\cdot 10^{-26}.$]

1.0082325861537202021533990065619^32 = 1.3000000000000000000000012246312.

(kad pakelti 1.0082325861537202021533990065619 32 laipsniu reikia ant kalkuliatoriaus 5 kartus paspasuti mygtuką $x^{2}.$)

Toliau apskaičiuosime ln(1.3). Kadangi

$\ln 1.3=-\ln 1.3^{-1}=-32\ln 1.3^{-1/32}=-32\ln {\frac {1}{z^{1/32}}}\approx -32\ln y_{5},$

tai tereikia apskaičiuoti $\ln y_{5}.$ Gauname

$\ln y_{5}=$ ln(0.9918346359096302599222167323178) = ln(1-0.0081653640903697400777832676822), x=-0.0081653640903697400777832676822.

Pagal Teiloro formulę ( https://lt.wikibooks.org/wiki/Matematika/Teiloro_eilutė#Natūrinio_logaritmo_liekamojo_nario_įvertinimas ) apskaičiuosime -32*ln(1+x) = -32*ln(1-0.0081653640903697400777832676822). Imsime n=9. Tada

$\ln(1+x)=x-{\frac {x^{2}}{2}}+{\frac {x^{3}}{3}}-{\frac {x^{4}}{4}}...+(-1)^{n-1}{\frac {x^{n}}{n}}+r_{n+1}(x);\quad (13)$

čia $|-32r_{n+1}(x)|={\Big |}-32{\frac {x^{n+1}}{(1+c)^{n+1}(n+1)}}{\Big |}\;,\;\;-0.008165364<c<0.$

Vietoje x įstatome -0.0081653640903697400777832676822 į (13) formulę ir su n=9, gauname

$\ln(1+x)=\ln(1-0.0081653640903697400777832676822)=$

= -0.0081653640903697400777832676822 - (-0.0081653640903697400777832676822)^2/2 + (-0.0081653640903697400777832676822)^3/3 - (-0.0081653640903697400777832676822)^4/4 + (-0.0081653640903697400777832676822)^5/5 - (-0.0081653640903697400777832676822)^6/6 + (-0.0081653640903697400777832676822)^7/7 - (-0.0081653640903697400777832676822)^8/8 + (-0.0081653640903697400777832676822)^9/9 =

(galima nuo skaičiaus "-0.0081653640903697400777832676822" iki ženklo "=" nukopijuoti visą eilutę su 9 dėmenim ir įstatyti (padaryti "Paste") į Windows 10 kalkuliatorių ir iškart bus gautas apskaičiuotas atsakymas)

= -0.00819888326460909537597654216045.

[-0.00819888326460909537597654216045 - ln(0.9918346359096302599222167323178) = 1.3273686782585303927157997035195e-22.]

Gavome paklaidą $\approx 1.327368678\cdot 10^{-22}<10^{-21}.$

Įvertinsime liekamąjį narį teoriškai.

$|-32r_{9+1}(x)|<{\Big |}-32\cdot {\frac {(-0.008165364)^{9+1}}{(1-0.008165364)^{9+1}(9+1)}}{\Big |}=$

=| -32*((-0.0081653640903697400777832676822)^10)/(10*(1-0.0081653640903697400777832676822)^10) | =

=| -32*((-0.0081653640903697400777832676822)^10)/(10*0.9918346359096302599222167323178^10) | =

(čia galima padaryti "Copy-Paste", viską nukopijavus, išskyrus modulio ženklus "|" ir paspaudus paskui kalkuliatoriuje ženklą "="; tokiu budu bus gautas žemiau esantis atsakymas)

=| -4.5762824470151633013165981309199e-21 | $<10^{-20}.$

Padauginsime anksčiau per Teiloro eilutę gautą atsakymą -0.00819888326460909537597654216045 iš -32, kad gauti ln(1.3). Taigi,

$\ln(1.3)\approx$ -32*(-0.00819888326460909537597654216045) = 0.2623642644674910520312493491344.

Atėmę apskaičiuotą ln(1.3) [apytikslią] reikšmę iš Windows 10 kalkuliatoriaus tikslios reikšmės, gauname paklaidą:

ln(1.3) - 0.2623642644674910520312493491344 = 0.26236426446749105203549598688095 - 0.2623642644674910520312493491344 =

= 4.2466377465543972041664561314341e-21.

Taigi, gavome paklaidą $4.246637746554\cdot 10^{-21}<|-4.576282447\cdot 10^{-21}|<10^{-20}.$

Ir čia viskas atitinka teoriją.

Palyginsime kiek reikia daugybos ir sudeties operacijų, 5 kartus traukiant kvadratinę šaknį su vieno karto trisdešimt antro laipsnio šaknies traukimo formulę. Priminsime 32-o laipsnio šaknies traukimo formulę (tiksliau 1/z^(1/32))

$y_{n+1}={\frac {y_{n}\left(33-zy_{n}^{32}\right)}{32}}.\quad (5)$

Kad $y_{n}$ pakelti 32 laipsniu, reikia tiesiog $y_{n}$ penkis kartus pakelti kvadratu. Tai jau yra 5 daugybos operacijos. Toliau dar reikia $y_{n}$ padauginti iš z. Jau turime 6 daugybos operacijas. Tada reikia $zy_{n}$ atimti iš 33. Gauname vieną atimties operaciją. Ir galiausiai viską padauginti iš $y_{n}$ ir iš 1/32. Gavosi iš viso 6+2=8 daugybos operacijos ir viena atimties operacija. Truputi tiksliau [nei mes] iš procesoriaus FPU lentelės pasirinkus $y_{1},$ reikšmę užteks 3 iteracijų, kad gauti 18 skaitmenų tikslumą. Todėl kad ištraukti 32-o laipsnio šaknį reikia 8*3=24 daugybos operacijos ir 1*3=3 atimties operacijų. Iš viso 24+3=27 operacijos (nors šiuolaikiniai procesoriai turi FMA (fuse multiply add) instrukcija, kuri leidžia per vieną operaciją sudauginti du skaičius ir pridėti trečią skaičių). Daugybos ir sudeties operacijos x87 FPU skaičiuojamos panašiu greičiu. Čia nereikia vienos dalybos iteracijos, kad atlikti $1/y_{4},$ nes skaičiuojant ln(1+x) be to apseinama kaip galima matyti iš aukščiau.

Išnagrinėsime kiek reikia daugybos ir sudeties operacijų, jei penkis kartus traukti kvadratinę šaknį iš skaičiaus z. Kvadtatinės šaknies (tiksliau 1/x^(1/2)) traukimo vienos iteracijos formulė yra tokia:

$y_{n+1}={\frac {y_{n}\left(3-xy_{n}^{2}\right)}{2}}.\quad (4)$

Kad pakelti $y_{n}$ kvadratu reikia vienos daugybos operacijos. Paskui dar vienos daugybos operacijos reikia, kad padauginti $y_{n}^{2}$ iš x. Jau turime dvi daugybos operacijas. Čia yra dar viena atimties operacija. Dar viena daugybos operacija reiškinio skliaustuose iš $y_{n}$ ir dar viena daugybos operacija visko padauginimas iš 0.5. Taigi, iš viso reikia 2+2=4 daugybos operacijų ir vienos atimties operacijos. Procesoriui gana tiksliai pasirinkus pirmą artinį $y_{1},$ užteks 3 iteracijų. Atlikus 3 iteracijas dar reikia padauginti gautą rezultatą iš x, kad gauti x^(1/2). Šiuolaikiniai procesoriai skaičiuojantis pagal Goldschmidt’o algorithm'ą moka suskaičiuoti vienu metu ir 1/x^(1/2) ir x^(1/2), todėl sutaupoma viena daugybos operacija (palyginus su Niutono metodu), bet nesutaupomi tranzistoriai. Trims iteracijoms reikia 4*3=12 daugybos operacijų ir 1*3=3 atimties operacijos. Taigi, kvadratinę šaknį reikia traukti 5 kartus. Tai iš viso reikia 12*5=60 daugybos operacijų ir 3*5=15 atimties operacijų. Gaunasi, kad iš viso reikia 60+15=75 operacijų, kad penkis kartus ištraukti kvadratinę šaknį (ir tai paskui pritaikyti ln(1+x) skaičiavimuose).

Taigi, kad apskaičiuoti 1/z^(1/32) reikia 27 operacijų. Gauname, kad penkių kartų skaičiavimui 1/x^(1/2) reikia 75/27=2.77778 karto daugiau operacijų negu 1/z^(1/32) skaičiavimui. Ne tiek jau ir daug ir todėl greičiausiai nieks neeikvos papildomų tranzistorių $y_{1}=1/z^{1/32}$ lentelėms bei nekurs naujų instrukcijų, kad greičiau apskaičiuoti ln(x). Tiesa, ln(x), kaip parodyta čia ( https://lt.wikibooks.org/wiki/Matematika/Teiloro_eilutė#e_pakelta_x_liekamojo_nario_įvertinimas ), taikomas skaičiaus kelimui bet kokiu laipsniu. Bet skaičiaus kelimas bet kokiu laipsniu mažai kur reikalingas. Dar yra ir kitokių gana konkurencingų ln(x) skaičiavimo būdų.

Geriau paskaičius, Goldschmidt’o algorithm'as gali būti, kad sutaupo kiekvienoje iteracijoje daugybą iš 0.5. Tai iš viso trijose iteracijose sutaupomos trys daugybos operacijos (arba trys dalybos operacijos iš 2; dvejetainėje sistemoje dalyba iš dviejų yra pastumimas visų bitų į dešinę per vieną bito poziciją; bet kas greičiau ar daugyba, ar pastumimas tai yra dar neaišku).

Va čia gerai paaiškinta apie Goldschmidt’o algorithm'ą: https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Goldschmidt’s_algorithm

• Apskaičiuosime ${\frac {1}{\sqrt[{32}]{z}}}={\frac {1}{\sqrt[{32}]{0.75}}}$ ir rasime kam lygus ${\sqrt[{32}]{0.75}}={\sqrt[{32}]{z}}.$ Tiksli ${\sqrt[{32}]{0.75}}$ reikšmė yra tokia:

${\sqrt[{32}]{0.75}}=$ 0.99105022504172762555415281377767.

O $y={\frac {1}{\sqrt[{32}]{0.75}}}=$ 1.0090305967670766121721228310046.

Mes paimsime, kad $y_{1}=1.009.$ Tada pagal (5) formulę

$y_{2}={\frac {y_{1}\left(33-zy_{1}^{32}\right)}{32}}={\frac {1.009\left(33-0.75\cdot 1.009^{32}\right)}{32}}=$

= 1.009*(33 - 0.75*1.3320401627616656376316134288889)/32 = 1.0090305814634409227725711418028.

[1.0090305967670766121721228310046 - 1.0090305814634409227725711418028 = 0.0000000153036356893995516892018 $\approx 1.53\cdot 10^{-8}.$]

Toliau pagal (5) formulę rasime $y_{3}$:

$y_{3}={\frac {y_{2}\left(33-zy_{2}^{32}\right)}{32}}=$

= 1.0090305814634409227725711418028*(33 - 0.75*1.0090305814634409227725711418028^32)/32 = 1.0090305967670727824366475789071.

[1.0090305967670766121721228310046 - 1.0090305967670727824366475789071 = 0.0000000000000038297354752520975 $\approx 3.8297\cdot 10^{-15}.$]

Toliau pagal (5) formulę rasime $y_{4}$:

$y_{4}={\frac {y_{3}\left(33-zy_{3}^{32}\right)}{32}}=$

= 1.0090305967670727824366475789071*(33 - 0.75*1.0090305967670727824366475789071^32)/32 = 1.0090305967670766121721228307648.

Atėmę gautą $y_{4}$ reikšmę iš tikslios $y={\frac {1}{\sqrt[{32}]{0.75}}}$ reikšmės, gauname:

1.0090305967670766121721228310046 - 1.0090305967670766121721228307648 = 0.0000000000000000000000000002398 $\approx 2.398\cdot 10^{-28}.$

Tik po 3 iteracijų paklaida yra mažesnė nei $10^{-27}.$ Taip yra todėl, nes parinkome $y_{1}=1.009$ su paklaida

1.0090305967670766121721228310046 - 1.009 = 0.0000305967670766121721228310046 $\approx 3.0596767\cdot 10^{-5}.$

Na o,

${\sqrt[{32}]{z}}={\sqrt[{32}]{0.75}}\approx 1/y_{4}=$ 1/1.0090305967670766121721228307648 = 0.99105022504172762555415281401323.

[0.99105022504172762555415281377767 - 0.99105022504172762555415281401323 = -0.00000000000000000000000000023556 =

$\approx -2.3556\cdot 10^{-28}.$ Čia kažkodėl iš viso, su pirmu nuliu prieš tašką, yra 33 skaitmenys, o ne 32 skaitmenys kaip prieš tai esančiame $y_{4}$ paklaidos atsakyme.]

Toliau apskaičiuosime ln(0.75). Kadangi

$\ln 0.75=-\ln 0.75^{-1}=-32\ln 0.75^{-1/32}=-32\ln {\frac {1}{z^{1/32}}}\approx -32\ln y_{4},$

tai tereikia apskaičiuoti $\ln y_{4}.$ Gauname

$\ln y_{4}=$ ln(1.0090305967670766121721228307648) = ln(1+0.0090305967670766121721228307648), x=0.0090305967670766121721228307648.

Pagal Teiloro formulę apskaičiuosime -32*ln(1+x) = -32*ln(1+0.0090305967670766121721228307648). Imsime n=9. Tada

$\ln(1+x)=x-{\frac {x^{2}}{2}}+{\frac {x^{3}}{3}}-{\frac {x^{4}}{4}}...+(-1)^{n-1}{\frac {x^{n}}{n}}+r_{n+1}(x);\quad (13)$

čia $|-32r_{n+1}(x)|={\Big |}-32{\frac {x^{n+1}}{(1+c)^{n+1}(n+1)}}{\Big |},\;\;0<c<0.009030596767.$

Vietoje x įstatome 0.0090305967670766121721228307648 į (13) formulę ir su n=9, gauname

$\ln(1+x)=\ln(1+0.0090305967670766121721228307648)=$

= 0.0090305967670766121721228307648 - 0.0090305967670766121721228307648^2/2 + 0.0090305967670766121721228307648^3/3 - 0.0090305967670766121721228307648^4/4 + 0.0090305967670766121721228307648^5/5 - 0.0090305967670766121721228307648^6/6 + 0.0090305967670766121721228307648^7/7 - 0.0090305967670766121721228307648^8/8 + 0.0090305967670766121721228307648^9/9 =

= 0.00899006476411815398283337192369.

[0.00899006476411815398283337192369 - ln(1.0090305967670766121721228307648) =

= 0.00899006476411815398283337192369 - 0.00899006476411815398247559369963 =

= 0.00000000000000000000035777822406.]

Gavome paklaidą $3.5777822406\cdot 10^{-22}<10^{-21}.$

Įvertinsime liekamąjį narį teoriškai.

$|-32r_{9+1}(x)|<{\Big |}-32\cdot {\frac {0.009030596767^{9+1}}{(1+0)^{9+1}(9+1)}}{\Big |}=$

=| -32*(0.0090305967670766121721228307648^10)/10 | =

(čia galima padaryti "Copy-Paste", viską nukopijavus į kalkuliatoriaus laukelį kur rodo skaičius, išskyrus modulio ženklus "|" ir paspaudus paskui kalkuliatoriuje ženklą "="; tokiu budu bus gautas žemiau esantis atsakymas)

=| -1.1542888083111938172050915484026e-20 | $<10^{-19}.$

Padauginsime anksčiau per Teiloro eilutę gautą atsakymą ln(1+x)=0.00899006476411815398283337192369 iš -32, kad gauti ln(0.75). Taigi,

$\ln(0.75)\approx$ -32*0.00899006476411815398283337192369 = -0.28768207245178092745066790155808.

Atėmę apskaičiuotą ln(0.75) [apytikslią] reikšmę iš Windows 10 kalkuliatoriaus tikslios reikšmės, gauname paklaidą:

ln(0.75) - (-0.28768207245178092745066790155808) = -0.28768207245178092743921900599383 + 0.28768207245178092745066790155808 =

= 0.00000000000000000001144889556425.

Taigi, gavome paklaidą $1.144889556425\cdot 10^{-20}<|-32r_{n+1}(x)|=|-1.1542888083\cdot 10^{-20}|<10^{-19}.$

Viskas atitinka teoriją.

Goldschmidt dalybos algoritmas

https://en.wikipedia.org/wiki/Division_algorithm#Goldschmidt_division

Goldschmidt'o algoritmas veikia taip, kad skaitiklis ir vardiklis dauginamas iš bendro skaičiaus F_i. Skaičius F_i pasirinktas taip, kad dauginant vardiklį kiekvieną kartą iš F_i, vardiklis konverguoja į vienetą. Taip kelis kartus dauginant skaitiklį ir vardiklį iš F_i skaitiklis konverguoja į Q=N/D; čia N yra dalinys (skaitiklis), D yra daliklis (vardiklis). Formulė Q skaičiuoti yra tokia:

$Q={\frac {N}{D}}{\frac {F_{1}}{F_{1}}}{\frac {F_{2}}{F_{2}}}{\frac {F_{\ldots }}{F_{\ldots }}}.$

Žingsniai Goldschmidt'o dalybai yra:

1. Generuoti F_i .
2. Padauginti dalinį ir daliklį iš F_i .
3. Jeigu daliklis pakankmai arti prie 1, pateikti dalmenį Q, priešintu atveju, kartoti nuo žingsnio 1.

Pirma reikia daliklį (vardiklį) padauginti iš 0.1, kad jisai būtų intervale nuo 0 iki 1 (0<D<1). Kiekvienas F_i randamas iš D_i (D₁ = D):

$F_{i+1}=2-D_{i}.$

Dauginant dalinį $N_{i}$ ir daliklį $D_{i}$ iš faktoriaus F_i duoda $N_{i+1}$ ir $D_{i+1}$:

${\frac {N_{i+1}}{D_{i+1}}}={\frac {N_{i}}{D_{i}}}{\frac {F_{i+1}}{F_{i+1}}}.$

Po pakankmo iteracijų skaičiaus k, $Q=N_{k}$.

Pavyzdžiai

• Apskaičiuosime Q=N/D=3/7. Padauginame 7 iš 0.1: $D_{1}=0.1\cdot 7=0.7.$ Tada

$F_{i+1}=2-D_{i},$

$F_{1+1}=F_{2}=2-D_{1}=2-0.7=1.3.$

Toliau padauginame $N_{1}=3$ ir $D_{1}=0.7$ iš $F_{2}$:

$N_{2}=N_{1}\cdot F_{2}=3\cdot 1.3=3.9,$

$D_{2}=D_{1}\cdot F_{2}=0.7\cdot 1.3=0.91.$

Toliau apskaičiuojame $F_{3}$:

$F_{3}=2-D_{2}=2-0.91=1.09.$

Toliau padauginame $N_{2}$ ir $D_{2}$ iš $F_{3}$:

$N_{3}=N_{2}\cdot F_{3}=3.9\cdot 1.09=4.251,$

$D_{3}=D_{2}\cdot F_{3}=0.91\cdot 1.09=0.9919.$

Toliau apskaičiuojame $F_{4}$:

$F_{4}=2-D_{3}=$ 2 - 0.9919 = 1.0081.

Toliau padauginame $N_{3}$ ir $D_{3}$ iš $F_{4}$:

$N_{4}=N_{3}\cdot F_{4}=$ 4.251*1.0081 = 4.2854331,

$D_{4}=D_{3}\cdot F_{4}=$ 0.9919*1.0081 = 0.99993439.

Toliau apskaičiuojame $F_{5}$:

$F_{5}=2-D_{4}=$ 2 - 0.99993439 = 1.00006561.

Toliau padauginame $N_{3}$ ir $D_{3}$ iš $F_{4}$:

$N_{5}=N_{4}\cdot F_{5}=$ 4.2854331*1.00006561 = 4.285714267265691,

$D_{5}=D_{4}\cdot F_{5}=$ 0.99993439*1.00006561 = 0.9999999956953279.

Taigi $Q=N_{5}\cdot 0.1=$ 4.285714267265691*0.1 = 0.4285714267265691. Tiksli Q=3/7 reikšmė yra tokia:

Q = 3/7 = 0.42857142857142857142857142857143.

[0.42857142857142857142857142857143 - 0.4285714267265691 = 0.00000000184485947142857142857143 $\approx 1.84485947\cdot 10^{-9}.$]

Po 4 iteracijų gavome paklaidą mažesnę nei $10^{-8}.$

Dar padarysime penktą įteraciją, apskaičiuodami $F_{6},\;$ $N_{6}$ ir $D_{6}$:

$F_{6}=2-D_{5}=$ 2 - 0.9999999956953279 = 1.0000000043046721;

$N_{6}=N_{5}\cdot F_{6}=$ 4.285714267265691*1.0000000043046721 = 4.2857142857142856348705633349211,

$D_{6}=D_{5}\cdot F_{6}=$ 0.9999999956953279*1.0000000043046721 = 0.99999999999999998146979811148159.

Taigi $Q=N_{6}\cdot 0.1=$ 4.2857142857142856348705633349211*0.1 = 0.42857142857142856348705633349211.

[0.42857142857142857142857142857143 - 0.42857142857142856348705633349211 = 0.00000000000000000794151509507932 =

$\approx 7.941515\cdot 10^{-18}.$]

Po 5 iteracijų gavome paklaidą mažesnę nei $10^{-17}.$

Nuorodos

• https://en.wikipedia.org/wiki/Methods_of_computing_square_roots

• http://www.techmat.vgtu.lt/konspektai/Dskaiciavimas/t_03.pdf Sekos riba.
• http://www.techmat.vgtu.lt/konspektai/Dskaiciavimas/t_04.pdf Konverguojančios sekos ir jų pagrindinės savybės.