MN02: Różnice pomiędzy wersjami

Z Studia Informatyczne
Przejdź do nawigacjiPrzejdź do wyszukiwania
Aneczka (dyskusja | edycje)
Przykry (dyskusja | edycje)
mNie podano opisu zmian
Linia 1: Linia 1:


==Rozwiązywanie równań nieliniowych==
=Równania nieliniowe=
 
Możesz zastanawiać się, jak w procesorach implementuje się działania
arytmetyczne, na przykład dzielenie. Okazuje się, że dzielenie <math>\displaystyle b/a</math> można
zaimplementować korzystając z uprzednio zaimplementowanych operacji
dodawania i mnożenia...


W wielu zadaniach, m.in. matematyki stosowanej, spotykamy się z problemem
W wielu zadaniach, m.in. matematyki stosowanej, spotykamy się z problemem
Linia 20: Linia 15:


Bardzo wiele modeli matematycznych wymaga rozwiązania równania z wielomianową
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
nieliniowością. [[sec:kwadratury|Dodaj link: Piękne kwadratury]] (Gaussa) opierają się na węzłach będących
zerami pewnego wielomianu. Wielomiany są bardzo szczególnymi funkcjami i dla
zerami pewnego wielomianu. Wielomiany są bardzo szczególnymi funkcjami i dla
nich istnieje szereg specjalizowanych metod znajdowania ich pierwiastków, m.in.
nich istnieje szereg specjalizowanych metod znajdowania ich pierwiastków, m.in.
Linia 26: Linia 21:
zaskakujące metody sprowadzające zadanie poszukiwania miejsc zerowych wielomianu
zaskakujące metody sprowadzające zadanie poszukiwania miejsc zerowych wielomianu
do zupełnie innego zadania matematycznego  --- o nich jednak
do zupełnie innego zadania matematycznego  --- o nich jednak
będzie mowa dopiero w wykładzie dotyczącym [[sec:eigenvalue|Uzupełnij: znajdowania wartości własnych
będzie mowa dopiero w wykładzie dotyczącym [[sec:eigenvalue|Dodaj link: znajdowania wartości własnych
macierzy]].
macierzy]].
* znajdowanie miejsc zerowych trójmianu kwadratowego:  
* znajdowanie miejsc zerowych trójmianu kwadratowego:  
Linia 43: Linia 38:
x^2 - a = 0</math></center>
x^2 - a = 0</math></center>


Jeden ze sposobów na implementację funkcji "<code>sqrt()</code>". Szybkie algorytmy
czyli sposób na implementację funkcji "<code>sqrt()</code>". Szybkie algorytmy
wyznaczania pierwiastka kwadratowego były znane już starożytnym. W wykładzie
wyznaczania pierwiastka kwadratowego były znane już starożytnym. W wykładzie
zrozumiemy, dlaczego metoda Herona,
zrozumiemy, dlaczego metoda Herona,
[[grafika:Heron.jpg|thumb|right|| Heron<br>  [[Biografia Heron|Zobacz biografi�]]]]
<center><math>\displaystyle  
<center><math>\displaystyle  
x_{k+1} = \frac{1}{2}\left(x_k + \frac{a}{x_k}\right)
x_{k+1} = \frac{1}{2}\left(x_k + \frac{a}{x_k}\right)
Linia 51: Linia 49:


daje bardzo dobre przybliżenie <math>\displaystyle \sqrt{a}</math> już po kilku iteracjach.
daje bardzo dobre przybliżenie <math>\displaystyle \sqrt{a}</math> już po kilku iteracjach.
* implementacja wyznaczania odwrotności liczby <math>\displaystyle a</math> (''bez'' dzielenia!):  
* implementacja wyznaczania odwrotności liczby <math>\displaystyle a</math> (<strong>bez</strong> dzielenia!):  


<center><math>\displaystyle f(x) \equiv
<center><math>\displaystyle f(x) \equiv
\frac{1}{x} - a = 0</math></center>
\frac{1}{x} - a = 0</math></center>


