Rješavanje jednadžbi

Razmotrit ćemo nekoliko iterativnih postupaka za rješavanje jednadžbe

$$f(x) = 0,$$ (3.1)

uz opću pretpostavku da je $f$ neprekidna funkcija. Pretpostavimo da su $\alpha$ i $\beta{}$ u domeni funkcije $f$ takvi da je

$$f(\alpha{})\,f(\beta{}) < 0.$$ (3.2)

Zbog neprekidnosti funkcije, prema Bolzanovom teoremu (neprekidnost funkcije na segmentu, v. [8, Teorem 4, str. 31]), postoji barem jedno rješenje $s$ jednadžbe (3.1) u segmentu $[\alpha{},\beta{}].$ Iterativni postupak je postupak, kojim nalazimo niz brojeva $x_n,\;n=0,1,2,\ldots{}$ tako, da je

$$\lim_{n\rightarrow{}\infty{}} x_n = s.$$

Član niza $x_n$ se zove $n$-ta aproksimacija rješenja $s.$ Naravno, možemo naći samo konačno mnogo članova niza. Tako se moramo zadovoljiti s približnim rješenjem. Koja će aproksimacija biti dovoljno dobra ovisi o tome kolika je greška dozvoljiva. Prema tome bit će nam važno znati ocijeniti grešku koju činimo kad pravo rješenje $s$ zamjenimo s $n$-tom aproksimacijom $x_n.$

Metoda polovljenja

Metoda polovljenja se sastoji u tome da se segment $[a,b],$ na kojem je ispunjen uvjet $f(a)f(b)<0,$ raspolovi, tj. nađe polovište $c.$ Ako je $f(c)=0,$ onda je $s=c.$ U protivnom se ponovi operacija na onom od segmenata $[a,c]$ ili $[c,b]$ na kojem je ispunjen uvjet (3.2), itd. Tako imamo sljedeći algoritam.

Algoritam 3   , Za $k=0,1,2,\ldots{}$ računamo

$$x_{k+2} = \frac{1}{2}\,(x_{k+1} + x_{k}),$$    ako je $$f(x_k)\,f(x_{k+1}) \leqslant{} 0,$$

$$x_{k+2} = \frac{1}{2}\,(x_{k+1} + x_{k-1}),$$    ako je $$f(x_k-1)\,f(x_{k+1}) \leqslant{} 0,$$

s tim da je

$$x_0 = a,\hspace{1cm}x_1 = b.$$

Ovaj algoritam omogućava sljedeći program u programskom paketu 'Mathematica'.

Mathematica program 1   (Metoda polovljenja)

f[t_]=; (* funkcija *)
a=;
b=;     (* pocetni interval *)
x=(a+b)/2.;
n=0;
Print[{n,a,x,b}];
        While[
        N[f[x]]!=0.,
        Print["   ",N[f[a]],"  ",N[f[x]],"   ",N[f[b]]];
                If[
                f[b]f[x]<0,
                a=x;x=(x+b)/2,
                b=x;x=(a+x)/2
                ];
        n=n+1;
        If[n>100,Break[]  (* prekid petlje nakon 100 koraka *)
        ];
Print[{n,a,x,b}]]

Metoda uvijek konvergira, ali vrlo sporo. Očito je

$$\vert x_2 - s\vert \leqslant{} \frac{1}{2}\,(b-a),\quad \vert x_3 - s\vert \leqslant{} \frac{1}{2^2}\,(b-a),\quad \ldots{}$$

Tako za $k$-tu aproksimaciju imamo ocjenu greške

$$\left\vert x_k - s\right\vert \leqslant{} \frac{1}{2^{k+1}}\,(b-a).$$ (3.3)

Ova ocjena se zove apriorna, jer se može izračunati prije nego smo postupak iteracije uopće započeli. Naravno, ova ocjena nam ne kaže kolika je greška, već samo to da greška nije veća od tog broja.

Primjer 3.1   Riješiti metodom polovljenja jednadžbu

$$x^3 - x + 1 = 0.$$

Rješenje. Računanjem vrijednosti polinoma $f(x)=x^3 - x + 1$ na skupu cijelih brojeva, dobivamo da je $f(-2)=-5,$ a $f(-1)=1.$ Tako na segmentu $[-2,-1]$ postoji bar jedno rješenje. Ako na ovaj problem primijenimo gornji program, dobivamo sljedeću tablicu.

$$k x_k x_{k+2} x_{k+1} f(x_k) f(x_{k+2}) f(x_{k+1})$$

$$0 -2 -1.5 -1 -5. -0.875 1.$$

$$1 -1.5 -1.25 -1 -0.875 0.296875 1.$$

$$2 -1.5 -1.375 -1.25 -0.875 -0.224609 0.296875$$

$$3 -1.375 -1.3125 -1.25 -0.224609 0.0515137 0.296875$$

$$4 -1.375 -1.34375 -1.3125 -0.224609 -0.0826111 0.0515137$$

$$5 -1.34375 -1.32813 -1.3125 -0.0826111 -0.014576 0.0515137$$

$$6 -1.32813 -1.32031 -1.3125 -0.014576 0.0187106 0.0515137$$

$$7 -1.32813 -1.32422 -1.32031 -0.014576 0.00212795 0.0187106$$

$$8 -1.32813 -1.32617 -1.32422 -0.014576 -0.00620883 0.00212795$$

