MN02LAB

From Studia Informatyczne


Równania nieliniowe skalarne

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

Oglądaj wskazówki i rozwiązania
Ukryj wskazówki i rozwiązania

Ć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
Enlarge
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.
Enlarge
Metoda bisekcji ma kłopoty, gdy funkcja zadana jest wzorem D.
Residuum też jest duże, gdy  zadana jest wzorem D.
Enlarge
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