Wciąż spotykane zadanie, np. tak można w praktyce poprawić precyzję funkcji
Wciąż spotykane zadanie, np. tak można w praktyce poprawić precyzję
wektorowych stosowanych w niektórych procesorach AMD, zob.
[http://www.amd.com/us-en/assets/content_type/white_papers_and_tech_docs/21928.pdf  funkcji
'''(Uzupełnij źródło: AMD-optimization-guide)'''. Instrukcja procesora służąca do obliczania
wektorowych stosowanych w niektórych procesorach AMD]. Instrukcja procesora służąca do obliczania
odwrotności sekwencji liczb umieszczonych w 128-bitowym rejestrze wektorowym
odwrotności sekwencji liczb umieszczonych w 128-bitowym rejestrze wektorowym
daje wynik z małą precyzją (oczywiście po to, by wykonywała się szybciej!).
daje wynik z małą precyzją (oczywiście po to, by wykonywała się szybciej!).
Linia 65: Linia 63:
korzystającą wyłącznie z (wektorowych) operacji mnożenia i dodawania.
korzystającą wyłącznie z (wektorowych) operacji mnożenia i dodawania.


===Metoda bisekcji===
==Bisekcja==


Najprostsza metoda rozwiązywania równania <math>\displaystyle f(x) = 0</math>.
<strong>Metoda bisekcji</strong>, czyli <strong>połowienia</strong>, często stosowana w innych
 
''Metoda bisekcji'', czyli ''połowienia'', często stosowana w innych
działach informatyki, jest dość  
działach informatyki, jest dość  
naturalną metodą obliczania zer skalarnych funkcji  
naturalną metodą obliczania zer skalarnych funkcji  
Linia 87: Linia 83:
o połowę przedział, w którym na pewno znajduje się  
o połowę przedział, w którym na pewno znajduje się  
zero funkcji.  
zero funkcji.  
<div class="thumb tright"><div><flash>file=bisekcja.swf</flash><div.thumbcaption>Pierwsze trzy kroki metody bisekcji</div></div></div>


Bisekcję realizuje następujący ciąg  
Bisekcję realizuje następujący ciąg  
Linia 94: Linia 88:
zera funkcji <math>\displaystyle f</math> z zadaną dokładnością <math>\displaystyle \epsilon</math>.  
zera funkcji <math>\displaystyle f</math> z zadaną dokładnością <math>\displaystyle \epsilon</math>.  


{{algorytm|{Metoda bisekcji}||
{{algorytm|Metoda bisekcji||
<pre>
<pre>


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


Metoda bisekcji jest więc zbieżna ''liniowo'' z  
Metoda bisekcji jest więc zbieżna <strong>liniowo</strong> z  
ilorazem <math>\displaystyle 1/2</math>. Choć ta zbieżność nie jest  
ilorazem <math>\displaystyle 1/2</math>. Choć ta zbieżność nie jest  
imponująca, bisekcja ma kilka istotnych zalet. Oprócz  
imponująca, bisekcja ma kilka istotnych zalet. Oprócz  
Linia 126: Linia 120:
takimi, że <math>\displaystyle f</math> przyjmuje w nich wartości przeciwnych znaków, to metoda bisekcji
takimi, że <math>\displaystyle f</math> przyjmuje w nich wartości przeciwnych znaków, to metoda bisekcji
z pewnością znajdzie miejsce zerowe funkcji, choćby początkowa długość
z pewnością znajdzie miejsce zerowe funkcji, choćby początkowa długość
przedziału <math>\displaystyle |b-a|</math> była bardzo duża: zbieżność metody bisekcji jest ''globalna''. Co ważniejsze, dla zbieżności metody bisekcji
przedziału <math>\displaystyle |b-a|</math> była bardzo duża: zbieżność metody bisekcji jest <strong>globalna</strong>. Co ważniejsze, dla zbieżności metody bisekcji
wystarcza jedynie ''ciągłość'' funkcji. Poza tym  
wystarcza jedynie <strong>ciągłość</strong> funkcji. Poza tym  
możemy łatwo kontrolować ''błąd bezwzględny aproksymacji miejsca zerowego''. Konsekwencją  
możemy łatwo kontrolować <strong>błąd bezwzględny aproksymacji miejsca zerowego</strong>. Konsekwencją  
powyższego oszacowania błędu jest bowiem następujący wniosek.  
powyższego oszacowania błędu jest bowiem następujący wniosek.  


Linia 142: Linia 136:
}}
}}


[[Image:MNXXX.png|thumb|300px|Zbieżność metody bisekcji dla ....]]
==Iteracja prosta Banacha==
 
===Metoda iteracji prostej Banacha===


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


