MN02

Uwaga: przekonwertowane latex2mediawiki; prawdopodobnie trzeba wprowadzi? poprawki

Rozwiązywanie równań nieliniowych

Możesz zastanawiać się, jak w procesorach implementuje się działania arytmetyczne, na przykład dzielenie. Okazuje się, że dzielenie $b/a$ można zaimplementować korzystając z uprzednio zaimplementowanych operacji dodawania i mnożenia...

W wielu zadaniach, m.in. matematyki stosowanej, spotykamy się z problemem rozwiązania skalarnego równania nieliniowego postaci $f(x)=0$:

• rozwiązywanie równania Keplera

$f(x)\equiv x-ϵ\sin (x)=0$

To równanie jest bardzo ważne w astronomii.

• znajdowanie miejsc zerowych wielomianu:

$f(x)\equiv a_n{x}^n+\dots +a_{1}x+a_0=0$

Bardzo wiele modeli matematycznych wymaga rozwiązania równania z wielomianową nieliniowością. Piękne kwadratury (Gaussa) opierają się na węzłach będących zerami pewnego wielomianu. Wielomiany są bardzo szczególnymi funkcjami i dla nich istnieje szereg specjalizowanych metod znajdowania ich pierwiastków, m.in. metoda Laguerre'a, metoda Bairstow'a (o nich tu nie będziemy mówić), a także zaskakujące metody sprowadzające zadanie poszukiwania miejsc zerowych wielomianu do zupełnie innego zadania matematycznego --- o nich jednak będzie mowa dopiero w wykładzie dotyczącym Uzupełnij: znajdowania wartości własnych macierzy.

• znajdowanie miejsc zerowych trójmianu kwadratowego:

$f(x)\equiv a_2{x}^2+a_{1}x+a_0=0.$

Jest to szczególny, ale oczywiście bardzo ważny (takie równania m.in. trzeba było kiedyś rozwiązywać w artylerii) przypadek poprzedniego zadania. Chociaż wydawać by się mogło, że to umiemy już robić (wszyscy znamy wzory "z deltą") ale --- jak wkrótce się przekonamy --- i tutaj mogą spotkać nas niespodzianki!

• obliczanie pierwiastka kwadratowego z zadanej liczby $a$:

$f(x)\equiv x^2-a=0$

Jeden ze sposobów na implementację funkcji "sqrt()". Szybkie algorytmy wyznaczania pierwiastka kwadratowego były znane już starożytnym. W wykładzie zrozumiemy, dlaczego metoda Herona,

$x_{k+1}=\frac{1}{2}(x_k+\frac{a}{x_k})$

daje bardzo dobre przybliżenie $\sqrt{a}$ już po kilku iteracjach.

• implementacja wyznaczania odwrotności liczby $a$ (bez dzielenia!):

$f(x)\equiv \frac{1}{x}-a=0$

Wciąż spotykane zadanie, np. tak można w praktyce poprawić precyzję funkcji wektorowych stosowanych w niektórych procesorach AMD, zob. \cite{AMD-optimization-guide}. Instrukcja procesora służąca do obliczania odwrotności sekwencji liczb umieszczonych w 128-bitowym rejestrze wektorowym daje wynik z małą precyzją (oczywiście po to, by wykonywała się szybciej!). Jeśli taka dokładność wyniku nie odpowiada nam, możemy ją --- zgodnie z manualem procesora --- poprawić, rozwiązując właśnie takie równanie jak powyżej metodą korzystającą wyłącznie z (wektorowych) operacji mnożenia i dodawania.

Metoda bisekcji

Najprostsza metoda rozwiązywania równania $f(x)=0$.

Metoda bisekcji, czyli połowienia, często stosowana w innych działach informatyki, jest dość naturalną metodą obliczania zer skalarnych funkcji ciągłych określonych na danym przedziale $[a,b]$ i zmieniających znak. Dokładniej, rozpatrzmy klasę funkcji

$F=\{f\in C([a,b]):f(a)\cdot f(b)<0\}.$

