MN02LAB

Równania nieliniowe skalarne

Ćwiczenie: 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.

Rozwiązanie

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

Ćwiczenie: 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 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.

Basen zbiezności metody Newtona w okolicy początku układu współrzędnych, dla równania $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 ....
	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!

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 $a\ge 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ą.

Wskazówka

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

Rozwiązanie

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

$x_{k+1}=\frac{1}{2}(x_k+\frac{a}{x_k}).$

Dla $a\ne 0$, spełnione są założenia twierdzenia o zbieżności metody Newtona, więc

$|x_{k+1}-x^{*}|\le 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$.

Ćwiczenie: 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,

$|a{x}_k-1|=|a{x}_{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|}\le ϵ$, gdy $\frac{|x_0-\frac{1}{a}|}{|a|}\le \gamma <1$,

Rozwiązanie

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) $|a{x}_0-1|<1$, to znaczy

$0<x_0<\frac{2}{a}.$

Ćwiczenie

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

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

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

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

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

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

Rozwiązanie

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 w ogóle nie zawierają prawdziwego miejsca zerowego.

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

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 przypuszczać, że implementacja jest poprawna.

Ćwiczenie

Wskaż 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,\dots$. Sprawdź eksperymentalnie, czy tak jest rzeczywiście.

Wskazówka

Aby znaleźć graniczne $X_0$, czyli takie, dla którego dostaniesz

oscylacje $X_0,-X_0,\dots$, musisz rozwiązać równanie nieliniowe:

$-X_0=X_0-\frac{\arctan (X_0)}{\frac{1}{1+X_0^2}}.$

To łatwe.

Rozwiązanie

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:

 
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