Najpierw nasze równanie nieliniowe  
Najpierw nasze równanie nieliniowe  
Linia 162: Linia 154:
</math></center>
</math></center>


Następnie, startując z pewnego przybliżenia  
Następnie, startując z pewnego przybliżenia  
początkowego <math>\displaystyle  x_0</math>, konstruujemy ciąg kolejnych  
początkowego <math>\displaystyle  x_0</math>, konstruujemy ciąg kolejnych  
przybliżeń <math>\displaystyle  x_k</math> według wzoru  
przybliżeń <math>\displaystyle  x_k</math> według wzoru  


<center><math>\displaystyle x_k\,=\,\phi( x_{k-1}),\qquad k\ge 1.
<center><math>\displaystyle x_k\,=\,\phi( x_{k-1}),\qquad k\ge 1.
</math></center>
</math></center>
[[grafika:Banach.jpg|thumb|right||Stefan Banach<br>  [[Biografia Banach|Zobacz biografi�]]]]


{{twierdzenie|Banacha, o zbieżności iteracji prostej||
{{twierdzenie|Banacha, o zbieżności iteracji prostej||
Linia 221: Linia 215:
Ciąg <math>\displaystyle \{ x_k\}_k</math> jest więc ciągiem Cauchy'ego.  
Ciąg <math>\displaystyle \{ x_k\}_k</math> jest więc ciągiem Cauchy'ego.  
Stąd istnieje granica  
Stąd istnieje granica  
<math>\displaystyle \vec\alpha=\lim_{k\to\infty} x_k</math>, która należy do  
<math>\displaystyle \alpha=\lim_{k\to\infty} x_k</math>, która należy do  
<math>\displaystyle D_0</math>, wobec domkniętości tego zbioru. Ponieważ  
<math>\displaystyle D_0</math>, wobec domkniętości tego zbioru. Ponieważ  
lipschitzowskość <math>\displaystyle \phi</math> implikuje jej ciągłość,  
lipschitzowskość <math>\displaystyle \phi</math> implikuje jej ciągłość,  
mamy też   
mamy też   


<center><math>\displaystyle \phi(\vec\alpha)\,=\,\phi\Big(\lim_{k\to\infty} x_k\Big)  
<center><math>\displaystyle \phi(\alpha)\,=\,\phi\Big(\lim_{k\to\infty} x_k\Big)  
   \,=\,\lim_{k\to\infty}\phi( x_k)
   \,=\,\lim_{k\to\infty}\phi( x_k)
   \,=\,\lim_{k\to\infty} x_k\,=\,\vec\alpha,
   \,=\,\lim_{k\to\infty} x_k\,=\,\alpha,
</math></center>
</math></center>


tzn. <math>\displaystyle \vec\alpha</math> jest punktem stałym odwzorowania <math>\displaystyle \phi</math>.  
tzn. <math>\displaystyle \alpha</math> jest punktem stałym odwzorowania <math>\displaystyle \phi</math>.  
Dla jednoznaczności zauważmy, że jeśliby istniał  
Dla jednoznaczności zauważmy, że jeśliby istniał  
drugi, różny od <math>\displaystyle \vec\alpha</math>, punkt stały <math>\displaystyle \vec\beta</math>,  
drugi, różny od <math>\displaystyle \alpha</math>, punkt stały <math>\displaystyle \beta</math>,  
to mielibyśmy  
to mielibyśmy  


<center><math>\displaystyle \|\vec\alpha-\vec\beta\|\,=\,
<center><math>\displaystyle \|\alpha-\beta\|\,=\,
   \|\phi(\vec\alpha)-\phi(\vec\beta)\|
   \|\phi(\alpha)-\phi(\beta)\|
   \,\le\,L\,\|\vec\alpha-\vec\beta\|.  
   \,\le\,L\,\|\alpha-\beta\|.  
</math></center>
</math></center>


Linia 248: Linia 242:


{{wniosek|||
{{wniosek|||
Przy założeniach [[twit|Uzupe�nij: twierdzenia Banacha]],  
Przy założeniach [[twit|Dodaj link: twierdzenia Banacha]],  
metoda iteracji prostych jest zbieżna co  
metoda iteracji prostych jest zbieżna co  
najmniej liniowo z ilorazem <math>\displaystyle L</math>, tzn.
najmniej liniowo z ilorazem <math>\displaystyle L</math>, tzn.
Linia 257: Linia 251:
}}
}}


{{przyklad|||
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="font-variant:small-caps;">Przykład</span>
<div class="solution">
 
Dla ilustracji, rozpatrzmy natępujące proste  
Dla ilustracji, rozpatrzmy natępujące proste  
równanie skalarne:  
równanie skalarne:  
Linia 277: Linia 274:
iteracji prostych, startując z dowolnego przybliżenia  
iteracji prostych, startując z dowolnego przybliżenia  
początkowego <math>\displaystyle  x_0\in [0,1]</math>.  
początkowego <math>\displaystyle  x_0\in [0,1]</math>.  
}}
</div></div>


Zaletą iteracji prostych jest fakt, że zbieżność  
Zaletą iteracji prostych jest fakt, że zbieżność  
Linia 297: Linia 294:
iteracji, niezależnie od <math>\displaystyle x_0</math>. Metody zbieżne dla  
iteracji, niezależnie od <math>\displaystyle x_0</math>. Metody zbieżne dla  
dowolnego przybliżenia początkowego, nazywamy  
dowolnego przybliżenia początkowego, nazywamy  
''zbieżnymi globalnie''. Obie przedstawione dotychczas metody: bisekcji i
<strong>zbieżnymi globalnie</strong>. Obie przedstawione dotychczas metody: bisekcji i
Banacha, przy rozsądnych
Banacha, przy rozsądnych
założeniach, są zbieżne globalnie.
założeniach, są zbieżne globalnie.
Linia 305: Linia 302:
gdy korzystać będziemy z metody Newtona.
gdy korzystać będziemy z metody Newtona.


===Metoda Newtona===
==Metoda Newtona==


Zarówno metoda Banacha jak i bisekcja są zbieżnie liniowo, co w praktyce może
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
okazać się zbieżnością dość powolną (np. dla metody zbieżnej liniowo z ilorazem
<math>\displaystyle \frac{1}{2}</math>, dopiero co 5 iteracji, dostajemy kolejną
<math>\displaystyle \frac{1}{2}</math>, dopiero po piątej iteracji, dostajemy kolejną
dokładną cyfrę wyniku). Wykorzystując więcej informacji o funkcji <math>\displaystyle f</math>, której
dokładną cyfrę wyniku). Wykorzystując więcej informacji o funkcji <math>\displaystyle f</math>, której
miejsca zerowego poszukujemy, możemy istotnie przyspieszyć zbieżność metody.
miejsca zerowego poszukujemy, możemy istotnie przyspieszyć zbieżność metody.
Linia 315: Linia 312:


Metoda Newtona i jej podobne należą do  
Metoda Newtona i jej podobne należą do  
grupy metod ''zbieżnych lokalnie''. Znaczy to, że  
grupy metod <strong>zbieżnych lokalnie</strong>. Znaczy to, że  
zbieżność ciągu <math>\displaystyle \{x_k\}_k</math> do zera danej funkcji <math>\displaystyle f</math>  
zbieżność ciągu <math>\displaystyle \{x_k\}_k</math> do zera danej funkcji <math>\displaystyle f</math>  
jest zapewniona jedynie wtedy, gdy przybliżenia początkowe  
jest zapewniona jedynie wtedy, gdy przybliżenia początkowe  
Linia 323: Linia 320:
uproszczenia, że dziedzina <math>\displaystyle D=R</math>.  
uproszczenia, że dziedzina <math>\displaystyle D=R</math>.  


Idea metody Newtona opiera się na popularnym wśród inżynierów pomyśle ''linearyzacji'': zamiast szukać miejsca zerowego skomplikowanej <math>\displaystyle f</math>, przybliżmy ją
Idea metody Newtona opiera się na popularnym wśród inżynierów pomyśle <strong>linearyzacji</strong>: zamiast szukać miejsca zerowego skomplikowanej <math>\displaystyle f</math>, przybliżmy ją
linią prostą, a dla niej już umiemy znaleźć miejsce zerowe!  
linią prostą, a dla niej już umiemy znaleźć miejsce zerowe!  
[[grafika:Newton.jpg|thumb|right||Isaac Newton<br> 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. [[Biografia Newton|Zobacz biografi�]]]]


Startując z pewnego przybliżenia  
Startując z pewnego przybliżenia  
Linia 339: Linia 343:
musimy założyć, że <math>\displaystyle f'(x_{k-1})</math> istnieje i nie  
musimy założyć, że <math>\displaystyle f'(x_{k-1})</math> istnieje i nie  
jest zerem.
jest zerem.
<div class="thumb tright"><div><flash>file=Newton.swf</flash><div.thumbcaption>Postęp iteracji Newtona</div></div></div>


Zauważmy, że metodę Newtona można traktować jako  
Zauważmy, że metodę Newtona można traktować jako  
Linia 413: Linia 419:
jednokrotnego <math>\displaystyle x^*</math> oraz <math>\displaystyle f''(x^*)\ne 0</math> mamy bowiem  
jednokrotnego <math>\displaystyle x^*</math> oraz <math>\displaystyle f''(x^*)\ne 0</math> mamy bowiem  


<center><math>\displaystyle (x_k-x^*)\,\approx\,(x-x_{k-1})^2\frac{f''(x^*)}{2f'(x^*)}.
<center><math>\displaystyle |x_k-x^*|\, \approx \,|x-x_{k-1}|^2 \frac{|f''(x^*)|}{2|f'(x^*)|}.
</math></center>
</math></center>


Mówimy, że zbieżność jest ''kwadratowa''. Jeśli zaś
Mówimy, że zbieżność metody Newtona, gdy <math>\displaystyle f(x^*)\neq 0</math> jest <strong>kwadratowa</strong>.  
<math>\displaystyle f''(x^*)=0</math> to zbieżnośc jest nawet szybsza. Z kolei dla  
 
zera <math>\displaystyle m</math>-krotnego zbieżność jest liniowa z ilorazem  
{{stwierdzenie|||
<math>\displaystyle (1-\frac{1}{m})</math>.  
Jeśli <math>\displaystyle f(x^*)\neq 0</math> oraz
<math>\displaystyle f''(x^*)=0</math> to zbieżność jest nawet szybsza. Z kolei dla  
zera <math>\displaystyle m</math>-krotnego (tzn. <math>\displaystyle f(x^*) = f'(x^*)= \ldots f^{(m)}(x^*)= 0</math>, <math>\displaystyle m>1</math>)
zbieżność jest liniowa z ilorazem <math>\displaystyle (1-\frac{1}{m})</math>.  
}}


Metoda Newtona jest pierwszą poznaną tutaj metodą  
Metoda Newtona jest pierwszą poznaną tutaj metodą  
iteracyjną, która jest (dla zer jednokrotnych) zbieżna  
iteracyjną, która jest (dla zer jednokrotnych) zbieżna  
szybciej niż liniowo. Dla takich metod wprowadza się  
szybciej niż liniowo. Dla takich metod wprowadza się  
pojęcie ''wykładnika zbieżności'', który jest  
pojęcie <strong>wykładnika zbieżności</strong>, który jest  
zdefiniowany następująco.  
zdefiniowany następująco.  
[[Image:MNstycznebisekcja.png|frame|400px|center|Porównanie zbieżności metody bisekcji i stycznych
dla równania <math>\displaystyle e^x - 1 = 0</math>. Błąd kolejnych przybliżeń wyświetlany jest w skali
logarytmicznej, dzięki czemu lepiej widać różnicę między zbieżnością liniową a
kwadratową.]]


Powiemy, że metoda iteracyjna <math>\displaystyle \phi</math> jest w klasie funkcji <math>\displaystyle F</math>  
Powiemy, że metoda iteracyjna <math>\displaystyle \phi</math> jest w klasie funkcji <math>\displaystyle F</math>  
Linia 461: Linia 476:
}}
}}


[[Image:MNXXX.png|thumb|300px|Zbieżność metody Newtona na tle metody bisekcji]]
[[Image:MNmultiplezeros.png|frame|400px|center|Zbieżność metody Newtona dla zer wielokrotnych <math>\displaystyle f(x)
 
= (x-1)^5</math> jest liniowa z ilorazem <math>\displaystyle \frac{4}{5}</math> (końcowe załamanie wykresu
[[Image:MNXXX.png|thumb|300px|Zbieżność metody Newtona dla zer wielokrotnych]]
spowodowane jest przypadkowym trafieniem w dokładne miejsce zerowe). Metoda bisekcji nie jest na to czuła i dalej zbiega z ilorazem
<math>\displaystyle \frac{1}{2}</math>.]]


===Metoda siecznych===
==Metoda siecznych==


Inną znaną i często używaną metodą iteracyjną, opartą na podobnym pomyśle
Inną znaną i często używaną metodą iteracyjną, opartą na podobnym pomyśle
linearyzacyjnych co metoda Newtona  
linearyzacyjnych co metoda Newtona  
jest ''metoda siecznych'', w której zamiast przybliżenia wykresu <math>\displaystyle f</math> przez
jest <strong>metoda siecznych</strong>, w której zamiast przybliżenia wykresu <math>\displaystyle f</math> przez
styczną,  stosuje się  przybliżenie sieczną.
styczną,  stosuje się  przybliżenie sieczną.
   
   
Linia 513: Linia 529:
wypadku okazuje się, że metoda stycznych --- choć wolniej zbieżna niż metoda
wypadku okazuje się, że metoda stycznych --- choć wolniej zbieżna niż metoda
stycznych --- dzięki temu, że
stycznych --- dzięki temu, że
jej iteracja wymaga jedynie wyznaczenia jednej wartości <math>\displaystyle f</math>, jest ''bardziej
jej iteracja wymaga jedynie wyznaczenia jednej wartości <math>\displaystyle f</math>, jest <strong>bardziej
efektywna'' od metody Newtona: koszt osiągnięcia zadanej dokładności jest w
efektywna</strong> od metody Newtona: koszt osiągnięcia zadanej dokładności jest w
takim przypadku mniejszy od analogicznego kosztu dla metody Newtona.
takim przypadku mniejszy od analogicznego kosztu dla metody Newtona.


Linia 521: Linia 537:
przy odejmowaniu.
przy odejmowaniu.


[[Image:MNXXX.png|thumb|300px|Zbieżność metody siecznych na tle metody Newtona]]
[[Image:MNstycznesiecznebisekcja.png|frame|400px|center|Porównanie zbieżności metody bisekcji,
stycznych i siecznych  
dla równania <math>\displaystyle e^x - 1 = 0</math>. Błąd kolejnych przybliżeń wyświetlany jest w skali
logarytmicznej.]]


===Metoda Brenta===
==Metoda Brenta==


Naturalnie, uważny student zaczyna zadawać sobie pytanie, czy nie można w jakiś
Naturalnie, uważny student zaczyna zadawać sobie pytanie, czy nie można w jakiś
Linia 532: Linia 551:
studentom gorszych uczelni).  
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
Okazuje się, że można to zrobić, wprowadzając metodę opartą na <strong>trzech</strong> punktach lokalizujących miejsce zerowe: dwóch odcinających zero tak jak
w metodzie bisekcji i trzecim, konstruowanym jak np. w metodzie stycznych. W
w metodzie bisekcji i trzecim, konstruowanym jak np. w metodzie stycznych. W
kolejnej iteracji konstruujemy wymieniamy jeden z punktów albo wedle metody
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ę
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
(aby zagwarantować sobie, że w wiadomym przedziale miejsce zerowe rzeczywiście
sie znajduje).
się znajduje).