Oczywiście, każda funkcja $f\in F$ ma co najmniej jedno zero w $[a,b]$. Startując z przedziału $[a,b]$, w kolejnych krokach metody bisekcji obliczamy informację o wartości $f$ w środku przedziału, co pozwala nam, w zależności od znaku obliczonej wartości, zmniejszyć o połowę przedział, w którym na pewno znajduje się zero funkcji.

<flash>file=</flash><div.thumbcaption>Pierwsze trzy kroki metody bisekcji

Bisekcję realizuje następujący ciąg poleceń, po wykonaniu którego $x$ jest przybliżeniem zera funkcji $f$ z zadaną dokładnością $ϵ$.

Algorytm {Metoda bisekcji}

---

xl = a; xr = b;
x = (a+b)/2;  e = (b-a)/2;
while (e > <math>\displaystyle \epsilon</math>) 
{
	if (f(x)*f(xl) < 0)	
		xr = x;
	else
		xl = x;
	x = (xl+xr)/2; e = e/2;
} 

Z konstrukcji metody łatwo wynika, że po wykonaniu $k$ iteracji (czyli po obliczeniu $k$ wartości funkcji) otrzymujemy $x$, które odległe jest od pewnego rozwiązania $x^{*}$ o co najwyżej

$|x-x^{*}|\le \left(\frac{1}{2}{)}^k\right(\frac{b-a}{2}).$

Metoda bisekcji jest więc zbieżna liniowo z ilorazem $1/2$. Choć ta zbieżność nie jest imponująca, bisekcja ma kilka istotnych zalet. Oprócz jej prostoty, należy podkreślić fakt, że bisekcja jest w pewnym sensie uniwersalna. Jeśli tylko dysponujemy dwoma punktami $a$ i $b$ takimi, że $f$ przyjmuje w nich wartości przeciwnych znaków, to metoda bisekcji z pewnością znajdzie miejsce zerowe funkcji, choćby początkowa długość przedziału $|b-a|$ była bardzo duża: zbieżność metody bisekcji jest globalna. Co ważniejsze, dla zbieżności metody bisekcji wystarcza jedynie ciągłość funkcji. Poza tym możemy łatwo kontrolować błąd bezwzględny aproksymacji miejsca zerowego. Konsekwencją (Uzupelnic: blbis ) jest bowiem następujący wniosek.

Wniosek

Dla znalezienia zera $x^{*}$ z dokładnością $ϵ>0$, wystarczy obliczyć w metodzie bisekcji

$k=k(ϵ)=\lceil {\log }_2\frac{(b-a)}{ϵ}\rceil -1$

wartości funkcji.

Zbieżność metody bisekcji dla ....

Dodajmy jeszcze, że bisekcja minimalizuje błąd najgorszy w klasie $F$ zdefiniowanej Uzupełnij: na początku tej sekcji, wśród wszystkich algorytmów korzystających z określonej liczby obliczeń wartości funkcji, zob. Uzupełnij: uwaga na końcu wykładu.

Uwaga

Metoda bisekcji jest optymalna w następującym sensie. Niech $A:F\to R$ będzie dowolną metodą (algorytmem) aproksymującą zero $x^{*}(f)$ funkcji $f$ z klasy $F$ zdefiniowanej w (Uzupelnic: dfkl ), korzystającą jedynie z obliczeń (informacji o) $f$ w $k$ punktach. Wtedy dla dowolnego $\gamma >0$ istnieje funkcja $f_{\gamma }\in F$ mająca tylko jedno zero $x^{*}(f_{\gamma })$ w $[a,b]$ i taka, że

$|A(f_{\gamma })-x^{*}(f_{\gamma })|\ge \left(\frac{1}{2}{)}^k\right(\frac{b-a}{2})+\gamma .$

Co więcej, można pokazać, że fakt ten zachodzi też w węższej klasie $F_1$ funkcji $f\in F$, które są dowolnie wiele razy różniczkowalne. Oczywiście, nie wyklucza to istnienia metod iteracyjnych (takich jak metoda Newtona), które dla $f\in F_1$ są zbieżne szybciej niż liniowo.

