|
|
| Linia 10: |
Linia 10: |
|
| |
|
| [[Image:Fraktal.png|400px|To jest tytuł]] | | [[Image:Fraktal.png|400px|To jest tytuł]] |
|
| |
| W wielu zadaniach, m.in. matematyki stosowanej, spotykamy się z problemem
| |
| rozwiązania skalarnego równania nieliniowego postaci <math>\displaystyle f(x) = 0</math>:
| |
| * rozwiązywanie równania Keplera
| |
|
| |
| <center><math>\displaystyle f(x) \equiv x - \epsilon \sin(x) = 0</math></center>
| |
|
| |
| To równanie jest bardzo ważne w astronomii.
| |
| * znajdowanie miejsc zerowych wielomianu:
| |
|
| |
| <center><math>\displaystyle f(x) \equiv a_nx^n + \ldots +
| |
| a_1x + a_0 = 0</math></center>
| |
|
| |
| 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 [[sec:eigenvalue|Uzupełnij: znajdowania wartości własnych
| |
| macierzy]].
| |
| * znajdowanie miejsc zerowych trójmianu kwadratowego:
| |
|
| |
| <center><math>\displaystyle f(x) \equiv a_2x^2 +
| |
| a_1x + a_0 = 0.
| |
| </math></center>
| |
|
| |
| 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 <math>\displaystyle a</math>:
| |
|
| |
| <center><math>\displaystyle f(x) \equiv
| |
| x^2 - a = 0</math></center>
| |
|
| |
| Jeden ze sposobów na implementację funkcji "<code>sqrt()</code>". Szybkie algorytmy
| |
| wyznaczania pierwiastka kwadratowego były znane już starożytnym. W wykładzie
| |
| zrozumiemy, dlaczego metoda Herona,
| |
| <center><math>\displaystyle
| |
| x_{k+1} = \frac{1}{2}\left(x_k + \frac{a}{x_k}\right)
| |
| </math></center>
| |
|
| |
| 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!):
| |
|
| |
| <center><math>\displaystyle f(x) \equiv
| |
| \frac{1}{x} - a = 0</math></center>
| |
|
| |
| 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 <math>\displaystyle f(x) = 0</math>.
| |
|
| |
| ''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 <math>\displaystyle [a,b]</math>
| |
| i zmieniających znak. Dokładniej, rozpatrzmy klasę
| |
| funkcji
| |
|
| |
| <center><math>\displaystyle
| |
| F\,=\,\{\,f\in C([a,b])\,:\;f(a)\cdot f(b) < 0\,\}.
| |
| </math></center>
| |
|
| |
| Oczywiście, każda funkcja <math>\displaystyle f\in F</math> ma co najmniej jedno
| |
| zero w <math>\displaystyle [a,b]</math>. Startując z przedziału <math>\displaystyle [a,b]</math>, w
| |
| kolejnych krokach metody bisekcji obliczamy informację
| |
| o wartości <math>\displaystyle f</math> 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.
| |
|
| |
| <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
| |
| poleceń, po wykonaniu którego <math>\displaystyle x</math> jest przybliżeniem
| |
| zera funkcji <math>\displaystyle f</math> z zadaną dokładnością <math>\displaystyle \epsilon</math>.
| |
|
| |
| {{algorytm|{Metoda bisekcji}||
| |
| <pre>
| |
|
| |
| xl <nowiki>=</nowiki> a; xr <nowiki>=</nowiki> b;
| |
| x <nowiki>=</nowiki> (a+b)/2; e <nowiki>=</nowiki> (b-a)/2;
| |
| while (e > <math>\displaystyle \epsilon</math>)
| |
| {
| |
| if (f(x)*f(xl) < 0)
| |
| xr <nowiki>=</nowiki> x;
| |
| else
| |
| xl <nowiki>=</nowiki> x;
| |
| x <nowiki>=</nowiki> (xl+xr)/2; e <nowiki>=</nowiki> e/2;
| |
| }
| |
| </pre>}}
| |
|
| |
| Z konstrukcji metody łatwo wynika, że po wykonaniu
| |
| <math>\displaystyle k</math> iteracji (czyli po obliczeniu <math>\displaystyle k</math> wartości funkcji)
| |
| otrzymujemy <math>\displaystyle x</math>, które odległe jest od pewnego
| |
| rozwiązania <math>\displaystyle x^*</math> o co najwyżej
| |
|
| |
| <center><math>\displaystyle
| |
| |x-x^*|\,\le\,\Big(\frac 12\Big)^k
| |
| \Big(\frac{b-a}{2}\Big).
| |
| </math></center>
| |
|
| |
| Metoda bisekcji jest więc zbieżna ''liniowo'' z
| |
| ilorazem <math>\displaystyle 1/2</math>. 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 <math>\displaystyle a</math> i <math>\displaystyle b</math>
| |
| 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ść
| |
| 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
| |
| wystarcza jedynie ''ciągłość'' funkcji. Poza tym
| |
| możemy łatwo kontrolować ''błąd bezwzględny aproksymacji miejsca zerowego''. Konsekwencją
| |
| ([[##blbis|Uzupelnic: blbis ]]) jest bowiem następujący wniosek.
| |
|
| |
| {{wniosek|||
| |
| Dla znalezienia zera <math>\displaystyle x^*</math> z dokładnością
| |
| <math>\displaystyle \epsilon>0</math>, wystarczy obliczyć w metodzie bisekcji
| |
|
| |
| <center><math>\displaystyle k\,=\,k(\epsilon)\,=\,
| |
| \Big\lceil{\log_2\frac{(b-a)}{\epsilon}}\Big\rceil - 1
| |
| </math></center>
| |
|
| |
| wartości funkcji.
| |
| }}
| |
|
| |
| [[Image:MN|400px|Zbieżność metody bisekcji dla ....]]
| |
|
| |
| Dodajmy jeszcze, że bisekcja minimalizuje błąd
| |
| najgorszy w klasie <math>\displaystyle F</math> zdefiniowanej [[dfkl|Uzupełnij: na początku tej sekcji]],
| |
| wśród wszystkich algorytmów korzystających
| |
| z określonej liczby obliczeń wartości funkcji,
| |
| zob. [[dwopt|Uzupełnij: uwaga na końcu wykładu]].
| |
|
| |
| {{uwaga|||
| |
|
| |
| Metoda bisekcji jest
| |
| optymalna w następującym sensie. Niech <math>\displaystyle A:F\to R</math>
| |
| będzie dowolną metodą (algorytmem) aproksymującą
| |
| zero <math>\displaystyle x^*(f)</math> funkcji <math>\displaystyle f</math> z klasy <math>\displaystyle F</math> zdefiniowanej
| |
| w ([[##dfkl|Uzupelnic: dfkl ]]), korzystającą jedynie z obliczeń
| |
| (informacji o) <math>\displaystyle f</math> w <math>\displaystyle k</math> punktach. Wtedy dla dowolnego
| |
| <math>\displaystyle \gamma>0</math> istnieje funkcja <math>\displaystyle f_\gamma\in F</math> mająca
| |
| tylko jedno zero <math>\displaystyle x^*(f_\gamma)</math> w <math>\displaystyle [a,b]</math> i taka, że
| |
|
| |
| <center><math>\displaystyle |A(f_\gamma)-x^*(f_\gamma)|\,\ge\,
| |
| \Big(\frac 12\Big)^k\Big(\frac{b-a}{2}\Big)+\gamma.
| |
| </math></center>
| |
|
| |
| Co więcej, można pokazać, że fakt ten zachodzi
| |
| też w węższej klasie <math>\displaystyle F_1</math> funkcji <math>\displaystyle f\in F</math>, które
| |
| są dowolnie wiele razy różniczkowalne. Oczywiście,
| |
| nie wyklucza to istnienia metod iteracyjnych (takich
| |
| jak metoda Newtona), które dla <math>\displaystyle f\in F_1</math> 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
| |
|
| |
| <center><math>\displaystyle
| |
| f(x) = 0
| |
| </math></center>
| |
|
| |
| przekształcamy (dobierając odpowiednią funkcję <math>\displaystyle \phi</math>) do równania równoważnego
| |
| (tzn. mającego te same rozwiązania)
| |
|
| |
| <center><math>\displaystyle
| |
| x\,=\,\phi( x).
| |
| </math></center>
| |
|
| |
| Następnie, startując z pewnego przybliżenia
| |
| początkowego <math>\displaystyle x_0</math>, konstruujemy ciąg kolejnych
| |
| przybliżeń <math>\displaystyle x_k</math> według wzoru
| |
|
| |
| <center><math>\displaystyle x_k\,=\,\phi( x_{k-1}),\qquad k\ge 1.
| |
| </math></center>
| |
|
| |
| {{twierdzenie|Banacha, o zbieżności iteracji prostej||
| |
| Niech <math>\displaystyle D_0</math> będzie domkniętym
| |
| podzbiorem dziedziny <math>\displaystyle D</math>,
| |
|
| |
| <center><math>\displaystyle \overline D_0\,=\,D_0\subset D,
| |
| </math></center>
| |
|
| |
| w którym <math>\displaystyle \phi</math> jest odwzorowaniem zwężającym.
| |
| To znaczy, <math>\displaystyle \phi(D_0)\subset D_0</math>, oraz istnieje stała
| |
| <math>\displaystyle 0\le L<1</math> taka, że
| |
|
| |
| <center><math>\displaystyle \|\phi( x)-\phi( y)\|\,\le\,L\,\| x- y\|,
| |
| \qquad\forall x, y\in D_0.
| |
| </math></center>
| |
|
| |
| Wtedy równanie ([[##rrw|Uzupelnic: rrw ]]) ma dokładnie jedno
| |
| rozwiązanie <math>\displaystyle x^*</math>, oraz
| |
|
| |
| <center><math>\displaystyle x^*\,=\,\lim_{k\to\infty} x_k,
| |
| </math></center>
| |
|
| |
| dla dowolnego przybliżenia początkowego
| |
| <math>\displaystyle x_0\in D_0</math>.
| |
| }}
| |
|
| |
| {{dowod|||
| |
| Wobec
| |
|
| |
| <center><math>\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</math></center>
| |
|
| |
| dla <math>\displaystyle k\ge s</math> mamy
| |
|
| |
| <center><math>\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</math></center>
| |
|
| |
| Ciąg <math>\displaystyle \{ x_k\}_k</math> jest więc ciągiem Cauchy'ego.
| |
| Stąd istnieje granica
| |
| <math>\displaystyle \vec\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ż
| |
| lipschitzowskość <math>\displaystyle \phi</math> implikuje jej ciągłość,
| |
| mamy też
| |
|
| |
| <center><math>\displaystyle \phi(\vec\alpha)\,=\,\phi\Big(\lim_{k\to\infty} x_k\Big)
| |
| \,=\,\lim_{k\to\infty}\phi( x_k)
| |
| \,=\,\lim_{k\to\infty} x_k\,=\,\vec\alpha,
| |
| </math></center>
| |
|
| |
| tzn. <math>\displaystyle \vec\alpha</math> jest punktem stałym odwzorowania <math>\displaystyle \phi</math>.
| |
| 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>,
| |
| to mielibyśmy
| |
|
| |
| <center><math>\displaystyle \|\vec\alpha-\vec\beta\|\,=\,
| |
| \|\phi(\vec\alpha)-\phi(\vec\beta)\|
| |
| \,\le\,L\,\|\vec\alpha-\vec\beta\|.
| |
| </math></center>
| |
|
| |
| Stąd <math>\displaystyle 1<L</math>, co jest sprzeczne z założeniem, że
| |
| <math>\displaystyle \phi</math> jest zwężająca. }}
| |
|
| |
| Z powyższych rozważań otrzymujemy natychmiastowy
| |
| wniosek dotyczący zbieżności iteracji prostych.
| |
|
| |
| {{wniosek|||
| |
| Przy założeniach [[twit|Uzupełnij: twierdzenia Banacha]],
| |
| metoda iteracji prostych jest zbieżna co
| |
| najmniej liniowo z ilorazem <math>\displaystyle L</math>, tzn.
| |
|
| |
| <center><math>\displaystyle \| x_k- x^*\|\,\le\,L^k\,\| x_0- x^*\|.
| |
| </math></center>
| |
|
| |
| }}
| |
|
| |
| {{przyklad|||
| |
| Dla ilustracji, rozpatrzmy natępujące proste
| |
| równanie skalarne:
| |
|
| |
| <center><math>\displaystyle
| |
| x\,=\,\cos(x), \qquad \mbox{dla} \qquad x\in D= R.
| |
| </math></center>
| |
|
| |
| W tym przypadku <math>\displaystyle \phi(x)=\cos(x)</math>. Zauważamy, że w
| |
| przedziale <math>\displaystyle [0,1]</math> funkcja <math>\displaystyle \phi</math> jest zwężająca ze
| |
| stałą
| |
|
| |
| <center><math>\displaystyle L\,=\,\max_{0\le x\le 1}|\cos'(x)|\,=\,\sin(1)\,<\,1.
| |
| </math></center>
| |
|
| |
| Stąd istnieje dokładnie jedno rozwiązanie naszego równania
| |
| w przedziale <math>\displaystyle [0,1]</math>. 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 <math>\displaystyle x_0\in [0,1]</math>.
| |
| }}
| |
|
| |
| Zaletą iteracji prostych jest fakt, że zbieżność
| |
| nie zależy od wymiaru <math>\displaystyle n</math> zadania, ale tylko od stałej
| |
| Lipschitza <math>\displaystyle L</math> (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 <math>\displaystyle \phi</math> jest zwężająca na całym
| |
| zbiorze <math>\displaystyle D</math>, tzn. <math>\displaystyle D_0=D</math>. Jeśli ponadto <math>\displaystyle D</math> ma
| |
| skończoną średnicę <math>\displaystyle \mbox{diam} (D)</math>, to dla
| |
| osiągnięcia <math>\displaystyle \epsilon</math>-aproksymacji zera funkcji <math>\displaystyle f</math>
| |
| wystarczy wykonać
| |
|
| |
| <center><math>\displaystyle k\,=\,k(\epsilon)\,=\,\Big\lceil\frac
| |
| {\log(\| x_0- x^*\|/\epsilon)}{\log(1/L)}\Big\rceil
| |
| \,=\,\Big\lceil\frac
| |
| {\log( \mbox{diam} (D)/\epsilon)}{\log(1/L)}\Big\rceil
| |
| </math></center>
| |
|
| |
| iteracji, niezależnie od <math>\displaystyle x_0</math>. 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
| |
| <math>\displaystyle \frac{1}{2}</math>, dopiero co 5 iteracji, dostajemy kolejną
| |
| 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.
| |
| 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 <math>\displaystyle \{x_k\}_k</math> do zera danej funkcji <math>\displaystyle f</math>
| |
| jest zapewniona jedynie wtedy, gdy przybliżenia początkowe
| |
| zostały wybrane dostatecznie blisko <math>\displaystyle x^*</math>.
| |
|
| |
| W dalszych rozważaniach będziemy zakładać dla
| |
| 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ą
| |
| linią prostą, a dla niej już umiemy znaleźć miejsce zerowe!
| |
|
| |
| Startując z pewnego przybliżenia
| |
| początkowego <math>\displaystyle x_0</math>, w kolejnych krokach metody, <math>\displaystyle k</math>-te
| |
| przybliżenie <math>\displaystyle x_k</math> jest punktem przecięcia stycznej do
| |
| wykresu <math>\displaystyle f</math> w punkcie <math>\displaystyle x_{k-1}</math>. Ponieważ równanie
| |
| stycznej wynosi <math>\displaystyle y(x)=f(x_{k-1})+f'(x_{k-1})(x-x_{k-1})</math>,
| |
| otrzymujemy wzór
| |
|
| |
| <center><math>\displaystyle x_k\,=\,x_{k-1}\,-\,\frac{f(x_{k-1})}{f'(x_{k-1})}.
| |
| </math></center>
| |
|
| |
| Oczywiście, aby metoda Newtona była dobrze zdefiniowana,
| |
| musimy założyć, że <math>\displaystyle f'(x_{k-1})</math> istnieje i nie
| |
| jest zerem.
| |
|
| |
| Zauważmy, że metodę Newtona można traktować jako
| |
| szczególny przypadek iteracji prostych, gdzie
| |
|
| |
| <center><math>\displaystyle \phi(x)\,=\,x-\frac{f(x)}{f'(x)}.
| |
| </math></center>
| |
|
| |
| Widać też, że nie jest ona zbieżna globalnie.
| |
|
| |
| Nawet jeśli pochodna w <math>\displaystyle x_{k-1}</math> się nie zeruje,
| |
| ciąg <math>\displaystyle \{x_k\}_k</math> może nie zbiegać do zera funkcji <math>\displaystyle f</math>.
| |
| Okazuje się jednak, że jeśli
| |
| wystartujemy dostatecznie blisko rozwiązania <math>\displaystyle x^*</math>, to
| |
| metoda Newtona jest zbieżna. Dokładniej, załóżmy
| |
| najpierw, że <math>\displaystyle f(x^*)=0</math> oraz
| |
|
| |
| <center><math>\displaystyle f'(x^*)\,\ne\,0.
| |
| </math></center>
| |
|
| |
| Ponadto załóżmy, że <math>\displaystyle f</math> jest dwukrotnie
| |
| różniczkowalna w sposób ciągły, <math>\displaystyle f\in C^2(D)</math>.
| |
| Rozwijając <math>\displaystyle \phi</math> w szereg Taylora w punkcie <math>\displaystyle x^*</math>
| |
| otrzymujemy
| |
|
| |
| <center><math>\displaystyle x_k-x^*\,=\,\phi(x_{k-1})-\phi(x^*)\,=\,
| |
| (x_{k-1}-x^*)\phi'(x^*)+(x_{k-1}-x^*)^2\phi''(\xi_k)/2,
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle \min(x^*,x_{k-1})\le\xi_k\le\max(x^*,x_{k-1})</math>.
| |
| Wobec tego, że <math>\displaystyle \phi'(x^*)=f(x)f''(x)/(f'(x))^2=0</math> i
| |
| <math>\displaystyle \phi''(\xi_k)=f''(\xi_k)/f'(\xi_k)</math>, mamy
| |
|
| |
| <center><math>\displaystyle
| |
| x_k-x^*\,=\,(x_{k-1}-x^*)^2\frac{f''(\xi_k)}{2f'(\xi_k)}.
| |
| </math></center>
| |
|
| |
| Zdefiniujmy liczbę
| |
|
| |
| <center><math>\displaystyle R_f\,=\,\sup_{r\ge 0}\sup_{\{x:|x-x^*|\le r\}}
| |
| \Big|\frac{2(x-x^*)f''(x)}{f'(x)}\Big|\,<\,1.
| |
| </math></center>
| |
|
| |
| Oczywiście <math>\displaystyle R_f>0</math>. Dla <math>\displaystyle x_{k-1}</math> spełniającego
| |
| <math>\displaystyle |x_{k-1}-x^*|\le R<R_f</math>, mamy z poprzedniej równości ([[##nrpdst|Uzupelnic: nrpdst ]])
| |
|
| |
| <center><math>\displaystyle |x_k-x^*|\,\le\,q\,|x_{k-1}-x^*|,
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle q<1</math> i <math>\displaystyle q</math> zależy tylko od <math>\displaystyle R</math>.
| |
|
| |
| Niech teraz <math>\displaystyle x^*</math> będzie zerem <math>\displaystyle m</math>-krotnym,
| |
|
| |
| <center><math>\displaystyle f(x^*)=f'(x^*)=\cdots =f^{(m-1)}(x^*)=0\ne f^{(m)}(x^*),
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle m\ge 2</math>, oraz niech <math>\displaystyle f</math> będzie <math>\displaystyle m</math>-krotnie
| |
| różniczkowalna w sposób ciągły. Wtedy
| |
|
| |
| <center><math>\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</math></center>
| |
|
| |
| o ile <math>\displaystyle x_{k-1}</math> jest "blisko" <math>\displaystyle x^*</math>.
| |
|
| |
| Metoda Newtona jest więc zbieżna lokalnie.
| |
| Z ([[##nrpdst|Uzupelnic: nrpdst ]]) i ([[##nrtp|Uzupelnic: nrtp ]]) można też wywnioskować,
| |
| jaki jest charakter zbieżności metody Newtona. Dla zera
| |
| 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^*)}.
| |
| </math></center>
| |
|
| |
| Mówimy, że zbieżność jest ''kwadratowa''. Jeśli zaś
| |
| <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
| |
| <math>\displaystyle (1-\frac{1}{m})</math>.
| |
|
| |
| 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 <math>\displaystyle \phi</math> jest w klasie funkcji <math>\displaystyle F</math>
| |
| rzędu co najmniej <math>\displaystyle p\ge 1</math>, gdy spełniony jest następujący
| |
| warunek. Niech <math>\displaystyle f\in F</math> i <math>\displaystyle f(x^*)=0</math>. Wtedy istnieje stała
| |
| <math>\displaystyle C<\infty</math> taka, że dla dowolnych przybliżeń początkowych
| |
| <math>\displaystyle x_0,\ldots,x_{s-1}</math> dostatecznie bliskich <math>\displaystyle x^*</math>, kolejne
| |
| przybliżenia <math>\displaystyle x_k=\phi(x_{k-1},\ldots,x_{k-s})</math> generowane
| |
| tą metodą spełniają
| |
|
| |
| <center><math>\displaystyle |x_k-x^*|\,\le\,C\,|x_{k-1}-x^*|^p.
| |
| </math></center>
| |
|
| |
| Ponadto, jeśli <math>\displaystyle p=1</math> to dodatkowo żąda się, aby <math>\displaystyle C<1</math>.
| |
|
| |
| {{definicja|||
| |
| Wykładnikiem zbieżności metody
| |
| iteracyjnej <math>\displaystyle \phi</math> w klasie <math>\displaystyle F</math> nazywamy liczbę <math>\displaystyle p^*</math>
| |
| zdefiniowaną równością
| |
|
| |
| <center><math>\displaystyle p^*\,=\,\sup\,\{\,p\ge 1:\,\phi
| |
| \mbox{ jest rzędu co najmniej } p\,\}.
| |
| </math></center>
| |
|
| |
| }}
| |
|
| |
| 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 <math>\displaystyle p^*=2</math> w klasie funkcji o zerach
| |
| jednokrotnych, oraz <math>\displaystyle p^*=1</math> w klasie funkcji o zerach
| |
| wielokrotnych.
| |
| }}
| |
|
| |
| [[Image:MN|400px|Zbieżność metody Newtona na tle metody bisekcji]]
| |
|
| |
| [[Image:MN|400px|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 <math>\displaystyle f</math> przez
| |
| styczną, stosuje się
| |
| przybliżenie sieczną.
| |
|
| |
| Metoda ta
| |
| wykorzystuje więc do konstrukcji <math>\displaystyle x_k</math> przybliżenia
| |
| <math>\displaystyle x_{k-1}</math> i <math>\displaystyle x_{k-2}</math>. Musimy również wybrać dwa różne
| |
| punkty startowe <math>\displaystyle x_0</math> i <math>\displaystyle x_1</math>. Ponieważ prosta interpolująca
| |
| <math>\displaystyle f</math> w <math>\displaystyle x_{k-1}</math> i <math>\displaystyle x_{k-2}</math> ma wzór
| |
|
| |
| <center><math>\displaystyle 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}),
| |
| </math></center>
| |
|
| |
| otrzymujemy
| |
|
| |
| <center><math>\displaystyle x_k\,=\,x_{k-1}\,-\,\frac{x_{k-1}-x_{k-2}}
| |
| {f(x_{k-1})-f(x_{k-2})}\,f(x_{k-1}).
| |
| </math></center>
| |
|
| |
| Zauważmy, że jeśli <math>\displaystyle x_{k-1}</math> i <math>\displaystyle x_{k-2}</math> są blisko
| |
| siebie, to <math>\displaystyle x_k</math> jest podobny do tego z metody Newtona,
| |
| bowiem wtedy iloraz różnicowy
| |
| <center><math>\displaystyle
| |
| \frac{f(x_{k-1})-f(x_{k-2})}{x_{k-1}-x_{k-2}} \approx f'(x_{k-1}).
| |
| </math></center>
| |
|
| |
| Nie wystarcza to
| |
| jednak, aby osiągnąć zbieżność z wykładnikiem
| |
| <math>\displaystyle 2</math>. Dokładniej, można pokazać, że wykładnik
| |
| zbieżności metody siecznych dla zer jednokrotnych dostatecznie gładkich funkcji
| |
| wynosi <math>\displaystyle p^*=\frac{1+\sqrt{5}}{2}=1.618\ldots</math>. 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, <math>\displaystyle f'(x_k)</math>, jest
| |
| tak samo, albo i bardziej kosztowne od wyznaczenia wartości <math>\displaystyle f(x_k)</math>. 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 <math>\displaystyle f</math>, 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.
| |
|
| |
| [[Image:MN|400px|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) <code>fzero</code> implementuje metodę Brenta.
| |
| Ciekawostką jest, że autorem implementacji w Octave jest ówczesny student
| |
| matematyki na Uniwersytecie Warszawskim, Łukasz Bodzon.
| |
|
| |
| [[Image:MN|400px|Porównanie zbieżności różnych metod dla równań nieliniowych: zero
| |
| jednokrotne]]
| |
|
| |
| [[Image:MN|400px|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.
| |
