MN02

Równania nieliniowe

<<< Powrót do strony głównej przedmiotu Metody numeryczne

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

Przykład

Równanie Keplera

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

jest bardzo ważne w astronomii, jego rozwiązanie pozwala wyznaczyć przyszłe położenie planety. Parametr $ϵ$ odpowiada ekscentryczności orbity i przyjmuje wartości z przedziału $[0,1]$. Poza paru prostymi przypadkami, w ogólności równanie Keplera nie daje się rozwiązać w terminach funkcji elementarnych.

Johann Kepler
Zobacz biografię

Przykład

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 wyspecjalizowanych 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 znajdowania wartości własnych macierzy.

Przykład

Obliczanie pierwiastka kwadratowego z zadanej liczby $a$, czyli sposób na implementację funkcji "sqrt()". Można to zadanie wyrazić, jako rozwiązywanie równania

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

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.

Przykład

Implementacja wyznaczania odwrotności liczby $a$ (bez dzielenia!) jest możliwa, gdy odpowiednią metodą będziemy poszukiwać rozwiązania równania

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

To zadanie jest ważne praktycznie, np. tak można poprawić precyzję funkcji wektorowych stosowanych w niektórych procesorach AMD. Okazuje się, że instrukcja procesora służąca do równoległego 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=\{f\in C([a,b]):f(a)\cdot f(b)<0\}$,

to znaczy $f\in F$ przyjmują w krańcach przedziału wartości przeciwnego znaku. Oczywiście, każda funkcja $f\in F$ ma, na mocy twierdzenia Darboux, 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 następnym kroku 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

---

lewy = a; prawy = b;
flewy = f(lewy); fprawy = f(prawy);
x = (a+b)/2; 	/* przybliżenie rozwiązania */
e = (b-a)/2; 	/* przedział lokalizujący rozwiązanie dokładne */
while (e > <math>\epsilon</math>) 
{
	fx = f(x);   	/* reszta */
	if ( abs(fx) == 0 ) /* trafiliśmy dokładnie w miejsce zerowe */
		return(x); 
	if ( sign(fx) != sign(flewy) ) /* tzn. f(lewy)*f(x) < 0 */
	{
		prawy = x;
		fprawy = fx;
	}
	else
	{
		lewy = x;
		flewy = fx;
	}
	x = (lewy+prawy)/2; /* najlepsze przybliżenie rozwiązania przy danym przedziale */
	e = e/2;
}
return(x);

Z konstrukcji metody łatwo wynika, że po wykonaniu $k$ obrotów pętli while (czyli po obliczeniu $k+2$ 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ą 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(ϵ)=\lceil {\log }_2\frac{(b-a)}{ϵ}\rceil -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. Dla większej ogólności, będziemy zakładać teraz, że $f:D\to R^N$ i $D$ jest otwartym, niepustym podzbiorem $R^N$.

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)$

Taki $x$, dla którego zachodzi powyższa równość, nazywamy punktem stałym odwzorowania $\varphi$.

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

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

Stefan Banach
Zobacz biografię

Twierdzenie Banacha, o kontrakcji

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

$x=\varphi (x)$

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

$\begin{array}{rl}\Vert x_k-x_{k-1}\Vert & =\Vert \varphi (x_{k-1})-\varphi (x_{k-2})\Vert \le L\Vert x_{k-1}-x_{k-2}\Vert \\ & \le & \cdots \le L^{k-1}\Vert x_1-x_0\Vert ,\end{array}$

dla $k\ge s$ mamy

$\begin{array}{rlr}\Vert x_k-x_s\Vert & & \le \sum _{j=s+1}^k\Vert x_j-x_{j-1}\Vert \le \sum _{j=s+1}^k{L}^{j-1}\Vert x_1-x_0\Vert \\ & & =L^s(1+L+\cdots +L^{k-s-1})\Vert x_1-x_0\Vert \le \frac{L^s}{1-L}\Vert x_1-x_0\Vert .\end{array}$

Ciąg $\{x_k{\}}_k$ jest więc ciągiem Cauchy'ego. Stąd istnieje granica $\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 (\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=\alpha$,

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

$\Vert \alpha -\beta \Vert =\Vert \varphi (\alpha )-\varphi (\beta )\Vert \le L\Vert \alpha -\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 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 równanie Keplera, gdy $0<ϵ<1$:

$x=M+ϵ\sin (x),\text{dla}x\in R$

Graficzna ilustracja równania Keplera dla $M=1$ i $ϵ=\frac{1}{4}$.