$$9 -1.32617 -1.3252 -1.32422 -0.00620883 -0.00203665 0.00212795$$

$$\vdots \vdots \vdots \vdots \vdots \vdots \vdots$$

$$50 -1.32472 -1.32472 -1.32472 -8.88178\times 10^{-16} -8.88178\times 10^{-16} 2.66454\times 10^{-15}$$

Računalo je tek u pedesetoj aproksimaciji došlo do granica svoje točnosti. Točnost ispisanih rezultata (pet znamanaka) se postiže u sedamnaestoj aproksimaciji, što se moglo unaprijed izračunati. Budući da je $b-a=1,$ greška $k$-te aproksimacije nije veća od $\frac{1}{2^{k+1}},$ pa zbog

$$\frac{1}{2^{17}} \approx 7.62939\times 10^{-6}, \qquad \frac{1}{2^{18}} \approx 3.8147\times 10^{-6},$$

rješenje na pet decimala točno se postiže u sedamnaestoj aproksimaciji.

Zadatak iz primjera 3.1 se može doduše egzaktno riješiti pomoću Cardanovih formula. Ipak te formule nisu tako jednostavne kao formula za rješenje kvadratne jednadžbe, pa se algebarske jednadžbe trećeg stupnja često rješavaju približnim metodama. U sljedećem primjeru se rješava jednadžba koju ne možemo egzaktno riješiti.

Primjer 3.2   Naći najmanje pozitivno rješenje jednadžbe (v. 2.19)

$${\rm tg}\,\lambda = \frac{1}{\lambda}$$

metodom polovljenja.

Rješenje. Iz slike 2.7 se vidi da se traženo rješenje nalazi u intervalu $\langle{}0,\frac{\pi}{2}\rangle{}.$ Neka je

$$f(\lambda{}) = {\rm tg}\,\lambda - \frac{1}{\lambda}.$$

Imamo

$$f\left(\frac{\pi}{3}\right) = {\sqrt{3}} - {\frac{3}{\pi }} > 0, \quad f\left(\frac{\pi}{4}\right) = {\frac{-4 + \pi }{\pi }} < 0.$$

Prema tome, za početni segment možemo uzeti segment $[\frac{\pi}{4},\frac{\pi}{3}].$

Izračunajmo sada koju aproksimaciju treba izračunati da bi se dobila točnost na četiri decimale. Iz formule (3.3) slijedi da treba naći $k$ tako da bude

$$\frac{1}{2^{k+1}}\,(b-a) = \frac{1}{2^{k+1}}\,\frac{\pi}{12} \leqslant \frac{1}{2^{k+1}}\,\frac{1}{3} \leqslant 0.5\times 10^{-6}.$$

Odatle slijedi $k\geqslant 15.0247,$ pa treba izračunati četrnaestu aproksimaciju. Koristeći program 1 dobivamo

$$k x_k x_{k+2} x_{k+1} f(x_k) f(x_{k+2}) f(x_{k+1})$$

$$0 0.785398 0.916298 1.0472 -0.27324 0.211877 0.777121$$

$$1 0.785398 0.850848 0.916298 -0.27324 -0.0350166 0.211877$$

$$2 0.850848 0.883573 0.916298 -0.0350166 0.086735 0.211877$$

$$3 0.850848 0.86721 0.883573 -0.0350166 0.0255155 0.086735$$

$$4 0.850848 0.859029 0.86721 -0.0350166 -0.00482682 0.0255155$$

$$5 0.859029 0.86312 0.86721 -0.0048268 0.0103241 0.0255155$$

$$6 0.859029 0.861075 0.86312 -0.0048268 0.0027437 0.0103241$$

$$7 0.859029 0.860052 0.861075 -0.0048268 -0.0010428 0.0027437$$

$$8 0.860052 0.860563 0.861075 -0.0010428 0.0008502 0.0027437$$

$$9 0.860052 0.860308 0.860563 -0.0010428 -0.000096 0.000850$$

$$10 0.860308 0.860435 0.860563 -0.0000963 0.0003769 0.000850$$

$$11 0.860308 0.860371 0.860435 -0.0000964 0.000140 0.0003769$$

$$12 0.860308 0.86034 0.860371 -0.0000964 0.0000219 0.000140$$

$$13 0.860308 0.860324 0.86034 -0.0000964 -0.000037 0.0000219$$

$$14 0.860324 0.860332 0.86034 -0.000037 -7.633771\times 10^{-6} 0.000021$$

Dakle rješenje, na četiri decimale točno, je

$$\lambda{} = 0.860332.$$

Metoda iteracije

Napišimo jednadžbu (3.1) u obliku

$$x = \varphi(x)$$ (3.4)

Riješiti jednadžbu sada znači naći takav $s\in [a,b]$ da vrijedi $\varphi(s)=s.$ Geometrijski to znači da tražimo presjek grafa funkcije i pravca $y=x.$

Imamo sljedeći algoritam

Algoritam 4   (Metoda iteracije) Izaberemo $x_0\in [a,b],$ i računamo niz $x_n,$ za $n=1,2,3,\ldots{}$ po formuli

$$x_n = \varphi(x_{n-1}).$$ (3.5)

Ako je $x_k = \varphi(x_{k})$ za neki $k,$ onda je $s=x_k.$ U protivnom nastavljamo računanje daljnih članova niza.

Ovaj algoritam dobro opisuje sljedeća slika.

Imamo sljedeći program za metodu iteracije

Mathematica program 2   (Metoda iteracije)