Metoda iteracji prostej Banacha

Zupełnie inne, i jak się okaże --- przy odrobinie sprytu bardzo skuteczne --- podejście do wyznaczania miejsca zerowego jest oparte na metodzie Banacha.

Najpierw nasze równanie nieliniowe

$f(x)=0$

przekształcamy (dobierając odpowiednią funkcję $\varphi$) do równania równoważnego (tzn. mającego te same rozwiązania)

$x=\varphi (x).$

Następnie, startując z pewnego przybliżenia początkowego $x_0$, konstruujemy ciąg kolejnych przybliżeń $x_k$ według wzoru

$x_k=\varphi (x_{k-1}),k\ge 1.$

Twierdzenie Banacha, o zbieżności iteracji prostej

Niech $D_0$ będzie domkniętym 

podzbiorem dziedziny $D$,

${\bar{D}}_0=D_0\subset D,$

w którym $\varphi$ jest odwzorowaniem zwężającym. To znaczy, $\varphi (D_0)\subset D_0$, oraz istnieje stała $0\le L<1$ taka, że

$\Vert \varphi (x)-\varphi (y)\Vert \le L\Vert x-y\Vert ,\mathrm{\forall }x,y\in D_0.$

Wtedy równanie (Uzupelnic: rrw ) ma dokładnie jedno rozwiązanie $x^{*}$, oraz

$x^{*}=\underset{k\to \mathrm{\infty }}{\lim }{x}_k,$

dla dowolnego przybliżenia początkowego $x_0\in D_0$.

Dowód

Wobec

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned \| x_k- x_{k-1}\| &= \|\phi( x_{k-1})-\phi( x_{k-2})\| \,\le\,L\,\| x_{k-1}- x_{k-2}\| \\ &\le &\cdots\;\le\;L^{k-1}\| x_1- x_0\|, \endaligned}

dla $k\ge s$ mamy

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned \| x_k- x_s\| &\le & \sum_{j=s+1}^k\| x_j- x_{j-1}\| \,\le\,\sum_{j=s+1}^k L^{j-1}\| x_1- x_0\| \\ &= L^s(1+L+\cdots+L^{k-s-1})\| x_1- x_0\| \,\le\,\frac{L^s}{1-L}\| x_1- x_0\|. \endaligned}

Ciąg $\{x_k{\}}_k$ jest więc ciągiem Cauchy'ego. Stąd istnieje granica $\overrightarrow{\alpha }=\underset{k\to \mathrm{\infty }}{\lim }{x}_k$, która należy do $D_0$, wobec domkniętości tego zbioru. Ponieważ lipschitzowskość $\varphi$ implikuje jej ciągłość, mamy też

$\varphi (\overrightarrow{\alpha })=\varphi \left(\underset{k\to \mathrm{\infty }}{\lim }{x}_k\right)=\underset{k\to \mathrm{\infty }}{\lim }\varphi (x_k)=\underset{k\to \mathrm{\infty }}{\lim }{x}_k=\overrightarrow{\alpha },$

tzn. $\overrightarrow{\alpha }$ jest punktem stałym odwzorowania $\varphi$. Dla jednoznaczności zauważmy, że jeśliby istniał drugi, różny od $\overrightarrow{\alpha }$, punkt stały $\overrightarrow{\beta }$, to mielibyśmy

$\Vert \overrightarrow{\alpha }-\overrightarrow{\beta }\Vert =\Vert \varphi (\overrightarrow{\alpha })-\varphi (\overrightarrow{\beta })\Vert \le L\Vert \overrightarrow{\alpha }-\overrightarrow{\beta }\Vert .$

Stąd $1<L$, co jest sprzeczne z założeniem, że

$\varphi$ jest zwężająca.

Z powyższych rozważań otrzymujemy natychmiastowy wniosek dotyczący zbieżności iteracji prostych.