W tym przypadku $\varphi (x)=M+ϵ\sin (x)$. Zauważmy, że w funkcja $\varphi$ jest zwężająca ze stałą

$L\le \underset{x}{\max }|{\varphi }^{\prime }(x)|\le ϵ<1$

Ponieważ obrazem prostej przy przekształceniu $\varphi$ jest odcinek $D=[M-ϵ,M+ϵ]$, to znaczy, że $\varphi$ --- ograniczona do $D$ --- spełnia założenia twierdzenia Banacha o kontrakcji. Stąd istnieje dokładnie jedno rozwiązanie naszego równania w przedziale $D$. 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 D$. Jednak, gdy $ϵ\approx 1$, zbieżność może być bardzo powolna... (Wkrótce przekonasz się, że są szybsze metody).

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)<+\mathrm{\infty }$, 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 do czynienia, gdy rozpatrzymy metodę 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 $1/2$ 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.

Plik:Newton.jpg

Isaac Newton
Zobacz biografię

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

Algorytm Metoda Newtona (stycznych)

---

for k = 1,2,...
	<math>x_k\,=\,x_{k-1}\,-\,\frac{f(x_{k-1})}{f'(x_{k-1})}</math>;

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

<flash>file=Newtononestep.swf|width=550|height=300</flash>

Postęp iteracji Newtona

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.

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^{*}$.

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\text{ 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

$|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

$\begin{array}{rl}{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)!}}\\ & =(x_{k-1}-x^{*})(1-\frac{1}{m}\frac{f^{(m)}({\eta }_k^{(1)})}{f^{(m)}({\eta }_k^{(2)})})\\ & \approx & (x_{k-1}-x^{*})(1-\frac{1}{m}),\end{array}$

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

Metoda Newtona jest więc zbieżna lokalnie. Gdy $x_0$ jest zbyt daleko od rozwiązania może zdarzyć się, że iteracja Newtona zacznie nas oddalać od miejsca zerowego, co ilustruje poniższy przykład:

<flash>file=Newtononestepdiv.swf|width=550|height=300</flash>

Metoda Newtona: jeśli startujemy zbyt daleko od miejsca zerowego $f$, zamiast przybliżać się do niego, zaczynamy się oddalać! (gdzie będzie $x_3$?...)

Z powyższego 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^{*})|}{2|f^{\prime }(x^{*})|}$

Mówimy, że zbieżność metody Newtona, gdy $f^{\prime }(x^{*})\ne 0$ jest kwadratowa.

Stwierdzenie

Jeśli $f^{\prime }(x^{*})\ne 0$ oraz $f^{″}(x^{*})=0$ to zbieżność jest nawet szybsza. Z kolei dla zera $m$-krotnego (tzn. $f(x^{*})=f^{\prime }(x^{*})=\dots f^{(m)}(x^{*})=0$, $m>1$) 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.