varphi[t_]=; (* funkcija *)
x=;          (* pocetna aproksimacija *)
n=0;
   While[
   N[f[x]]!=x,
   x=varphi[x];
   n=n+1;
   If[n>100,Break[]]; (* prekid procesa nakon 100 iteracija *)
   Print[N[x]]
   ]

Može se dogoditi da ne postoji rješenje,

da postoji ali da ga ne možemo dostići ovom metodom,

ili da postoji i da nam je dostupno.

Evo još nekoliko slika, koje pokazuju različito ponašanje iteracijskog niza.

Da bi ova metoda funkcionirala, moraju brojevi $x_1,x_2,x_3,\ldots{}$ biti u $[a,b]$ za bilo koju početnu aproksimaciju. To će svakako biti ispunjeno, ako je $\varphi([a,b])\subset [a,b].$ A kako je funkcija $\varphi$ još i neprekidna, onda postoji bar jedna točka presjeka.

Teorem 23   Neka je

1.

$\varphi:[a,b]\rightarrow{}R$ funkcija klase $C^1$ na $[a,b],$

2.

$\varphi(x)\in{}[a,b]$ za svaki $x\in{}[a,b],$

3.

$\left\vert\varphi'(x)\right\vert\leqslant{}L<1,\hspace{1cm}\forall x\in{}[a,b].$

Tada za proizvoljan $x_0\in{}[a,b]$ niz $(x_n)$ definiran algoritmom (3.5) konvergira k jedinstvenom rješenju jednadžbe $x=\varphi(x).$

Dokaz. Dokažimo najprije da postoji barem jedno rješenje. Stavimo $g(x)=x-\varphi(x).$ Prema uvjetima teorema, $g(a)\leqslant{}0,$ i $g(b)\geqslant{}0.$ Budući da je $g$ neprekidna funkcija na segmentu $[a,b],$ slijedi da postoji bar jedan $s\in [a,b]$ takav da je $g(s)=0,$ tj. $s-\varphi(s)=0,$ dakle $s$ je rješenje.
Dokažimo sada da ne postoji više od jednog rješenja. Pretpostavimo da su $s_1$ i $s_2$ dva međusobno različita rješenja. Tada, prema Lagrangeovom teoremu srednje vrijednosti,

$$\vert s_1-s_2\vert = \vert\varphi(s_1)-\varphi(s_2)\vert \leqslan... ...vert\,\vert s_1-s_2\vert \leqslant{}L\,\vert s_1-s_2\vert < \vert s_1-s_2\vert.$$

što je kontradikcija. Dakle problem ima jedno i samo jedno rješenje $s.$
Preostaje dokazati da za bilo koju početnu aproksimaciju $x_0$ niz $(x_n)$ konvergira k rješenju $s.$ Imamo

$$\vert x_n - s\vert = \vert\varphi(x_{n-1}) - s\vert = \vert\varph... ...= \vert\varphi'(\xi)\,(x_{n-1} - s)\vert \leqslant{} L\,\vert x_{n-1} - s\vert.$$

Odatle

$$\vert x_n - s\vert \leqslant{} L\,\vert x_{n-1} - s\vert \leqslan... ...,\vert x_{n-2} - s\vert \leqslant{} \ldots \leqslant{} L^n\,\vert x_0 - s\vert.$$

Budući da je $0\leqslant{}L<1,$ slijedi

$$\lim_{n\rightarrow{}\infty} L^n = 0,$$

pa je prema tome

$$\lim_{n\rightarrow{}\infty} \vert x_n - s\vert = 0,$$

tj.

$$\lim_{n\rightarrow{}\infty{}} x_n = s.$$

$\heartsuit$

Apriornu ocjenu greške dobivamo na sljedeći način. Iz

$$\vert x_{n+1}-x_n\vert=\vert\varphi(x_n)-\varphi(x_{n-1})\vert = \vert\varphi'(\sigma{})\,(x_{n}-x_{n-1})\vert \leqslant{}L\,\vert x_n-x_{n-1}\vert$$

slijedi

$$\vert x_{n+1}-x_n\vert \leqslant{}L\,\vert x_n-x_{n-1}\vert \leqslant{}L^2\,\vert x_{n-1}-x_{n-2}\vert$$

$$\leqslant{}L^3\,\vert x_{n-2}-x_{n-3}\vert \leqslant{}\ldots\leqslant{}L^n\,\vert x_1-x_0\vert.$$

Za proizvoljni prirodni broj $m>n,$ imamo

$$\vert x_m-x_n\vert\leqslant{}\vert x_m-x_{m-1}\vert+\vert x_{m-1}-x_{m-2}\vert+\cdots +\vert x_{n+1}-x_n\vert.$$

Iz gornje nejednakosti slijedi

$$\vert x_m-x_n\vert\leqslant{}\left(L^{m-1}+L^{m-2}+\cdots +L^n\right)\vert x_{1}-x_0\vert.$$

Zbog $0<L<1,$

$$L^{m-1}+L^{m-2}+\cdots +L^n=L^n\left(1+L+L^2+\cdots +L^{m-n-1}\right)$$

$$\leqslant{}L^n\left(1+L+L^2+\cdots\right)=L^n\,\frac{1}{1-L},$$

po formuli za sumu geometrijskog reda. Tako je

$$\vert x_m-x_n\vert\leqslant{}L^n\,\frac{1}{1-L}\,\vert x_1-x_0\vert.$$