Wniosek

Przy założeniach Uzupełnij: twierdzenia Banacha, metoda iteracji prostych jest zbieżna co najmniej liniowo z ilorazem $L$, tzn.

$\Vert x_k-x^{*}\Vert \le L^k\Vert x_0-x^{*}\Vert .$

Przykład

Dla ilustracji, rozpatrzmy natępujące proste równanie skalarne:

$x=\cos (x),\text{dla}x\in D=R.$

W tym przypadku $\varphi (x)=\cos (x)$. Zauważamy, że w przedziale $[0,1]$ funkcja $\varphi$ jest zwężająca ze stałą

$L=\underset{0\le x\le 1}{\max }|{\cos }^{\prime }(x)|=\sin (1)<1.$

Stąd istnieje dokładnie jedno rozwiązanie naszego równania w przedziale $[0,1]$. Rozwiązanie to może być aproksymowane z dowolnie małym błędem przy pomocy iteracji prostych, startując z dowolnego przybliżenia początkowego $x_0\in [0,1]$.

Zaletą iteracji prostych jest fakt, że zbieżność nie zależy od wymiaru $n$ zadania, ale tylko od stałej Lipschitza $L$ (jednak w praktyce czasem sama stała Lipschitza może zależeć od wymiaru zadania...). Metoda Banacha ma szczególne zastosowanie w przypadku, gdy funkcja $\varphi$ jest zwężająca na całym zbiorze $D$, tzn. $D_0=D$. Jeśli ponadto $D$ ma skończoną średnicę $\text{diam}(D)$, to dla osiągnięcia $ϵ$-aproksymacji zera funkcji $f$ wystarczy wykonać

$k=k(ϵ)=\lceil \frac{\log (\Vert x_0-x^{*}\Vert /ϵ)}{\log (1/L)}\rceil =\lceil \frac{\log (\text{diam}(D)/ϵ)}{\log (1/L)}\rceil$

iteracji, niezależnie od $x_0$. Metody zbieżne dla dowolnego przybliżenia początkowego, nazywamy zbieżnymi globalnie. Obie przedstawione dotychczas metody: bisekcji i Banacha, przy rozsądnych założeniach, są zbieżne globalnie.

Okazuje się, że metoda iteracji prostej może być --- w bardzo szczególnych przypadkach --- zbieżna szybciej niż liniowo. Z taką sytuacją będziemy mieli, gdy korzystać będziemy z metody Newtona.

Metoda Newtona

Zarówno metoda Banacha jak i bisekcja są zbieżnie liniowo, co w praktyce może okazać się zbieżnością dość powolną (np. dla metody zbieżnej liniowo z ilorazem $\frac{1}{2}$, dopiero co 5 iteracji, dostajemy kolejną dokładną cyfrę wyniku). Wykorzystując więcej informacji o funkcji $f$, której miejsca zerowego poszukujemy, możemy istotnie przyspieszyć zbieżność metody. Ceną, jaką przyjdzie nam zapłacić, będzie utrata globalnej zbieżności.

Metoda Newtona i jej podobne należą do grupy metod zbieżnych lokalnie. Znaczy to, że zbieżność ciągu $\{x_k{\}}_k$ do zera danej funkcji $f$ jest zapewniona jedynie wtedy, gdy przybliżenia początkowe zostały wybrane dostatecznie blisko $x^{*}$.

W dalszych rozważaniach będziemy zakładać dla uproszczenia, że dziedzina $D=R$.

Idea metody Newtona opiera się na popularnym wśród inżynierów pomyśle linearyzacji: zamiast szukać miejsca zerowego skomplikowanej $f$, przybliżmy ją linią prostą, a dla niej już umiemy znaleźć miejsce zerowe!

Startując z pewnego przybliżenia początkowego $x_0$, w kolejnych krokach metody, $k$-te przybliżenie $x_k$ jest punktem przecięcia stycznej do wykresu $f$ w punkcie $x_{k-1}$. Ponieważ równanie stycznej wynosi $y(x)=f(x_{k-1})+f^{\prime }(x_{k-1})(x-x_{k-1})$, otrzymujemy wzór