Porównanie zbieżności metody bisekcji i stycznych dla równania $e^x-1=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 $\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 O rzędzie zbieżności metody Newtona

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)=(x-1{)}^5$ jest liniowa z ilorazem $\frac{4}{5}$ (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 $\frac{1}{2}$.

Metoda siecznych

Inną znaną i często używaną metodą iteracyjną, opartą na podobnym pomyśle linearyzacyjnym 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ż sieczna dla $f$ w punktach $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

Algorytm Metoda siecznych

---

for k = 1,2,...
	<math>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>;
end

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 przybliża pochodną $f$,

$\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$. Można pokazać, że przy podobnych założeniach o funkcji, wykładnik zbieżności metody siecznych dla zer jednokrotnych i 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.

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

Niewątpliwą zaletą metody siecznych jest jednak to, że nie wymaga obliczania pochodnej funkcji (bywa, że dokładne wyznaczenie pochodnej jest niemożliwe, gdy np. funkcja jest zadana zewnętrzną procedurą, do której kodu źródłowego nie mamy dostępu; zwykle też koszt obliczenia wartości pochodnej jest wyższy od kosztu obliczenia wartości funkcji). Jest to również istotne w pakietach numerycznych, gdzie czasem nie chcemy wymagać od użytkownika czegokolwiek, oprócz podania wzoru na funkcję 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 siecznych --- choć wolniej zbieżna niż metoda Newtona --- 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.

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.

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 np. jak w metodzie stycznych. W kolejnej iteracji wymieniamy jeden z punktów albo wedle metody siecznych (i wtedy zapewne szybciej zbliżamy się do zera), albo wykonując bisekcję (aby zagwarantować sobie, że w wiadomym przedziale miejsce zerowe rzeczywiście się znajduje).

Ten prosty pomysł metody hybrydowej wymaga jednak subtelnego dopracowania. Zostało to zrobione w 1973 roku przez Richarda Brenta. Funkcja MATLABa (i Octave'a) fzero implementują właśnie metodę Brenta. Autorem implementacji w Octave jest ówczesny student matematyki na Uniwersytecie Warszawskim, Łukasz Bodzon. Fortranowski kod metody Brenta można znaleźć także w Netlibie. Inną funkcją Octave'a służącą rozwiązywaniu równań nieliniowych jest fsolve:

octave:1> [X, MSG, INFO] = fsolve ('cos', 1) X = 1.5708 MSG = 1 INFO = solution converged within specified tolerance octave:2> cos(X) ans = 6.1230e-17

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

Niektóre z poznanych metod można łatwo rozszerzyć na przypadek układu $N$ równań z $N$ niewiadomymi, to znaczy

$F(x)=0$,

gdzie $F:R^N\to R^N$.

Metoda Banacha

Jak pamiętamy, metodę Banacha sformułowaliśmy od razu dla zagadnienia wielowymiarowego. Analiza i własności metody są zatem już omówione.

Wielowymiarowa metoda Newtona

Okazuje się, że metodę Newtona można uogólnić na przypadek układu $N$ równań nieliniowych z $N$ niewiadomymi. Zapiszmy wzór na skalarną metodę Newtona odrobinę inaczej:

$x_{k+1}=x_k-[F^{\prime }(x_k){]}^{-1}F(x_k)$

Niezwykłe jest, że taki wzór nie tylko ma sens w przypadku, gdy $F:R^N\to R^N$ (wtedy $F^{\prime }(x_k)$ jest macierzą Jakobianu $F$ w punkcie $x_k$), ale dodatkowo ta metoda zachowuje wszystkie własności metody stycznych dla przypadku skalarnego:

Twierdzenie O zbieżności wielowymiarowej metody Newtona

Załóżmy, że $F:R^N\to R^N$ i istnieje $x^{*}\in R^N$ taki, że

$F(x^{*})=0$

Załóżmy ponadto, że $F$ jest różniczkowalna, a jej pochodna $F^{\prime }:R^N\to R^{N\times N}$ jest lipschitzowska i dodatkowo

$F^{\prime }(x^{*})\text{ jest nieosobliwa}$

Wówczas, jeśli tylko $x_0$ jest dostatecznie blisko rozwiązania $x^{*}$, to ciąg kolejnych przybliżeń $x_k$, generowany wielowymiarową metodą Newtona, jest zbieżny do $x^{*}$. Co więcej, szybkość zbieżności jest kwadratowa.

Implementacja wielowymiarowej metody Newtona

Implementując wielowymiarową metodę Newtona, musimy dysponować nie tylko funkcją obliczającą $N$ współrzędnych wektora wartości $F$, ale także funkcją wyznaczającą $N^2$ elementów macierzy pochodnej $F$ w zadanym punkcie $x\in R^N$. Zwróćmy uwagę na to, że w implementacji metody nie trzeba wyznaczać $F^{\prime }(x_k{)}^{-1}$, tylko rozwiązać układ równań:

Algorytm Wielowymiarowa metoda Newtona

---

while (!stop)
{
	rozwiąż (względem <math>s</math>) układ równań liniowych <math>F'(x_k)\, s = -F(x_k)</math>;
	<math>x_{k+1}</math> = <math>x_k</math> + <math>s</math>;
}

O tym, jak skutecznie rozwiązywać układy równań liniowych, dowiesz się z kolejnych wykładów. Dowiesz się także, dlaczego nie należy w implementacji korzystać z wyznaczonej explicite macierzy odwrotnej do macierzy Jakobianu.

Literatura

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj rozdział 3 w

• D. Kincaid, W. Cheney Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.

Rozdziały 3.5 i 3.6 nie są obowiązkowe.

Wiele wariantów metod rozwiązywania układów równań nieliniowych jest przedstawionych w znakomitej monografii

• C.T.Kelley, Iterative Solution of Systems of Linear and Nonlinear Equations, SIAM, 1995.

Opis metody Brenta znajdziesz w książce

• R. Brent, Algorithms for Minimization Without Derivatives, Prentice-Hall, 1973.