Tekst źródłowy strony MN02LAB

<!-- 
Konwertowane  z pliku LaTeX przez latex2mediawiki, zob. http://www.ii.uj.edu.pl/&nbsp;pawlik1/latex2mediawiki.php.
Niezb�dne rozszerzenia i modyfikacje oryginalnego latex2mediawiki 
wprowadzi� przykry@mimuw.edu.pl
-->
 
=Równania nieliniowe skalarne=

{{powrot |Metody numeryczne | do strony głównej 
przedmiotu <strong>Metody numeryczne</strong>}}

<div class="mw-collapsible mw-made=collapsible mw-collapsed">
Oglądaj wskazówki i rozwiązania __SHOWALL__<br>
Ukryj wskazówki i rozwiązania __HIDEALL__
</div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie: Metoda Newtona może być zbieżna globalnie</span> 
<div class="exercise">

Wykaż, że jeśli <math>f</math> jest rosnąca i wypukła na <math>[a,b]</math> oraz dla <math>x^*\in [a,b]f(x^*) = 0</math>, to metoda Newtona startująca z <math>x_0>x^*</math> jest zbieżna.
</div></div>

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">  
Powtarzając rachunki z dowodu lokalnej zbieżności, możemy łatwo stwierdzić, że kolejne
iteracje dają <math>x_k</math>, które także spełniają warunek <math>x_k - x^* > 0</math>. Ponadto
nietrudno sprawdzić, że <math>x_{k+1} < x_k</math>, a więc mamy ciąg malejący i
ograniczony. Jego granicą musi być <math>x^*</math> (dlaczego?).
</div></div></div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie: Fraktale</span> 
<div class="exercise">

Ciekawy zbiór o charakterze fraktalnym powstaje w wyniku zastosowania metody
Newtona do rozwiązania równania <math>z^n-1=0</math> w dziedzinie zespolonej. Punkt <math>z_0 =
(x_0,y_0)</math> należy do <strong>basenu zbieżności metody</strong>, jeśli startująca z niego metoda Newtona jest zbieżna do (jakiegokolwiek) pierwiastka w/w równania. Kolor odpowiadający <math>z_0</math> 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.

[[Image:MNfraktal.png|thumb|550px|center|Basen zbieżności metody Newtona w okolicy początku układu
współrzędnych, dla równania <math>z^5-1=0</math>]]

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
<div style="font-size:smaller; background-color:#f9fff9; padding: 1em"> Bardzo wygodnie napisać taki program w postaci skryptu Octave. Pamiętaj jednak, by uniknąć w nim pętli w rodzaju

<code style="color: #006">for xpixels = 1:1024, for ypixels = 1:768,... wyznacz kolor pixela ....</code>

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! 
 </div>
</div></div>

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
<div style="font-size:smaller; background-color:#f9fff9; padding: 1em"> Kluczem do efektywnej implementacji w Octave jest przeprowadzenie iteracji Newtona na <strong>wszystkich pikselach jednocześnie</strong>. W tym celu musisz skorzystać z prowadzenia działań na całych macierzach. </div>
</div></div>

</div></div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie: Pierwiastkowanie</span> 
<div class="exercise">

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

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
<div style="font-size:smaller; background-color:#f9fff9; padding: 1em"> 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 <strong>tej konkretnej funkcji, <math>F(x) = x^2-a</math></strong>, będziesz miał szansę dokonać większych optymalizacji. </div>
</div></div>

</div></div>

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">  
Nasza iteracja, którą chcemy zaprogramować, to
<center><math>
x_{k+1} = \frac{1}{2}\left( x_k + \frac{a}{x_k}\right)</math></center>

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

<center><math>|x_{k+1}-x^*|\leq K |x_{k}-x^*|^2
</math></center>

co właśnie tłumaczy się jako podwajanie liczby dokładnych cyfr znaczących wyniku po każdej iteracji.