$x_k=x_{k-1}-\frac{f(x_{k-1})}{f^{\prime }(x_{k-1})}.$

Oczywiście, aby metoda Newtona była dobrze zdefiniowana, musimy założyć, że $f^{\prime }(x_{k-1})$ istnieje i nie jest zerem.

Zauważmy, że metodę Newtona można traktować jako szczególny przypadek iteracji prostych, gdzie

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

Widać też, że nie jest ona zbieżna globalnie.

Nawet jeśli pochodna w $x_{k-1}$ się nie zeruje, ciąg $\{x_k{\}}_k$ może nie zbiegać do zera funkcji $f$. Okazuje się jednak, że jeśli wystartujemy dostatecznie blisko rozwiązania $x^{*}$, to metoda Newtona jest zbieżna. Dokładniej, załóżmy najpierw, że $f(x^{*})=0$ oraz

$f^{\prime }(x^{*})\ne 0.$

Ponadto załóżmy, że $f$ jest dwukrotnie różniczkowalna w sposób ciągły, $f\in C^2(D)$. Rozwijając $\varphi$ w szereg Taylora w punkcie $x^{*}$ otrzymujemy

$x_k-x^{*}=\varphi (x_{k-1})-\varphi (x^{*})=(x_{k-1}-x^{*}){\varphi }^{\prime }(x^{*})+(x_{k-1}-x^{*}{)}^2{\varphi }^{″}({\xi }_k)/2,$

gdzie $\min (x^{*},x_{k-1})\le {\xi }_k\le \max (x^{*},x_{k-1})$. Wobec tego, że ${\varphi }^{\prime }(x^{*})=f(x)f^{″}(x)/(f^{\prime }(x){)}^2=0$ i ${\varphi }^{″}({\xi }_k)=f^{″}({\xi }_k)/f^{\prime }({\xi }_k)$, mamy

$x_k-x^{*}=(x_{k-1}-x^{*}{)}^2\frac{f^{″}({\xi }_k)}{2f^{\prime }({\xi }_k)}.$

Zdefiniujmy liczbę

$R_f=\underset{r\ge 0}{\sup }\underset{\{x:|x-x^{*}|\le r\}}{\sup }\left|\frac{2(x-x^{*})f^{″}(x)}{f^{\prime }(x)}\right|<1.$

Oczywiście $R_f>0$. Dla $x_{k-1}$ spełniającego $|x_{k-1}-x^{*}|\le R<R_f$, mamy z poprzedniej równości (Uzupelnic: nrpdst )

$|x_k-x^{*}|\le q|x_{k-1}-x^{*}|,$

gdzie $q<1$ i $q$ zależy tylko od $R$.

Niech teraz $x^{*}$ będzie zerem $m$-krotnym,

$f(x^{*})=f^{\prime }(x^{*})=\cdots =f^{(m-1)}(x^{*})=0\ne f^{(m)}(x^{*}),$

gdzie $m\ge 2$, oraz niech $f$ będzie $m$-krotnie różniczkowalna w sposób ciągły. Wtedy

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned x_k-x^* &= (x_{k-1}-x^*)\,-\,\frac{(x_{k-1}-x^*)^m \frac{f^{(m)} (\eta_k^{(1)})}{m!}}{(x_{k-1}-x^*)^{m-1} \frac{f^{(m-1)}(\eta_k^{(2)})}{(m-1)!}} \nonumber \\ &= (x_{k-1}-x^*)\left(1-\frac 1m\frac {f^{(m)}(\eta_k^{(1)})}{f^{(m)}(\eta_k^{(2)})} \right) \nonumber \\ &\approx & (x_{k-1}-x^*)\Big( 1-\frac 1m\Big), \endaligned}

o ile $x_{k-1}$ jest "blisko" $x^{*}$.