Desna strana ove nejednakosti ne ovisi o $m,$ pa prema tome

$$\lim_{m\rightarrow{}\infty{}}\vert x_m-x_n\vert=\vert s-x_n\vert \leqslant{} \frac{L^n}{1-L}\,\vert x_1-x_0\vert.$$

Dakle apriorna ocjena greške $n$-te aproksimacije je

$$\vert x_n-s\vert\leqslant \frac{L^n}{1-L}\,\vert x_1-x_0\vert.$$ (3.6)

Aposteriorna ocjena greške je ocjena koja se računa pomoću dobivenih aproksimacija. U ovom slučaju ona je

$$\left\vert x_k - s\right\vert \leqslant{} \frac{L}{1-L}\left\vert x_k - x_{k-1}\right\vert.$$ (3.7)

Doista,

$$\left\vert x_k - x_{k+p}\right\vert = \left\vert x_k - x_{k+1} + x_{k+1} - x_{k+2} + \cdots{} + x_{k+p-1} - x_{k+p}\right\vert$$

$$\leqslant{} \left\vert x_k - x_{k+1}\right\vert + \left\vert x_{k+1} - x_{k+2}\right\vert + \cdots{} + \left\vert x_{k+p-1} - x_{k+p}\right\vert$$

$$\leqslant{} L\,\left\vert x_{k-1} - x_k\right\vert + L^2\,\left\vert x_{k-1} - x_k\right\vert + \cdots{} + L^p\,\left\vert x_{k-1} - x_k\right\vert$$

$$= \left(L + L^2 + \cdots{} + L^p\right)\,\left\vert x_{k-1} - x_k\right\vert = L\,\frac{1 - L^p}{1 - L}\,\left\vert x_{k-1} - x_k\right\vert.$$

$$\lim_{p\rightarrow{}\infty{}} \left\vert x_k - x_{k+p}\right\vert = \left\vert x_k - s\right\vert.$$

$$\lim_{p\rightarrow{}\infty{}} L\,\frac{1 - L^p}{1 - L}\,\left\vert x_{k-1} - x_k\right\vert = \frac{L}{1-L}\left\vert x_k - x_{k-1}\right\vert.$$

Tako je

$$\left\vert x_k - s\right\vert \leqslant{} \frac{L}{1-L}\left\vert x_k - x_{k-1}\right\vert.$$

Primjer 3.3   Riješiti metodom iteracije jednadžbu u primjeru 3.1.

Rješenje. Iz

$$x^3 - x + 1 = 0$$

slijedi

$$x = \sqrt[3]{x-1}.$$

Tako je

$$\varphi(x) = \sqrt[3]{x-1},\qquad \varphi'(x) = {\frac{1}{3\,{{\left( x-1 \right) }^{{\frac{2}{3}}}}}}$$

Rješenje postoji na segmentu $[-2,-1]$ (v. primjer 3.1). $\varphi'$ je pozitivna, pa funkcija $\varphi$ raste.

$$\varphi(-2) = \sqrt[3]{-3} = -\sqrt[3]{3}\in [-2,-1],\quad \varphi(-1) = \sqrt[3]{-2} = -\sqrt[3]{2}\in [-2,-1].$$

Osim toga $\varphi'$ raste na $[-2,-1],$ pozitivna je, pa najveću vrijednost ima u $-1,$ i to

$$\varphi'(-1) = \frac{1}{3\,\sqrt[3]{4}} < \frac{1}{3}.$$

Ova diskusija pokazuje da su uvjeti teorema 23 ispunjeni, pa će iteracijski proces konvergirati bez obzira na to koji broj iz $[-2,-1]$ uzmemo kao početnu aproksimaciju. Ujedno nam ocjena

$$L=\max_{x\in [-2,-1]}\vert\varphi'(x)\vert\leqslant{}\varphi'(-1) \leqslant{}\frac{1}{3}$$

može poslužiti da bismo apriorno ocijenili grešku $k$-te aproksimacije.

Imamo na pr.

$$x_0 = -1, \qquad x_1 = \varphi(-1) = -\sqrt[3]{2} = -1.25992,$$

pa greška $k$-te aproksimacije nije veća od

$$\frac{\frac{1}{3^k}}{1-\frac{1}{3}}\,\vert x_1-x_0\vert = {\frac{{3^{1 - k}}}{2}}\,0.25992 = 0.38988\times 3^{-k}.$$

Ako želimo rješenje točno na pet decimala, izlazi da mora biti $k>12.349,$ dakle trinaesta aproksimacija daje traženu točnost. Pomoću programa 2 nalazimo da je

$$x_{13} = -1.32472.$$

Primjer 3.4   Riješiti zadatak u primjeru 3.2 metodom iteracije.

Rješenje. Jednadžbu treba napisati u obliku

$$\lambda{} = \varphi(\lambda{}).$$

Ako stavimo

$$\lambda{} = \frac{1}{{\rm tg}\,\lambda{},}$$

onda je

$$\varphi(\lambda) = \frac{1}{{\rm tg}\,\lambda},$$

$$\varphi'(\lambda) = -\frac{1}{\sin^2 \lambda},$$

pa je

$$\vert\varphi'(\lambda)\vert \geq 1,$$   za svaki $$\lambda.$$

To ne odgovara uvjetima teorema 23. Jednadžbu možemo i drukčije napisati

$$\lambda{} = {\rm Arctg}\,\frac{1}{\lambda{}}.$$

