MN02LAB

Równania nieliniowe skalarne

Ćwiczenie: Metoda Newtona może być zbieżna globalnie

Wykaż, że jeśli $\displaystyle f$ jest rosnąca i wypukła na $\displaystyle [a,b]$ oraz dla $\displaystyle x^{*}\in [a,b]\displaystyle f(x^{*})=0$, to metoda Newtona startująca z $\displaystyle x_{0}>x^{*}$ jest zbieżna.

Rozwiązanie

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

Ćwiczenie: Fraktale

Ciekawy zbiór o charakterze fraktalnym powstaje w wyniku zastosowania metody Newtona do rozwiązania równania $\displaystyle z^{n}-1=0$ w dziedzinie zespolonej. Punkt $\displaystyle z_{0}=(x_{0},y_{0})$ należy do basenu zbieżności metody, jeśli startująca z niego metoda Newtona jest zbieżna do (jakiegokolwiek) pierwiastka w/w równania. Kolor odpowiadający $\displaystyle 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, takie jak na rysunku poniżej.

Basen zbieżności metody Newtona w okolicy początku układu współrzędnych, dla równania $\displaystyle z^{5}-1=0$

Wskazówka

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

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!

Wskazówka

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

Ćwiczenie: Pierwiastkowanie

Niech $\displaystyle a\geq 0$. Aby wyznaczyć $\displaystyle {\sqrt {a}}$, można skorzystać z metody Newtona dla równania $\displaystyle 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 $\displaystyle a$, dla którego to nie będzie prawdą.

Wskazówka

Na pewno kusi Cię, by zaprogramować najpierw ogólnie funkcję "metoda Newtona", a następnie przekazać jej jako parametr naszą funkcję. Jednak gdy prostu wyliczysz na karteczce wzór na iterację dla tej konkretnej funkcji, $\displaystyle F(x)=x^{2}-a$, będziesz miał szansę dokonać większych optymalizacji.

Rozwiązanie

Nasza iteracja, którą chcemy zaprogramować, to

$\displaystyle x_{k+1}={\frac {1}{2}}\left(x_{k}+{\frac {a}{x_{k}}}\right).$

Jest to, znana już w starożytności, metoda Herona. Teraz jest jasne, dlaczego zbiega tak szybko do $\displaystyle {\sqrt {a}}$. Dla $\displaystyle a\neq 0$, spełnione są założenia twierdzenia o zbieżności metody Newtona, więc

$\displaystyle |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 po każdej iteracji.

Ale jeśli $\displaystyle a=0$, to metoda zbieżna jest już tylko liniowo (z ilorazem 0.5 --- bo jaki wtedy jest wzór na iterację?); dlatego kolejną dokładną cyfrę wyniku $\displaystyle x^{*}=0$ będziemy dostawali co mniej więcej trzy kroki... Zbieżność będzie powolna! (Na szczęście, gdy $\displaystyle a=0$, nie ma żadnego problemu z podaniem dokładnego rozwiązania...)

Nie dyskutujemy tutaj sposobu doboru właściwego punktu startowego dla metody: na razie wiadomo, że musi być "dostatecznie blisko" rozwiązania. Czy jesteś w stanie sprawdzić --- wypisując bardziej precyzyjne oszacowanie błędu (zob. dowód twierdzenia o zbieżności) --- jakie $\displaystyle x_{0}$ na pewno jest "dostatecznie blisko" rozwiązania?

Ćwiczenie: Odwrotność bez dzielenia

Aby wyznaczyć $\displaystyle 1/a$ dla $\displaystyle a>0$ bez dzielenia(!), można zastosować metodę Newtona do funkcji $\displaystyle f(x)=1/x-a$. Pokaż, że na $\displaystyle k$-tym kroku iteracji,

$\displaystyle |ax_{k}-1|=|ax_{k-1}-1|^{2}.$

Dla jakich $\displaystyle x_{0}$ metoda będzie zbieżna do $\displaystyle {\frac {1}{a}}$, a dla jakich nie? Oceń, ile iteracji potrzeba do spełnienia warunku $\displaystyle {\frac {|x_{k}-{\frac {1}{a}}|}{|a|}}\leq \epsilon$, gdy $\displaystyle {\frac {|x_{0}-{\frac {1}{a}}|}{|a|}}\leq \gamma <1$,

Rozwiązanie

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

$\displaystyle 0<x_{0}<{\frac {2}{a}}.$

Ćwiczenie

Zaimplementuj metodę bisekcji. Sprawdź, jak będzie działać m.in. dla funkcji

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

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

Metoda bisekcji ma kłopoty, gdy funkcja zadana jest wzorem D.

Residuum też jest duże, gdy $\displaystyle f$ zadana jest wzorem D.

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

Rozwiązanie

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

Zapewne zauważyłeś, ż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 w ogóle nie zawierają prawdziwego miejsca zerowego.

Wyjaśnieniem tego fenomenu jest zjawisko redukcji cyfr przy odejmowaniu, o którym mowa w następnym wykładzie.

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

Tak więc, można przypuszczać, że implementacja jest poprawna....

Ćwiczenie

Wskaż wszystkie wartości $\displaystyle x_{0}$, dla jakich metoda Newtona będzie zbieżna do rozwiązania $\displaystyle x^{*}$ równania

$\displaystyle \arctan(x)=0.$

Wyznacz wartość $\displaystyle X_{0}$, z którego startując powinieneś dostać ciąg oscylujący $\displaystyle X_{0},-X_{0},\ldots$. Sprawdź eksperymentalnie, czy tak jest rzeczywiście.

Wskazówka

Aby znaleźć graniczne $\displaystyle X_{0}$, czyli takie, dla którego dostaniesz oscylacje $\displaystyle X_{0},-X_{0},\ldots$, musisz rozwiązać równanie nieliniowe: $\displaystyle -X_{0}=X_{0}-{\frac {\arctan(X_{0})}{\frac {1}{1+X_{0}^{2}}}}.$

To łatwe.

Rozwiązanie

Dla $\displaystyle -X_{0}<x_{0}<X_{0}$ mamy zbieżność. Dla $\displaystyle |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 $\displaystyle X_{0}$. Na przykład, funkcja fzero w Octave wyznaczyła (przy standardowych ustawieniach) numeryczne przybliżenie $\displaystyle X_{0}$, dla którgo metoda Newtona dała następujący ciąg:

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