Metoda Newtona jest więc zbieżna lokalnie. Z (Uzupelnic: nrpdst ) i (Uzupelnic: nrtp ) można też wywnioskować, jaki jest charakter zbieżności metody Newtona. Dla zera jednokrotnego $x^{*}$ oraz $f^{″}(x^{*})\ne 0$ mamy bowiem

$(x_k-x^{*})\approx (x-x_{k-1}{)}^2\frac{f^{″}(x^{*})}{2f^{\prime }(x^{*})}.$

Mówimy, że zbieżność jest kwadratowa. Jeśli zaś $f^{″}(x^{*})=0$ to zbieżnośc jest nawet szybsza. Z kolei dla zera $m$-krotnego zbieżność jest liniowa z ilorazem $(1-\frac{1}{m})$.

Metoda Newtona jest pierwszą poznaną tutaj metodą iteracyjną, która jest (dla zer jednokrotnych) zbieżna szybciej niż liniowo. Dla takich metod wprowadza się pojęcie wykładnika zbieżności, który jest zdefiniowany następująco.

Powiemy, że metoda iteracyjna $\varphi$ jest w klasie funkcji $F$ rzędu co najmniej $p\ge 1$, gdy spełniony jest następujący warunek. Niech $f\in F$ i $f(x^{*})=0$. Wtedy istnieje stała $C<\mathrm{\infty }$ taka, że dla dowolnych przybliżeń początkowych $x_0,\dots ,x_{s-1}$ dostatecznie bliskich $x^{*}$, kolejne przybliżenia $x_k=\varphi (x_{k-1},\dots ,x_{k-s})$ generowane tą metodą spełniają

$|x_k-x^{*}|\le C|x_{k-1}-x^{*}{|}^p.$

Ponadto, jeśli $p=1$ to dodatkowo żąda się, aby $C<1$.

Definicja

Wykładnikiem zbieżności metody iteracyjnej $\varphi$ w klasie $F$ nazywamy liczbę $p^{*}$ zdefiniowaną równością

$p^{*}=\sup \{p\ge 1:\varphi \text{ jest rzędu co najmniej }p\}.$

Możemy teraz sformułować następujące twierdzenie, które natychmiast wynika z poprzednich rozważań.

Twierdzenie

Wykładnik zbieżności metody Newtona (stycznych) wynosi $p^{*}=2$ w klasie funkcji o zerach jednokrotnych, oraz $p^{*}=1$ w klasie funkcji o zerach wielokrotnych.

Zbieżność metody Newtona na tle metody bisekcji

Zbieżność metody Newtona dla zer wielokrotnych

Metoda siecznych

Inną znaną i często używaną metodą iteracyjną, opartą na podobnym pomyśle linearyzacyjnych co metoda Newtona jest metoda siecznych, w której zamiast przybliżenia wykresu $f$ przez styczną, stosuje się

przybliżenie sieczną.

Metoda ta wykorzystuje więc do konstrukcji $x_k$ przybliżenia $x_{k-1}$ i $x_{k-2}$. Musimy również wybrać dwa różne punkty startowe $x_0$ i $x_1$. Ponieważ prosta interpolująca $f$ w $x_{k-1}$ i $x_{k-2}$ ma wzór

$y(x)=\frac{x-x_{k-2}}{x_{k-1}-x_{k-2}}f(x_{k-1})+\frac{x-x_{k-1}}{x_{k-2}-x_{k-1}}f(x_{k-2}),$

otrzymujemy

$x_k=x_{k-1}-\frac{x_{k-1}-x_{k-2}}{f(x_{k-1})-f(x_{k-2})}f(x_{k-1}).$

Zauważmy, że jeśli $x_{k-1}$ i $x_{k-2}$ są blisko siebie, to $x_k$ jest podobny do tego z metody Newtona, bowiem wtedy iloraz różnicowy

$\frac{f(x_{k-1})-f(x_{k-2})}{x_{k-1}-x_{k-2}}\approx f^{\prime }(x_{k-1}).$