Tada je

$$\varphi(\lambda) = {\rm Arctg}\,\frac{1}{\lambda},$$

$$\varphi'(\lambda) = -\frac{1}{1 + \lambda^2},$$

pa je

$$\vert\varphi'(\lambda)\vert \leq 1,$$   za svaki $$\lambda.$$

Specijalno, $\vert\varphi'(\lambda)\vert=1$ samo za $\lambda=0.$ Domena od $\varphi$ je $R\setminus\{0\}.$ Budući da je $\varphi'(\lambda)<0,$ funkcija pada na svakom od intervala $\langle -\infty,0\rangle,$ $\langle{}0,\infty{}\rangle{}.$ Nas interesiraju samo pozitivna rješenja, pa nam je interesantan samo interval $\langle{}0,\infty{}\rangle{}.$ Na tom intervalu

$$\lim_{\lambda{}\rightarrow{}0+} {\rm... ...ad \lim_{\lambda{}\rightarrow{}\infty{}} {\rm Arctg}\, \frac{1}{\lambda{}} = 0.$$

Dakle $\varphi$ preslikava $\langle{}0,\infty{}\rangle{}$ na $\langle{}0,\frac{\pi}{2}\rangle{}.$ Također

$${\rm Arctg}\,1 = \frac{\pi}{4},\qquad {\rm Arctg}\,\frac{\pi}{4} = 0.905023 < 1,$$

pa $\varphi$ preslikava segment $[\frac{\pi}{4},1]$ u samog sebe. Apsolutna vrijednost derivacije je najveća na lijevom rubu. Tako možemo staviti

$$L = \frac{1}{1+\left(\frac{\pi}{4}\right)^2} = \frac{4}{4+\pi^2} \leqslant{} \frac{4}{13} \leqslant{} \frac{1}{3}.$$

To znači da će metoda iteracije konvergirati uzmemo li bilo koji broj iz segmenta $[\frac{\pi}{4},1]$ kao početnu aproksimaciju. Pomoću programa 2 nalazimo da je zaokruženo na šest decimala, uz $x_0=1,$

$$x_{21} = 0.860332,\quad x_{22} = 0.860334,\quad x_{23} = 0.860333,$$

$$x_{24} = 0.860334,\quad x_{25} = 0.860333,\quad x_{26} = 0.860334,$$

i dalje se ovaj broj ponavlja.

Newtonova metoda

Neka je $f:[a,b]\rightarrow{}R$ klase $C^2$ na $[a,b].$ Želimo riješiti jednadžbu

$$f(x) = 0.$$

Newtonova metoda, koja se često zove Newton-Raphsonova metoda ili metoda tangente, sastoji se u tome da se $n+1$-va aproksimacija $x_{n+1}$ odredi kao sjecište tangente na graf funkcije $f$ u točki s apscisom $x_n$ s osi $x.$

Jednadžba tangente na graf funkcije $f$ u točki $(x_n,f(x_n))$ je

$$y = f(x_n) + f'(x_n)\,(x-x_n).$$

Kad stavimo $y=0,$ dobijemo njezino sjecište s osi

$$x=x_n-\frac{f(x_n)}{f'(x_n)}.$$

Dakle $n+1$-va aproksimacija se računa po formuli

$$x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}.$$

Time smo dobili sljedeći algoritam.

Algoritam 5   (Newtonova metoda) Izaberemo $x_0\in [a,b],$ i računamo niz $x_n,$ za $n=1,2,3,\ldots{}$ po formuli

$$x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)},\hspace{1cm}n=1,2,3,\ldots\ .$$ (3.8)

Ako je $f(x_{k}) = 0$ za neki $k,$ onda je $s=x_k.$ U protivnom nastavljamo računanje daljnih članova niza.

Mathematica program 3   (Newtonova metoda)

f[t_]=; (* funkcija *)
x=;     (* pocetna aproksimacija *)
n=0;
   While[
   N[f[x]]!=0.,
   x=x-f[x]/f'[x];
   n=n+1;
   If[n>100,Break[]]; (* prekid nakon 100 iteracija *)
   Print[N[x]]
   ]

Primjer 3.5   Newtonovom metodom riješiti jednadžbu iz primjera 3.1.

Rješenje. Pomoću programa 3 s početnom aproksimacijom $x_0=-1$ dobivamo

$$x_1 = -1.5,\quad x_2 = -1.34783,\quad x_3 = -1.3252,\quad x_4 = -1.32472,\quad x_5 = -1.32472,$$

$$x_6 = -1.32472,\quad x_7 = -1.32472,\quad x_8 = -1.32472,\quad x_9 = -1.32472,\quad x_{10} = -1.32472$$

i to je sve, jer se nakon 10 iteracija postiže strojna točnost. To pokazuje da Newtonova metoda iteracije u ovom slučaju vrlo brzo konvergira.

Primjer 3.6   Riješiti problem iz primjera 3.2 Newtonovom metodom.

Rješenje. Pomoću programa 3 s početnom aproksimacijom $x_0=1$ dobivamo

$$x_1 = 0.874047,\quad x_2 = 0.8604,\quad x_3 = 0.860334,\quad x_4 = 0.860334.$$

Program je napravljen tako da ponavlja iteracije 100 puta ili završava postupak kad postigne strojnu točnost. Ovo pokazuje da je u ovom primjeru strojna točnost postignuta već u četvrtoj aproksimaciji.

