MN03LAB

\section*{Ćwiczenia: równania nieliniowe skalarne}

\begin{exe}[Metoda Newtona może być zbieżna globalnie] Wykaż, że jeśli $f$ jest rosnąca i wypukła na $[a,b]$ oraz dla $x^*\in [a,b]$ $f(x^*) = 0$, to metoda Newtona startująca z $x_0>x^*$ jest zbieżna. \end{exe}

\begin{sol} Powtarzając rachunki z dowodu lokalnej zbieżności, możemy łatwo stwierdzić, że kolejne iteracje dają $x_k$ które także spełniają warunek $x_k - x^* > 0$. Ponadto nietrudno sprawdzić, że $x_{k+1} < x_k$, a więc mamy ciąg malejący i ograniczony. Jego granicą musi być $x^*$ (dlaczego?). \end{sol}

\begin{exe}[Fraktale]

Ciekawy zbiór o charakterze fraktalnym powstaje w wyniku zastosowania metody Newtona do rozwiązania równania $z^n-1=0$ w dziedzinie zespolonej. Punkt $z_0 = (x_0,y_0)$ należy do {\em basenu zbiezności metody}, jeśli metoda Newtona jest zeń zbieżna do (jakiegokolwiek) pierwiastka w/w równania. Kolor odpowiadający $z_0$ jest określany na podstawie liczby iteracji potrzebnych metodzie do zbieżności.

Zupełnie miłym (i estetycznie wartościowym) doświadczeniem, jest napisanie programu w Octave, który wyświetla baseny zbieżności metody Newtona jak poniżej.

\rysunek{fraktal.png}{Basen zbiezności metody Newtona w okolicy początku układu współrzędnych, dla równania $z^5-1=0$}