Nie wystarcza to jednak, aby osiągnąć zbieżność z wykładnikiem $2$. Dokładniej, można pokazać, że wykładnik zbieżności metody siecznych dla zer jednokrotnych dostatecznie gładkich funkcji wynosi $p^{*}=\frac{1+\sqrt{5}}{2}=1.618\dots$. Jako wariant metody Newtona, metoda siecznych jest również zbieżna lokalnie.

Niewątpliwą zaletą metody siecznych jest jednak to, że nie wymaga ona obliczania pochodnej funkcji (co w praktyce jest często bardzo trudne, a niekiedy nawet niemożliwe), a tylko jej wartości. Jest to również istotne w pakietach numerycznych, gdzie czasem nie chcemy wymagać od użytkownika czegokolwiek ponad podanie funkcji i przybliżonej lokalizacji miejsca zerowego.

Ponadto, często zdarza się, że wyznaczenie wartości pochodnej, $f^{\prime }(x_k)$, jest tak samo, albo i bardziej kosztowne od wyznaczenia wartości $f(x_k)$. W takim wypadku okazuje się, że metoda stycznych --- choć wolniej zbieżna niż metoda stycznych --- dzięki temu, że jej iteracja wymaga jedynie wyznaczenia jednej wartości $f$, jest bardziej efektywna od metody Newtona: koszt osiągnięcia zadanej dokładności jest w takim przypadku mniejszy od analogicznego kosztu dla metody Newtona.

Jednak, gdy żądane przez użytkownika dokładności są bardzo wielkie, a sama funkcja "złośliwa", metoda siecznych może cierpieć z powodu redukcji cyfr przy odejmowaniu.

Zbieżność metody siecznych na tle metody Newtona

Metoda Brenta

Naturalnie, uważny student zaczyna zadawać sobie pytanie, czy nie można w jakiś sposób połączyć globalnej zbieżności metody bisekcji z szybką zbieżnością metody siecznych tak, by uzyskać metodę zbieżną globalnie, a jednocześnie istotnie szybciej niż liniowo. (Wariant odwrotny: opracowanie metody łączącej wolną zbieżność bisekcji z lokalną zbieżnością siecznych, pozostawiamy studentom gorszych uczelni).

Okazuje się, że można to zrobić, wprowadzając metodę opartą na trzech punktach lokalizujących miejsce zerowe: dwóch odcinających zero tak jak w metodzie bisekcji i trzecim, konstruowanym jak np. w metodzie stycznych. W kolejnej iteracji konstruujemy wymieniamy jeden z punktów albo wedle metody siecznych (i wtedy zapewne szybciej zbliżamy się do zera), albo robiąc bisekcję (aby zagwarantować sobie, że w wiadomym przedziale miejsce zerowe rzeczywiście sie znajduje).

Ten prosty pomysł metody hybrydowej wymaga jednak subtelnego dopracowania, co zostało zrobione w 1973 roku przez Richarda Brenta, który twórczo rozwinął wcześniejsze idee Dekkera, van Wijngaardena i Dijkstry.

Funkcja MATLABa (i Octave'a) fzero implementuje metodę Brenta. Ciekawostką jest, że autorem implementacji w Octave jest ówczesny student matematyki na Uniwersytecie Warszawskim, Łukasz Bodzon.

Porównanie zbieżności różnych metod dla równań nieliniowych: zero jednokrotne

Porównanie zbieżności różnych metod dla równań nieliniowych: zero wielokrotne

Ciekawostki

O ile metoda Brenta jest oczywiście jego autorstwa, o tyle przypisywanie metody stycznych Newtonowi jest pewną przesadą. Metodę Newtona taką jaką znamy (z pochodną w mianowniku) zaproponował w 1740 roku Simpson (ten od kwadratury), a więc kilknaście lat po śmierci Newtona. Żeby było jeszcze zabawniej, odkrywcą metody siecznych zdaje się być... Newton! Więcej na ten temat przeczytasz w artykule T.Ypma w SIAM Review 37, 1995.