S druge strane ako s istim programom (Newtonovom metodom) pokušamo riješiti jednadžbu

$$(x-1)^{21} = 0,$$

počevši s $x_0=0,$ dobivamo $x_{100}=0.992396.$ To pokazuje da konvergencija nije uvijek tako brza kao u prethodnim primjerima. Dapače, ako pokušamo na isti način riješiti jednadžbu

$${\rm Arctg}\,x = 0$$

počevši s $x_0=2,$ dobivamo redom $x_1=-3.53574,$ $x_2=13.951,$ $x_3=-279.344,$ $x_4=122017,$ $x_5=-2.3386\times 10^{10},$ itd. Aproksimacije se sve više udaljavaju od rješenja koje je očito $x=0.$ Sljedeći teorem je jedan od mnogih koji daju dovoljne uvjete da Newtonova metoda konvergira. Uvjerite se da su uvjeti teorema ispunjeni za primjer 3.5 na segmentu $[-2,-1].$

Teorem 24   Neka je $f:[a,b]\rightarrow{}R$ klase $C^2[a,b],$ i neka vrijedi

1.

$f(a)\,f(b)<0,$

2.

$f'(x)\neq 0$ za svaki $x\in{}[a,b],$

3.

$f''(x)\leqslant{}0$ ili $f''(x)\geqslant{}0$ za svaki $x\in{}[a,b],$

4.

ako je $c$ ona rubna točka segmenta $[a,b],$ u kojoj je $\vert f'(x)\vert$ manji, onda je

