Rozwiązywanie równań nieliniowych

W poniższych zadaniach przetestujemy jak działają cztery metody rozwiązywania równań nieliniowych skalarnych, tzn. metoda bisekcji, metoda Newtona, metoda siecznych i metoda iteracji prostych oraz
dwie metody rozwiązywania układów równań nieliniowych, czyli wielowymiarowa metoda Newtona i metoda iteracji prostych w najprostszej formie.

Zapoznamy się również z funkcjami octave'a fzero() oraz fsolve() służącymi do rozwiązywania równań nieliniowych skalarnych i układów równań nieliniowych.

  1. (#)
    Zaimplementuj metodę bisekcji w skrypcie - dla rozwiązania równania $ \cos(x)=0 $ na odcinku $ [0,2] $. Przetestuj, czy funkcja znajdzie dobre przybliżenie $ \pi/2 $.

  2. Przetestuj metodę bisekcji z poprzedniego zadania do rozwiązania innych problemów:
    \begin{eqnarray*}
    x^2-2=0 \\ \exp(x)=2 \\ \cos(10*x)=0
    \end{eqnarray*}

    startując z odcinka $ [0,110] $.

  3. Napisz funkcję octave'a:\\
     function [y,iter,kod]=bisekcja(FCN,a,b,...
                                    tola,maxit),

    której parametrami będą

    • wskaźnik do funkcji FCN (inaczej: function handle),
    • $ a,b $ - końce przedziału,
    • tola - żądana tolerancja bezwzględna - błąd powinien być mniejszy od tola,
    • maxit - maksymalna ilość iteracji.

    Funkcja ma zwrócić

    • y - przybliżone rozwiązanie,
    • iter - ilość iteracji,
    • kod - kod wyniku:
      • $ 0 $ metoda zbiegła,
      • $ 1 $ - metoda zatrzymała się z powodu przekroczenia maksymalnej ilości iteracji,
      • $ 2 $ - wartości w końcach odcinka funkcji mają ten sam znak.

    Funkcja ma działać również jeśli podamy tylko trzy lub cztery pierwsze argumenty - wtedy maksymalna ilość iteracji domyślnie powinna wynosić sto, a tolerancja $ 10^{-5} $.

  4. Zapoznaj się z pomocą dla funkcji octave'a fzero() - rozwiązującą skalarne równania nieliniowe, przetestuj ją na kilku prostych przykładach skalarnych np.
    • $ \cos(x)=0 $,
    • $ x^2-2=0 $,
    • $ x^{10}-2=0 $,
    • $ x-\sin(x)=0 $

    i innych prostych równaniach nieliniowych.

    Proszę przetestować tę funkcję dla różnych wartości początkowego przedziału $ x_0=[a,b] $.

  5. Zapoznaj się z pomocą dla funkcji octave'a fsolve() - rozwiązującą układy równań nieliniowych. Przetestuj ją na kilku prostych przykładach skalarnych, np. na tych z poprzedniego zadania.
  6. (#)Zaimplementuj metodę Newtona w octavie i przetestuj jej zbieżność dla następujących funkcji:
      [ f1(x)=]
    1. $ x*x-2 $ z $ x_0=2 $,
    2. $ x*x*x-27 $ z $ x_0=27 $,
    3. $ \exp(x)-2 $ z $ x_0=10,-10,1e3 $,
    4. $ \sin(x) $ dla $ x_0=2 $,
    5. $ \cos(x) $ z $ x_0=2 $
    6. $ (x-2)^k $ x $ x_0=3 $ dla $  k=2,4,8,16 $,
    7. $ x*x-2 $ z $ x_0=10^6 $, sprawdź, czy metoda Newtona zbiega, a jeśli tak - to czy zbiega ona kwadratowo,
    8. $ 1/x-a $ dla danych $ a=0.5,2,4,100 $ (oczywiście implementując bez dzielenia) jak dobrać $ x_0 $?

    Dla wszystkich tych funkcji znamy rozwiązania, więc można wyświetlać równocześnie na ekranie błąd

    $$e_n=x_n-r$$

    (tutaj $ r $ - znane rozwiązanie) i badać eksperymentalnie rząd zbieżności, tzn. drukować na ekranie stosunki błędu bieżącego do poprzedniego w odpowiedniej potędze:

    $$|e_n|/|e_{n-1}|^p$$

    dla $ p=1,2,3 $.

    Proszę powtórzyć testy z innymi różnymi wartościami startowymi $ x_0 $.

  7. Powtórz zadanie  [link] zastępując metodę Newtona przez metodę siecznych.

    Przetestuj, czy

    $$e_n/(e_{n-1}e_{n-2})$$

    asymptotycznie zbiega do stałej, i czy

    $$|e_n|/|e_{n-1}|^p$$

    dla $ p=(1+ \sqrt{5})/2 $ zbiega do stałej np. dla $ x*x-2 $ z $ x_0=2 $ lub innych równań z poprzedniego zadania.

    Za $ x_1 $ możemy przyjąć $ x_0+10^{-7} $.

  8. Sprawdź, czy metoda iteracji prostych
    $$x_n=\cos(x_{n-1})$$

    zbiega do $ x^*=\cos(x^*) $.

    Zbadaj eksperymentalnie, czy zbieżność jest liniowa, tzn. czy

    $$\frac{|x^*-x_n|}{|x^*-x_{n-1}|}$$

    zbiega do stałej. Za dobre przybliżenie $ x^* $ można przyjąć rozwiązanie równania obliczone za pomocą wywołania funkcji octave'a: fzero().

  9. Zaimplementuj przybliżoną metodę Newtona, w której pochodną przybliżamy ilorazem różnicowym, tzn.
    $$<br />
  x_{n+1}=x_n - \frac{f(x_n)}{\frac{f(x_n+h)-f(x_n)}{h}}=<br />
   x_n - \frac{f(x_n)*h}{f(x_n+h)-f(x_n)}<br />
$$

    dla ustalonego $ h $.

    Przetestuj różne $ h $, np. $ h=10^{-4},10^{-7},10^{-10} $ itp.

    Porównać zbieżność z dokładną metodą Newtona (szczególnie ostatnie iteracje) dla funkcji z zadania  [link].

  10. \textrm{Rozwikływanie funkcji:}

    Dla funkcji $ y(x) $ zadanej w sposób uwikłany równaniem

    $$g(x,y)=2x^2+3y^2-3=0$$

    znajdź przybliżone wartości $ y_k=y(x_k) $ dla $ x_k=k*h $ dla $ k=-N,\ldots,N $ i $ h=1/N $, tzn. spełniające $ g(x_k,y_k)=0 $.

    Testuj dla różnych wartości $ N $, np. $ N= 10,20,40,80,... $.

    oblicz $ y_k $ rozwiązując równanie

    $$g(x_k,y_k)=0.$$

    korzystając z odpowiedniej funkcji octave'a lub własnej implementacji metody Newtona.

    Jak w kolejnych krokach dobierać przybliżenie startowe w metodzie Newtona?

  11. \textrm{Odwracanie funkcji}

    Rozpatrzmy daną funkcję, np.

    $$f(x)=\sin(x)+2*x.$$

    Znajdź wartości funkcji odwrotnej $ f^{-1} $ na odcinku $ [0,5] $ na siatce $ k*h $ dla $ k=0,..,100 $. Narysuj wykres funkcji $ f^{-1} $.

    Sam wykres można narysować dużo prościej bez wyliczania wartości $ f^{-1} $. Jak to zrobić w octave'ie?

  12. (#) \textrm{Wielowymiarowa metoda Newtona}

    Zaimplementuj wielowymiarową metodę Newtona w wersji z dokładnym Jakobianem i w wersji, gdy Jakobian przybliżamy różnicami dzielonymi z parametrem $ h $.

    Zastosuj tę metodę do rozwiązania układu
    \begin{eqnarray*}
    f_1(x_1,x_2)&=&x_1+2x_2=1\\
    f_2(x_1,x_2)&=&3*x_1^2+x_2^2=1
    \end{eqnarray*}
    dla różnych przybliżeń początkowych.

    W przypadku zbieżności policz

    $$\frac{\|x-x_n\|_p}{\|x-x_{n-1}\|_p^q}$$

    dla $ p=1,2,\infty $ oraz $ q=1,2,3 $. Czy ten iloraz zbiega
    do stałej dla jakiegoś z tych $ q $?

    Funkcję $ F(x):=(f_1(x),f_2(x))^T $ w octavie można zdefiniować następująco:

    function y=F(x)
    #x- wektor pionowy

       y=[[1 2]*x-1;3.*x(1)*x(1)+x(2)*x(2)-1];

    endfunction

  13. \textrm{Wielowymiarowa metoda iteracji prostych}

    Zastosuj metodę iteracji prostych

    $$x_n=x_n-aF(x)$$

    do układu z zadania  [link], tzn. przyjmijmy, że $ F(x):=(f_1(x),f_2(x))^T $, a parametr $ a\not=0 $ np. $ a=1,-1,10^{-1} $ itp.

    Czy metoda zawsze zbiega, a jeśli tak - to jak szybko?
    W przypadku zbieżności policz

    $$\frac{\|x-x_n\|_p}{\|x-x_{n-1}\|_p}$$

    dla $ p=1,2,\infty $.

  14. Do układu równań z zadania  [link] zastosuj funkcję fsolve() - porównaj z implementacją metody Newtona w zadaniu  [link], która była zastosowana do rozwiązania tego układu.

    Czy obie metody zbiegają do tego samego rozwiązania dla tych samych wartości wektorów startowych?

    Sprawdź - za pomocą funkcji tic i toc - czas potrzebny do rozwiązania tego układu równań za pomocą obu metod.