MN02LAB

Równania nieliniowe skalarne

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

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

Ć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),
• Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle \displaystyle f(x) = (x-2)^{13}} (wzór C),
• Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle \displaystyle f(x) = x^{13} - 26x^{12} + 312x^{11} - 2288x^{10} + ... - 8192} (wzór D),

gdy tolerancję błędu zadasz na poziomie Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle \displaystyle 10^{-10}} .

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

Residuum też jest duże, gdy Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle \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 Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle \displaystyle \sin(x)} , dawać komunikat o błędzie dla Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle \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 Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle \displaystyle x_0} , dla jakich metoda Newtona będzie zbieżna do rozwiązania Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle \displaystyle x^*} równania

Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle \displaystyle \arctan(x) = 0. }

Wyznacz wartość Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle \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