Ten prosty pomysł metody hybrydowej wymaga jednak subtelnego dopracowania, co
Ten prosty pomysł metody hybrydowej wymaga jednak subtelnego dopracowania, co
Linia 543: Linia 562:
Dekkera, van Wijngaardena i Dijkstry.
Dekkera, van Wijngaardena i Dijkstry.


Funkcja MATLABa (i Octave'a) <code>fzero</code> implementuje metodę Brenta.
[[grafika:Brent.jpg|thumb|right||Richard Brent<br>  [[Biografia Brent|Zobacz biografi�]]]]
 
Funkcja MATLABa (i Octave'a) <code>fzero</code> implementuje właśnie metodę Brenta.
Ciekawostką jest, że autorem implementacji w Octave jest ówczesny student
Ciekawostką jest, że autorem implementacji w Octave jest ówczesny student
matematyki na Uniwersytecie Warszawskim, Łukasz Bodzon.
matematyki na Uniwersytecie Warszawskim, Łukasz Bodzon.


[[Image:MNXXX.png|thumb|300px|Porównanie zbieżności różnych metod dla równań nieliniowych: zero
==Metody dla układów równań nieliniowych==
jednokrotne]]
 
[[Image:MNXXX.png|thumb|300px|Porównanie zbieżności różnych metod dla równań nieliniowych: zero
wielokrotne]]

Wersja z 17:30, 1 wrz 2006

Równania nieliniowe

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)xϵsin(x)=0

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

  • znajdowanie miejsc zerowych wielomianu:
f(x)anxn++a1x+a0=0

Bardzo wiele modeli matematycznych wymaga rozwiązania równania z wielomianową nieliniowością. Dodaj link: 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 Dodaj link: znajdowania wartości własnych macierzy.

  • znajdowanie miejsc zerowych trójmianu kwadratowego:
f(x)a2x2+a1x+a0=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)x2a=0

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