\hint{Bardzo wygodnie napisać taki program w postaci skryptu Octave. Pamiętaj jednak, by uniknąć w nim pętli w rodzaju

for xpixels = 1:1024\\ for ypixels = 1:768\\ ... wyznacz kolor pixela ....\\ end\\ end\\

Taka podwójnie zagnieżdżona pętla może skutecznie zniechęcić Cię do Octave'a --- obrazek będzie liczyć się wieki --- zupełnie niepotrzebnie! }

\hint{Kluczem do efektywnej implementacji w Octave jest przeprowadzenie iteracji Newtona na {\em wszystkich pikselach jednocześnie}. W tym celu musisz skorzystać z prowadzenia działań na całych macierzach.}

\end{exe}

\begin{exe}[Pierwiastkowanie] Niech $a\geq 0$. Aby wyznaczyć $\sqrt{a}$, można skorzystać z metody Newtona dla równania $x^2 = a$. Zaprogramuj tę metodę i sprawdź, jak wiele dokładnych cyfr wyniku dostajesz w kolejnych krokach. Czy to możliwe, by liczba dokładnych cyfr znaczących z grubsza podwajała się na każdej iteracji? Wskaż takie $a$ dla którego to nie będzie prawdą.

\hint{Na pewno kusi Cię, by zaprogramować najpierw ogólnie funkcję ,,metoda Newtona, a następnie przekazać jej jako parametr naszą funkcję. To oczywiście będzie działać, i to praktycznie tak samo szybko jak drugi spsoób, w którym po prostu wyliczysz na karteczce wzór na iterację dla {\em tej konkretnej funkcji, $F(x) = x^2-a$} i potem go zaimplementujesz. Jednak drugie podejście jest tutaj godne polecenia, bo ma na sobie znak charakterystyczny dla numeryki: niech wszystko, co programujesz będzie tak proste, jak tylko możliwe (ale nie prostsze)!} \end{exe}

\begin{sol} Nasza iteracja, którą chcemy zaprogramować, to $$ x_{k+1} = \frac{1}{2}\left( x_k + \frac{a}{x_k}\right). $$

Dla $a\neq 0$, spełnione są założenia twierdzenia o zbieżności metody Newtona, więc \[ |x_{k+1}-x^*|\leq K |x_{k}-x^*|^2 \] co właśnie tłumaczy się jako podwajanie liczby dokładnych cyfr znaczących wyniku.

Ale jeśli $a=0$ to metoda zbieżna jest już tylko liniowo (z ilorazem 0.5 --- bo jaki wtedy jest wzór na iterację?), więc co mniej więcej trzy kroki będziemy dostawali kolejne zero dokładnego wyniku $x^*=0$.

\end{sol}

\begin{exe}[Odwrotność bez dzielenia] Aby wyznaczyć $1/a$ dla $a>0$ bez dzielenia(!), można zastosować metodę Newtona do funkcji $f(x) = 1/x - a$. Pokaż, że na $k$-tym kroku iteracji, $$ |ax_k - 1| = |ax_{k-1} - 1|^2. $$ Dla jakich $x_0$ metoda będzie zbieżna do $\frac{1}{a}$, a dla jakich nie? Oceń, ile iteracji potrzeba do spełnienia warunku $\frac{|x_k-\frac{1}{a}|}{|a|} \leq \epsilon$, gdy $\frac{|x_0-\frac{1}{a}|}{|a|} \leq \gamma < 1$, \end{exe}

\begin{sol} Proste ćwiczenie z analizy. Warunkiem zbieżności jest oczywiście (ze względu na równość w wyrażeniu na błąd iteracji powyżej) $|ax_0 - 1| < 1$, to znaczy $$0 < x_0 < \frac{2}{a}.$$ \end{sol}

\begin{exe} Zaimplementuj metodę bisekcji. Sprawdź, jak będzie działać m.in. dla funkcji

\begin{enumerate} \item $f(x) = \sin(x)$, \item $f(x) = \sin^2(x)$, \item $f(x) = (x-1)^5$ (wzór A), \item $f(x) = (x-1)*(x^4 - 4x^3 + 6x^2 - 4x + 1)$ (wzór B), \item $f(x) = (x-2)^13$ (wzór C), \item $f(x) = x^13 - 26*x^12 + 312*x^11 - 2288*x^10 + ... - 8192$ (wzór D), \end{enumerate}

gdy tolerancję błędu zadasz na poziomie $10^{-10}$.

Jak wyjaśnić te wyniki? Czy możesz już być pewna, że masz dobrą implementację?

\end{exe}

\begin{sol}

Jeśli nie pomylisz się, metoda powinna zbiegać bez problemów do zera funkcji $\sin$, dawać komunikat o błędzie dla $\sin^2$ (bo nie ma przeciwnych znaków).

Zapewne zauważyłaś, że wzory A i B są matematycznie równoważne, podobnie zresztą jak wzory C i D. Niestety, tylko wzór A i C dają w efekcie dobre przybliżenia miejsca zerowego, podczas gdy wzory B i D prowadzą do oszacowania na miejsce zerowe, które {\em w ogóle nie zawierają} prawdziwego miejsca zerowego.

Wyjaśnieniem tego fenomenu jest zjawisko redukcji cyfr przy odejmowaniu.

\rysunek{bisekcjawblad.png}{Metoda bisekcji ma kłopoty, gdy funkcja zadana jest wzorem D.} \rysunek{bisekcjawresid.png}{Residuum też jest duże, gdy $f$ zadana jest wzorem D.}

Jeśli chodzi o pewność... No cóż, sprawdziłaś, że działa w przypadkach, gdy spodziewałaś się, że będzie działać. Że tam, gdzie spodziewałaś się kłopotów, lub komunikatu o błędzie --- tak rzeczywiście się stało. Wreszcie, potwierdziłaś, że zachowanie metody jest zgodne z jej znanymi właściwościami.

Tak więc, można {\em przypuszczać}, że implementacja jest poprawna.

\end{sol}

\begin{exe} Wskaż {\em wszystkie} wartości $x_0$, dla jakich metoda Newtona będzie zbieżna do rozwiązania $x^*$ równania \[ \arctan(x) = 0. \] Wyznacz wartość $X_0$, z którego startując powinieneś dostać ciąg oscylujący $X_0, -X_0,\ldots$. Sprawdź eksperymentalnie, czy tak jest rzeczywiście.

\hint{Aby znaleźć graniczne $X_0$, czyli takie, dla którego dostaniesz oscylacje $X_0, -X_0,\ldots$, musisz rozwiązać równanie nieliniowe: \[ -X_0 = X_0 - \frac{\arctan(X_0)}{\frac{1}{1+X_0^2}}. \] To łatwe. }

\end{exe}

\begin{sol} Dla $-X_0 < x_0 < X_0$ mamy zbieżność. Dla $|x_0| = |X_0|$ oscylacje, dla większych wartości --- rozbieżność. Jednak w eksperymencie wszystko zależy od tego, jak dokładnie wyznaczyliśmy $X_0$. Na przykład, mój solver w Octave wyznaczył wartość $X_0$, dla której metoda Newtona dała następujący ciąg: \begin{outoct} Iteracja: 1, x = -1.39175 Iteracja: 2, x = 1.39175 Iteracja: 3, x = -1.39175 Iteracja: 4, x = 1.39175 Iteracja: 5, x = -1.39175

.. itd ..

Iteracja: 25, x = -1.39176 Iteracja: 26, x = 1.39178 Iteracja: 27, x = -1.39183 Iteracja: 28, x = 1.39197 Iteracja: 29, x = -1.39235 Iteracja: 30, x = 1.39333 Iteracja: 31, x = -1.39594 Iteracja: 32, x = 1.40283 Iteracja: 33, x = -1.42117 Iteracja: 34, x = 1.47061 Iteracja: 35, x = -1.60867 Iteracja: 36, x = 2.03161 Iteracja: 37, x = -3.67722 Iteracja: 38, x = 15.2779 Iteracja: 39, x = -337.619 Iteracja: 40, x = 178376 Iteracja: 41, x = -4.99792e+10 Iteracja: 42, x = 3.92372e+21 Iteracja: 43, x = -2.41834e+43 Iteracja: 44, x = 9.18658e+86 Iteracja: 45, x = -1.32565e+174 Iteracja: 46, x = inf Iteracja: 47, x = nan \end{outoct}

\end{sol}