$$\left\vert\frac{f(c)}{f'(c)}\right\vert\leqslant{}b-a.$$

Tada Newtonova metoda konvergira k jedinstvenom rješenju jednadžbe $f(x)=0$ za bilo koju početnu aproksimaciju $x_0\in{}[a,b].$

Dokaz. Kombinirajući uvjete iz teorema, mogli bismo promatrati razne slučajeve, no može se pokazati da se svaki od njih, uzimajući $-f$ umjesto $f$ i/ili $-x$ umjesto $x,$ može svesti na slučaj

$$f(a)<0,\hspace{1cm}f(b)>0,\hspace{1cm}f''(x)\leqslant{}0,\hspace{1cm}c=b.$$

Iz 2. uvjeta i neprekidnosti funkcije $f'$ slijedi $f'(x)>0$ ili $f'(x)<0$ za svaki $x\in{}[a,b].$ Prema tome funkcija $f$ na segmentu $[a,b]$ raste ili pada. Kako je $f(b)>f(a),$ funkcija $f$ raste. Osim toga iz $f''(x)\leqslant{}0$ slijedi da je $f$ konkavna. Dakle graf funkcije $f$ izgleda kao na slici

Osim toga za proizvoljan $\xi\in[a,b]$ vrijedi

$$f(x) = f(\xi)+f'(\xi)\,(x-\xi)+\frac{f''(x_{\xi})}{2!}\,(x-\xi),\hspace{1cm} \forall{}x\in[a,b],$$

pa je

$$f(x)\leqslant{}f(\xi)+f'(\xi)\,(x-\xi),\hspace{1cm}\forall{}x\in[a,b],$$ (3.9)

tj. tangenta u bilo kojoj točki se nalazi iznad grafa funkcije $f.$
Razmotrimo najprije slučaj $a\leqslant{}x_0\leqslant{}s,$ gdje je $s$ jedinstveno rješenje jednadžbe $f(x)=0.$ Za takav $x_0$ je $f(x_0)\leqslant{}0,$ i $f'(x_0)>0.$ Prema tome je

$$x_1 = x_0-\frac{f(x_0)}{f'(x_0)} \geqslant{} x_0.$$

Iz (3.9) slijedi da je

$$f(x_0)+f'(x_0)\,(s-x_0)\geqslant{}f(s)=0,$$

pa je $x_0\leqslant{}x_1\leqslant{}s.$ Polazeći od $x_1,$ na isti način možemo zaključiti da je $x_1\leqslant{}x_2\leqslant{}s,$ odnosno $x_0\leqslant{}x_1\leqslant{}x_2\leqslant{}s.$ Nastavljajući tako vidimo da je $x_0\leqslant{}x_1\leqslant{}x_2\leqslant{}x_3\leqslant{} \ldots \leqslant{}x_n\leqslant{}\ldots \leqslant{}s.$ Tako smo dobili rastući niz aproksimacija odozgo ograničen. Takav niz ima limes

$$\lim_{n\rightarrow{}\infty{}} x_n = x^{\star} \leqslant{}s.$$

Zbog neprekidnosti funkcije $f$ i njezine derivacije $f',$ imamo

$$x^{\star} = \lim_{n\rightarrow{}\infty{}} x_{n+1}=\lim_{n\rightar... ..._{n\rightarrow{}\infty{}} x_n)} = x^{\star}-\frac{f(x^{\star})}{f'(x^{\star})}.$$

Odatle slijedi $f(x^{\star})=0,$ pa je tako $x^{\star}=s.$
Neka je sada $x_0>s.$ Za $x_0=b$ imamo

$$x_1(b) = b - \frac{f(b)}{f'(b)} = b - \left\vert\frac{f(b)}{f'(b)}\right\vert \geqslant{} b - (b - a) \geqslant{} a.$$

Ako je $s<x_0\leqslant{}b,$ onda iz

$$x_1 = x_0 - \frac{f(x_0)}{f'(x_0)},$$

i (3.9) slijedi

$$0 = f(x_0)+f'(x_0)\,(x_1-x_0)\leqslant{}f(b) + f'(b)\,(x_0-b) + f'(x_0)\,(x_1-x_0)$$

$$\leqslant{}f(b) + f'(b)\,(x_0-b) + f'(b)\,(x_1-x_0) = f(b) + f'(b)\,(x_1-b).$$

Prema tome tangenta u točki $b$ prolazi iznad točke $x_1,$ pa se njezino sjecište $x_1(b)$ s osi $x$ nalazi lijevo od $x_1.$ Tako je

$$a\leqslant{}x_1(b)\leqslant{}x_1\leqslant{}s.$$

Sada smo u situaciji iz prvog dijela dokaza, pa dalje dokaz ide kao tamo. $\heartsuit$

Primjer 3.7   Neka je $k>1$ prirodan broj, i neka je $c>0.$ Nađimo, pomoću Newtonove metode, približnu vrijednost pozitivnog $k$-tog korijena iz $c.$

Rješenje. Izračunati $k$-ti korijen iz broja $c$ znači riješiti po $x$ jednadžbu

$$x^k - c = 0.$$

Ovdje je $f(x)=x^k - c,$ pa Newtonova metoda daje

$$x_{n+1} = x_n - \frac{x_n^k-c}{k\,x_{n}^{k-1}},$$

odnosno

$$x_{n+1} = \frac{1}{k}\,\left((k-1)\,x_n + \frac{c}{x_n^{k-1}}\right).$$

Što se tiče izbora početne aproksimacije $x_0$ i konvergencije, primijetimo sljedeće. Za $0<a<\sqrt[k]{c}<b$ imamo $f(a)f(b)<0.$ Zatim, zbog $f'(x)=kx^{k-1},$ $f''(x)=k(k-1)x^{k-2},$ za svaki $x\in [a,b]$ je $f'(x)>0,$ i $f''(x)>0.$ Na kraju, iz pozitivnosti funkcije $f''$ slijedi rast funkcije $f',$ pa je $a$ onaj rub segmenta u kojem $f'$ ima manju vrijednost. Da bi vrijedilo

$$\left\vert\frac{f(a)}{f'(a)}\right\vert \leqslant{}b-a$$

mora biti

$$b\geqslant{}\frac{a\,\left(1-\frac{c}{a^k} + k \right)}{k}.$$

Dakle za tako veliki $b$ ispunjeni su svi uvjeti teorema 24. Kako $b$ smijemo uzeti još veći, i kako $a$ može biti bilo koji broj veći od $0,$ i manji od $\sqrt[k]{c},$ slijedi da postupak konvergira za svaki $x_0>0.$

Specijalno kad je $k=2,$ imamo formulu za približno računanje drugog korijena

$$x_{n+1} = \frac{1}{2}\,\left(x_n + \frac{c}{x_n}\right).$$ (3.10)

U programskom paketu Mathematica se ovaj postupak programira vrlo jednostavno

Mathematica program 4   ($k$-ti korijen)

Map[Nest[((k-1) # + c/#^(k-1))/k&,x0,#1]&,Range[p]]//N

gdje je p broj aproksimacija koje želimo, a broj x0 je početna aproksimacija. Naravno k je broj korijena koji se vadi.

Uvjerite se na primjerima kako je formula za računanje drugog korijena efikasna.

Primjer 3.8   Za dani pozitivan broj $c$ naći približnu vrijednost njemu recipročnog broja bez dijeljenja.

Rješenje. To je isto kao približno riješiti jednadžbu

$$f(x) = \frac{1}{x} - c = 0.$$

Newtonova metoda daje

$$x_{n+1} = x_n + \left( {\frac{1}{x_n} -c} \right) \,{x_n^2} = x_n\,\left( 2 - c\,x_n \right).$$

U ovoj formuli nema dijeljenja, i mi smo riješili zadatak, ako postoji interval u kojem možemo birati početnu aproksimaciju tako da ovaj postupak konvergira. Ispitajmo uvjete teorema 24 u ovom slučaju.

$$f'(x) = -{x^{-2}} < 0,\hspace{1cm}f''(x) = 2\,x^{-3} > 0.$$

Neka su sada $a,b$ takvi da je $a<c^{-1}<b.$ Time je uvjet $f(a)f(b)<0$ ispunjen. Iz pozitivnosti $f''$ slijedi da $f'$ raste. No $f'$ ima negativne vrijednosti, pa iz rasta $f'$ slijedi da $\vert f'\vert$ pada. Tako $\vert f'\vert$ prima manju vrijednost u rubnoj točki $b.$

$$\left\vert\frac{f(b)}{f'(b)}\right\vert = b\,\left(b\,c -1 \right),$$

pa $b$ određujemo iz kvadratne nejednadžbe

$$b\,\left(b\,c -1 \right) \leqslant{} b-a.$$

Izlazi da se $b$ mora nalaziti između

$$\frac{1 - \sqrt{1 - a\,c}}{c} ,$$   i$$\hspace{1cm}\frac{1 + \sqrt{1 - a\,c}}{c}.$$

Budući da $a>0$ možemo uzeti po volji malen, slijedi da se za $b$ može uzeti bilo koji broj manji od $2c^{-1}.$ Tako se za početnu aproksimaciju može uzeti bilo koji $x_0$ takav da je

$$0 < x_0 < 2\,c^{-1}.$$

Ako želimo ocijeniti grešku, primijetimo da je Newtonova metoda zapravo metoda iteracije, ako stavimo

$$\varphi(x) = x-\frac{f(x)}{f'(x)}.$$

Tada je

$$\varphi'(x) = 1-\frac{f'(x)}{f'(x)} + \frac{f(x)\,f''(x)}{f'(x)^2} = \frac{f(x)\,f''(x)}{f'(x)^2}.$$

Neka je

$$\left\vert\frac{f(x)\,f''(x)}{f'(x)^2}\right\vert\leqslant L < 1,\hspace{1cm}\forall x\in [a,b].$$

Tada je apriorna ocjena greške dana formulom (3.6), a aposteriorna formulom (3.7).

Primjer 3.9   Naći $L$ u slučaju približnog računanja drugog korijena po formuli (3.10).

Rješenje. Treba naći približnu vrijednost od $\sqrt{c}.$ Neka je

$$n^2 \leqslant{} c < (n+1)^2,\qquad n\in N$$

Prirodno je početnu aproksimaciju uzeti u segmentu $[n,n+1].$ Tada je

$$L \leqslant{} \max{}\left\vert\frac{f(x)\,f''(x)}{f'(x)^2}\right\vert, \quad \forall{}x\in [n,n+1].$$

U ovom slučaju je

$$f(x) = x^2 - c,\quad f'(x) = 2\,x,\quad f''(x) = 2,$$

pa je

$$L \leqslant{} \max_{x\in [n,n+1]}\left\vert\frac{x^2-c}{2\,x^2}\r... ...ight\vert \leqslant{} \frac{1}{2}\,\left\vert 1 - \frac{n}{(n+1)^2}\right\vert.$$

Odavde se vidi da je uvijek $L<\frac{1}{2}.$

Modifikacije Newtonove metode

Whittakerova metoda.

Jedan od problema u Newtonovoj metodi je potreba računanja derivacije u svakom koraku. Najjednostavniji način da se to izbjegne je da se umjesto derivacije $f'(x_n)$ u formulu algoritma uvrsti konstanta

$$x_{n+1}=x_n-\frac{f(x_n)}{m}.$$

Ovisno o broju $m$ ova korekcija Newtonove metode može usporiti konvergenciju. Zato se ona radije koristi nakon što se s nekoliko koraka Newtonovom metodom stiglo blizu rješenja.

Metoda `Regula falsi' ili metoda sekante.

Ideja metode je da se $n+1$-va aproksimacija odredi kao sjecište osi $x$ i sekante kroz točke na grafu funkcije $f,$ čije apscise su prethodne aproksimacije takve, da su vrijednosti funkcije u njima različitog znaka. Dakle, možemo koristiti Newtonovu metodu u kojoj $f'(x_n)$ zamjenimo s koeficijentom sekante

$$\frac{f(x_n)-f(x_{l})}{x_n-x_{l}}.$$

Tako imamo formulu postupka

$$x_{n+1} = x_n - \frac{(x_n-x_{l})\,f(x_n)}{f(x_n)-f(x_{l})} = \frac{x_{l}\,f(x_n) - x_n\,f(x_{l})}{f(x_n) - f(x_{l})}.$$

Pri tom uzimamo $x_0=a,$ $x_1=b,$ pa $x_2$ izračunamo iz formule. Da bismo izračunali $x_3,$ stavimo u formulu $x_2$ umjesto $x_n,$ a umjesto $x_{l}$ stavimo onaj od brojeva $x_0,x_1$ u kojem funkcija prima vrijednost suprotnog znaka od onog u točki $x_2.$ Tako nastavljamo dalje. Da bismo izračunali $x_{n+2},$ trebamo uzeti $x_{n+1},$ i kao $x_l$ uzeti onu prethodnu aproksimaciju u kojoj funkcija ima suprotan znak nego u $x_{n+1}.$

Tako imamo sljedeći algoritam.

Algoritam 6   Stavimo $x_0=a$ i $x_1=b.$ Zatim računamo niz $(x_n), n=2,3,4,\ldots\$ po formuli

$$x_{k+1} = \frac{x_l\,f(x_k) - x_k\,f(x_l)}{f(x_k) - f(x_l)},$$

gdje je $x_l$ ona prethodna aproksimacija u kojoj funkcija ima suprotan znak nego u $x_{k}.$

Mathematica program 5   Metoda sekante

  f[t_]=; (* funkcija *)
a=;
b=;   (* pocetni interval (zadati kao realne brojeve) *)
If[f[a]f[b]<=0,
x=(a f[b]-b f[a])/(f[b]-f[a]);
n=0;
Print[{"   ",n,N[a],N[x],N[b]}];
        While[
        N[f[x]]!=0.,
        Print[N[f[a]],"  ",N[f[x]],"   ",N[f[b]]];
                If[
                f[b]f[x]<0,
                a=x;x=N[(x f[b]-b f[x])/(f[b]-f[x])],
                b=x;x=N[(a f[x]-x f[a])/(f[x]-f[a])]
                ];
        n=n+1;
        If[n>100,Break[]
        ];
Print[{"   ",n,N[a],N[x],N[b]}]
],"Na odabranom segmentu nije ispunjen nuzan uvjet postojanja rjesenja
"]

Metoda uvijek konvergira. Konvergencija je brža nego kod metode polovljenja, ali sporija nego kod Newtonove metode.