Plik:Heron.jpg
Heron
Zobacz biografi�
xk+1=12(xk+axk)

daje bardzo dobre przybliżenie a już po kilku iteracjach.

  • implementacja wyznaczania odwrotności liczby a (bez dzielenia!):
f(x)1xa=0

Wciąż spotykane zadanie, np. tak można w praktyce poprawić precyzję [http://www.amd.com/us-en/assets/content_type/white_papers_and_tech_docs/21928.pdf funkcji wektorowych stosowanych w niektórych procesorach AMD]. 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.

Bisekcja

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={fC([a,b]):f(a)f(b)<0}.

Oczywiście, każda funkcja fF 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.

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

|xx*|(12)k(ba2).

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 |ba| 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ą powyższego oszacowania błędu jest bowiem następujący wniosek.

Wniosek

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

k=k(ϵ)=log2(ba)ϵ1

wartości funkcji.

Iteracja prosta 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ę ϕ) do równania równoważnego (tzn. mającego te same rozwiązania)

x=ϕ(x).
Następnie, startując z pewnego przybliżenia 
początkowego x0, konstruujemy ciąg kolejnych 

przybliżeń xk według wzoru

xk=ϕ(xk1),k1.
Stefan Banach
Zobacz biografi�

Twierdzenie Banacha, o zbieżności iteracji prostej