Ale jeśli <math>a=0</math>, 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 <math>x^*=0</math> będziemy dostawali co mniej więcej trzy kroki... Zbieżność będzie ''powolna!'' (Na szczęście, gdy <math>a=0</math>, 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 <math>x_0</math> na pewno jest "dostatecznie blisko" rozwiązania?

</div></div></div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie: Odwrotność bez dzielenia</span> 
<div class="exercise">

Aby wyznaczyć <math>1/a</math> dla <math>a>0</math> bez dzielenia(!), można zastosować metodę Newtona
do funkcji <math>f(x) = 1/x - a</math>. Pokaż, że na <math>k</math>-tym kroku iteracji,
<center><math>
|ax_k - 1| = |ax_{k-1} - 1|^2</math></center>

Dla jakich <math>x_0</math> metoda będzie zbieżna do <math>\frac{1}{a}</math>, a dla jakich nie?
Oceń, ile iteracji potrzeba do spełnienia warunku <math>\frac{|x_k-\frac{1}{a}|}{|a|} \leq
\epsilon</math>, gdy <math>\frac{|x_0-\frac{1}{a}|}{|a|} \leq \gamma < 1</math>,
</div></div>

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">  
Proste ćwiczenie z [[Analiza_matematyczna_1|analizy matematycznej]]. Warunkiem zbieżności jest oczywiście (ze względu na
równość w wyrażeniu na błąd iteracji powyżej) <math>|ax_0 - 1| < 1</math>, to znaczy 
<center><math>0 < x_0 < \frac{2}{a}</math>.</center>

</div></div></div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie</span> 
<div class="exercise">

Zaimplementuj metodę bisekcji. Sprawdź, jak będzie działać m.in. dla funkcji
* <math>f(x) = \sin(x)</math>,
* <math>f(x) = \sin^2(x)</math>,
* <math>f(x) = (x-1)^5</math> (wzór A),
* <math>f(x) = (x-1)\cdot (x^4 - 4x^3 + 6x^2 - 4x + 1)</math> (wzór B),
* <math>f(x) = (x-2)^{13}</math> (wzór C),
* <math>f(x) = x^{13} - 26x^{12} + 312x^{11} - 2288x^{10} + ... - 8192</math> (wzór D),
 
gdy tolerancję błędu zadasz na poziomie <math>10^{-10}</math>.

[[Image:MNbisekcjawblad.png|thumb|550px|center|Metoda bisekcji ma kłopoty, gdy funkcja zadana jest
wzorem D.]]
[[Image:MNbisekcjawresid.png|thumb|550px|center|Residuum też jest duże, gdy <math>f</math> zadana jest
wzorem D.]]

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

</div></div>

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">  

Jeśli nie pomylisz się, metoda powinna zbiegać bez problemów do zera funkcji
<math>\sin(x)</math>, dawać komunikat o błędzie dla <math>\sin^2(x)</math> (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 <strong>w ogóle nie zawierają</strong> 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 <strong>działa</strong> 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 <strong>przypuszczać</strong>, że implementacja jest poprawna....

</div></div></div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie</span> 
<div class="exercise">

Wskaż <strong>wszystkie</strong> wartości <math>x_0</math>, dla jakich metoda Newtona będzie zbieżna
do rozwiązania <math>x^*</math> równania

<center><math>\arctan(x) = 0</math></center>

Wyznacz wartość <math>X_0</math>, z którego startując powinieneś dostać ciąg oscylujący <math>X_0, -X_0,\ldots</math>.
Sprawdź eksperymentalnie, czy tak jest rzeczywiście.

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
<div style="font-size:smaller; background-color:#f9fff9; padding: 1em"> Aby znaleźć graniczne <math>X_0</math>, czyli takie, dla którego dostaniesz oscylacje <math>X_0, -X_0,\ldots</math>,  musisz rozwiązać równanie nieliniowe:

<center><math>-X_0 = X_0 - \frac{\arctan(X_0)}{\frac{1}{1+X_0^2}}</math></center>

To łatwe.
 </div>
</div></div>

</div></div>

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">  
Dla <math>-X_0 < x_0 < X_0</math> mamy zbieżność. Dla <math>|x_0| = |X_0|</math> oscylacje, dla
większych wartości --- rozbieżność. Jednak w eksperymencie wszystko zależy od
tego, jak dokładnie wyznaczyliśmy <math>X_0</math>. Na przykład, funkcja <code style="color: #006">fzero</code> w Octave
wyznaczyła (przy standardowych ustawieniach) numeryczne przybliżenie <math>X_0</math>, dla którgo metoda Newtona dała następujący ciąg:
<div style="font-family: monospace; white-space: pre; border-style: dashed; border-width: thin; border-color: black; margin: 1em; padding:1em; color: #444; background-color:#fdfdfd;"><nowiki>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
</nowiki></div>
 
</div></div></div>