Niech D0 będzie domkniętym podzbiorem dziedziny D,

D0=D0D,

w którym ϕ jest odwzorowaniem zwężającym. To znaczy, ϕ(D0)D0, oraz istnieje stała 0L<1 taka, że

ϕ(x)ϕ(y)Lxy,x,yD0.

Wtedy równanie

x=ϕ(x).

ma dokładnie jedno rozwiązanie x*, oraz

x*=limkxk,

dla dowolnego przybliżenia początkowego x0D0.

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 ks 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 {xk}k jest więc ciągiem Cauchy'ego. Stąd istnieje granica α=limkxk, która należy do D0, wobec domkniętości tego zbioru. Ponieważ lipschitzowskość ϕ implikuje jej ciągłość, mamy też

ϕ(α)=ϕ(limkxk)=limkϕ(xk)=limkxk=α,

tzn. α jest punktem stałym odwzorowania ϕ. Dla jednoznaczności zauważmy, że jeśliby istniał drugi, różny od α, punkt stały β, to mielibyśmy

αβ=ϕ(α)ϕ(β)Lαβ.

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

ϕ jest zwężająca.

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

Wniosek

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

xkx*Lkx0x*.

Przykład

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

x=cos(x),dlaxD=R.

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

L=max0x1|cos(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 x0[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 ϕ jest zwężająca na całym zbiorze D, tzn. D0=D. Jeśli ponadto D ma skończoną średnicę diam(D), to dla osiągnięcia ϵ-aproksymacji zera funkcji f wystarczy wykonać

k=k(ϵ)=log(x0x*/ϵ)log(1/L)=log(diam(D)/ϵ)log(1/L)

iteracji, niezależnie od x0. 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 12, dopiero po piątej 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 {xk}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!

Plik:Newton.jpg
Isaac Newton
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. Zobacz biografi�

Startując z pewnego przybliżenia początkowego x0, w kolejnych krokach metody, k-te przybliżenie xk jest punktem przecięcia stycznej do wykresu f w punkcie xk1. Ponieważ równanie stycznej wynosi y(x)=f(xk1)+f(xk1)(xxk1), otrzymujemy wzór

xk=xk1f(xk1)f(xk1).

Oczywiście, aby metoda Newtona była dobrze zdefiniowana, musimy założyć, że f(xk1) istnieje i nie jest zerem.

<flash>file=Newton.swf</flash><div.thumbcaption>Postęp iteracji Newtona

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

ϕ(x)=xf(x)f(x).

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

Nawet jeśli pochodna w xk1 się nie zeruje, ciąg {xk}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(x*)0.

Ponadto załóżmy, że f jest dwukrotnie różniczkowalna w sposób ciągły, fC2(D). Rozwijając ϕ w szereg Taylora w punkcie x* otrzymujemy

xkx*=ϕ(xk1)ϕ(x*)=(xk1x*)ϕ(x*)+(xk1x*)2ϕ(ξk)/2,

gdzie min(x*,xk1)ξkmax(x*,xk1). Wobec tego, że ϕ(x*)=f(x)f(x)/(f(x))2=0 i ϕ(ξk)=f(ξk)/f(ξk), mamy

xkx*=(xk1x*)2f(ξk)2f(ξk).

Zdefiniujmy liczbę

Rf=supr0sup{x:|xx*|r}|2(xx*)f(x)f(x)|<1.

Oczywiście Rf>0. Dla xk1 spełniającego |xk1x*|R<Rf, mamy z poprzedniej równości

|xkx*|q|xk1x*|,

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

Niech teraz x* będzie zerem m-krotnym,

f(x*)=f(x*)==f(m1)(x*)=0f(m)(x*),

gdzie m2, 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 xk1 jest "blisko" x*.

Metoda Newtona jest więc zbieżna lokalnie. Z powyższego można też wywnioskować, jaki jest charakter zbieżności metody Newtona. Dla zera jednokrotnego x* oraz f(x*)0 mamy bowiem

|xkx*||xxk1|2|f(x*)|2|f(x*)|.

Mówimy, że zbieżność metody Newtona, gdy f(x*)0 jest kwadratowa.

Stwierdzenie

Jeśli f(x*)0 oraz f(x*)=0 to zbieżność jest nawet szybsza. Z kolei dla zera m-krotnego (tzn. f(x*)=f(x*)=f(m)(x*)=0, m>1) zbieżność jest liniowa z ilorazem (11m).

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.

Porównanie zbieżności metody bisekcji i stycznych dla równania ex1=0. Błąd kolejnych przybliżeń wyświetlany jest w skali logarytmicznej, dzięki czemu lepiej widać różnicę między zbieżnością liniową a kwadratową.

Powiemy, że metoda iteracyjna ϕ jest w klasie funkcji F rzędu co najmniej p1, gdy spełniony jest następujący warunek. Niech fF i f(x*)=0. Wtedy istnieje stała C< taka, że dla dowolnych przybliżeń początkowych x0,,xs1 dostatecznie bliskich x*, kolejne przybliżenia xk=ϕ(xk1,,xks) generowane tą metodą spełniają

|xkx*|C|xk1x*|p.

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

Definicja

Wykładnikiem zbieżności metody iteracyjnej ϕ w klasie F nazywamy liczbę p* zdefiniowaną równością

p*=sup{p1:ϕ 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 dla zer wielokrotnych f(x)=(x1)5 jest liniowa z ilorazem 45 (końcowe załamanie wykresu spowodowane jest przypadkowym trafieniem w dokładne miejsce zerowe). Metoda bisekcji nie jest na to czuła i dalej zbiega z ilorazem 12.

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 xk przybliżenia xk1 i xk2. Musimy również wybrać dwa różne punkty startowe x0 i x1. Ponieważ prosta interpolująca f w xk1 i xk2 ma wzór

y(x)=xxk2xk1xk2f(xk1)+xxk1xk2xk1f(xk2),

otrzymujemy

xk=xk1xk1xk2f(xk1)f(xk2)f(xk1).

Zauważmy, że jeśli xk1 i xk2 są blisko siebie, to xk jest podobny do tego z metody Newtona, bowiem wtedy iloraz różnicowy

f(xk1)f(xk2)xk1xk2f(xk1).

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*=1+52=1.618. 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(xk), jest tak samo, albo i bardziej kosztowne od wyznaczenia wartości f(xk). 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.

Porównanie zbieżności metody bisekcji, stycznych i siecznych dla równania ex1=0. Błąd kolejnych przybliżeń wyświetlany jest w skali logarytmicznej.

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 się 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.

Plik:Brent.jpg
Richard Brent
Zobacz biografi�

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

Metody dla układów równań nieliniowych