Matematyka Obliczeniowa

W tym skrypcie omówimy, jak implementować i testować metody numeryczne przy wykorzystaniu środowiska octave, czyli pakietu obliczeń numeryczno-naukowych. Chodzi nam o metody, które są omawiane w czasie standardowego semestralnego wykładu z Matematyki Obliczeniowej na wydziale Matematyki, Informatyki i Mechaniki Uniwersytetu Warszawskiego.

Głównym celem zajęć w laboratorium komputerowym, odbywających się równolegle z wykładem i ćwiczeniami z Matematyki Obliczeniowej, jest przetestowanie metod, które są zaimplementowane w pakiecie octave. Omawiamy też implementację metod rozwiązywania numerycznych zadań z wykorzystaniem octave'a. Wykorzystanie pakietu/środowiska matlab jest alternatywą dla pakietu octave. Pakiet ten jest komercyjny i w dużym stopniu kompatybilny z octavem. Szczegółowiej omawiamy kwestię pakietów w rozdziale  [link]. Zapoznanie się ze środowiskiem octave'a lub matlaba jest ubocznym celem laboratorium.

Warto dodać, że zdecydowana większość omawianych w tym skrypcie kodów octave'a będzie działała również w matlabie.

Wszystkie zadania i algorytmy omawiane w tym skrypcie są opisane w podstawowych podręcznikach do matematyki obliczeniowej.

Na wydziale Matematyki, Informatyki i Mechaniki Uniwersytetu Warszawskiego w semestrze letnim 2010/11 podstawowym podręcznikiem do wykładu z Matematyki Obliczeniowej była książka [Kincaid:2006:AN]. W tej pozycji można znaleźć również sporo zadań komputerowych, które mogą stanowić ciekawe uzupełnienie tego skryptu.

Metody numeryczne omawiane w tym skrypcie są szczegółowo omawiane w innych podręcznikach dotyczących analizy numerycznej, czyli inaczej matematykie obliczeniowej. Spośród nich warto również polecić: [Jankowska:1981:MN],[Jankowska:1982:MN], [Bjorck:1983:MN], [Ralston:1983:WAN] i [Stoer:1980:WMN]. Oczywiście w języku angielskim opublikowano wiele podręczników do matematyki obliczeniowej. Niektóre z nich zawierają zadania komputerowe. Szczególnie polecam książki: [Quarteroni:2003:SCM], [Quarteroni:2000:NM] zawierające więcej materiału, niż obejmuje program standardowego kursu Matematyki Obliczeniowej na wydziale Matematyki, Informatyki i Mechaniki Uniwersytetu Warszawskiego. Warto też polecić inne podręczniki i monografie: [Golub:1993:SC], [Golub:1992:SCD], [Demmel:1997:ANL], [Ortega:2000:ISN], [Ortega:1990:NA], [Kelley:2003:SNE], [Trefethen:1997:NLA], [Kielbasinski:1992:NAL], [Hackbusch:1994:ISL] i inne. Zazwyczaj treść tych książek tylko w jakiejś części pokrywa się z zakresem materiału przedstawionym w tym skrypcie. Co do samego octave'a - warto przeczytać manual, tzn. [eaton:2002], dostępny w dziale dokumentacji na stronach octave'a on-line pod adresem: http://www.gnu.org.

Kolejność zagadnień wykładanych w wykładzie z Matematyki Obliczeniowej może się różnić w kolejnych edycjach wykładu. W tym skrypcie przyjęliśmy kolejność zaprezentowaną w wykładzie z Matematyki Obliczeniowej z semestru letniego 2010/11.

W rozdziale  [link] omówimy krótko czym są w ogólności pakiety octave i matlab, oraz skąd można pobrać darmowe dystrybucje pakietu octave. W rozdziale Octave - podstawy omówimy podstawowe obiekty pakietu octave, czyli macierze i operacje na nich, oraz podstawowe struktury programistyczne: pętle, instrukcje warunkowe, operacje macierzowe i wektorowe, skrypty, operatory logiczne, itd. Rozdział Arytmetyka zmiennopozycyjna jest poświęcony arytmetyce zmiennopozycyjnej, czyli błędom zaokrągleń i ich wpływowi na obliczenia numeryczne. W kolejnym rozdziale Interpolacja wielomianowa i splajnowa omawiamy zagadnienie interpolacji wielomianowej i splajnowej. W rozdziale Obliczanie całek. Kwadratury mówimy o całkowaniu numerycznym. Rozdział Rozwiązywanie równań nieliniowych dotyczy metod rozwiązywania nieliniowych równań skalarnych. Rozdział Rozwiązywanie równań liniowych dotyczy numerycznej algebry liniowej, w szczególności rozwiązywania układów równań liniowych, a rozdział Liniowe zadania najmniejszych kwadratów dotyczy kolejnego zadania z numerycznej algebry liniowej, a mianowicie liniowego zadania najmniejszych kwadratów.

Octave i matlab - czym są te pakiety?

Pakiet octave jest zarówno środowiskiem obliczeń numerycznych, jak i językiem programowania.
Umożliwia on proste rozwiązywanie podstawowych zadań numerycznych takich jak: numeryczne obliczenie całki, rozwiązanie równań liniowych i nieliniowych, całkowanie przybliżone równań różniczkowych zwyczajnych itp. Można go używać z linii komend - jak kalkulator naukowy, pisać skrypty, czy m-pliki (czyli pliki z implementacjami funkcji octave'a).

Pakiet octave jest programem ogólnodostępnym jako wolne oprogramowanie.

Rozprowadzany jest na zasadach licencji GNU GPL.
Pakiet octave jest wolnym odpowiednikiem środowiska matlab, które jest programem komercyjnym, szeroko stosowanym do obliczeń numerycznych.

Zaletą octave'a jest to, że jest dostępny bez opłat. Octave jest dystrybuowany
w wersjach binarnych zarówno pod różne dystrybucje Linuxa, jak i w wersji binarnej pod system Windows.

Poniżej załączamy linki do stron octave'a:

  1. http://www.gnu.org/software/octave/index.html{Główna strona octave'a}
  2. http://www.gnu.org/software/octave/doc/interpreter/{Rozbudowany podręcznik do octave'a} - w języku angielskim dostępny on-line
  3. http://www.gnu.org/software/octave/download.html{Strona z linkami do plików z octavem}
  4. http://octave.sourceforge.net/{Strona octave-forge'a - czyli rozszerzeń octave'a}

Tutaj podajemy link to strony producenta matlaba:
http://www.mathworks.com/.
Jest to pakiet komercyjny.

Octave - podstawy

Pierwsza sesja w octave'ie

W systemie unix np. w linuxie najpierw musimy uruchomić terminal
(np. xterminal).
mo-Screenshot-term.jpg
Otworzy się wówczas okienko z linią komend z promptem - zazwyczaj migającym, w której możemy wpisywać komendy unixa, por. rysunek  [link].

Wpisujemy komendę octave i potwierdzamy Enter. W okienku terminala wyświetli się
nagłówek z numerem wersji i z linią komend octave'a, por. rysunek  [link].
mo-Screenshot-octave.jpg

Dalej, w linii komend octave'a wpisujemy polecenia octave'a.

Octave jako kalkulator, podstawowe zmienne

Najprostsze zastosowanie octave'a to kalkulator naukowy.
W linię komend wpisujemy np.: 234+76. Otrzymujemy:
ans=310 - zmiennej octave'a ans zostaje przypisana wartość $ 310 $, którą można dalej wykorzystać np. ans-1. Wówczas otrzymamy $ 309 $.

W octavie są funkcje takie, jak pierwiastek kwadratowy sqrt(), $ \sin() $, $ \cos() $, $ \exp() $ itp.

Jeśli chcemy zapamiętać wynik - musimy użyć zmiennej np.:

> a=12+34;
> sin(a)
ans =  0.90179

Typów zmiennych nie deklarujemy - octave domyślnie przyjmuje typ reprezentowany przez daną zmienną.

W octavie są np. zmienne typu zespolonego:

> a= sqrt(-1)
a =  0 + 1i
> b=1+4*i
b =  1 + 4i
> c=a*b
c = -4 + 1i
> d=a+b
d =  1 + 5i

Zachęcam do samodzielnego testowania.

Proszę zauważyć, że jeśli wywołamy jakąś funkcję albo operator, który zwraca wartość, np.: 2+3 zwraca wartość $ 5 $, i nie przypiszemy tej wartości zmiennej, to octave automatycznie przypisze tę zwracaną wartość do zmiennej ans, którą później można wykorzystać.

Wszystkie zmienne, jakie są w pamięci sesji octave'a, możemy wyświetlić poleceniem
who, czy whos - wersją poprzedniej funkcji, zwracającą
dokładny opis zmiennych.

Zmienne możemy usuwać z pamięci octave'a poleceniem clear np.

> a= sqrt(-1);
> b=1+4*i
b =  1 + 4i
> who
Variables in the current scope:

A    B    C    a    ans  b    c    d    f

> clear  a
> who
Variables in the current scope:

A    B    C    ans  b    c    d    f

Zmienna $ a $ zniknęła z sesji octave'a.
Podanie średnika po poleceniu spowoduje, że wynik nie zostanie wyświetlony w terminalu.

Pomoc do octave'a wywołujemy poleceniem help, pomoc
do konkretnej funkcji wywołujemy poleceniem help nazwa_funkcji np.
help sin:

> help sin
`sin' is a built-in function

 -- Mapping Function:  sin (X)
     Compute the sine for each element of X in radians.

     See also: asin, sind, sinh

Pomoc jest w języku angielskim.

Macierze i operacje na macierzach

Podstawowym typem zmiennym są macierze.

Macierz możemy zdefiniować wprost podając jej wartości:

A=[1,2,3,4;2 -4 6 7];

- proszę zauważyć, że elementy macierzy w wierszu oddzielamy przecinkiem albo spacją, a kolejne kolumny - średnikiem.

Do elementu macierzy o współrzędnych $ (i,j) $ odwołujemy się następująco:
A(i,j).

Jeśli zdefiniujemy element macierzy spoza zakresu, np. w tym przypadku

A(3,3)=7;

octave rozszerzy macierz, przyjmując domyślnie pozostałe niezdefiniowane wartości macierzy za zera.

Ważną rolę pełni operator średnik:

a:h:b

Generuje on wektor

$$[a, a+h,\ldots,a+k*h]$$

z $ a+k*h\leq b $.
Wywołanie:

a:b

działa jak a:1:b, czyli z $ h $ równym jeden.

Z macierzy możemy tworzyć inne macierze wycinając podmacierze.
Rozpatrzymy macierz $ A $ wymiaru $ M\times N $. Weźmiemy teraz dwa wektory indeksów
$ i=[i(1),\ldots,i(K)] $ oraz $ j=[j(1),\ldots,j(L)] $ dla $ 1\leq i(s)\leq M,1\leq j(s)\leq N $. Wtedy polecenie:

B=A(i,j);

stworzy macierz wymiaru $ K\times L $ taką, że
B(r,s)=A(i(r),j(s)), np.

> A=[3 3 3; 4 5 6]
A =

   3   3   3
   4   5   6

>  i=[ 2 1]; j=[2 3];
>  B=A(i,j)
B =

   5   6
   3   3

W ten sposób można wycinać też podmacierze lub np. kolumny, czy wiersze macierzy:

 A=[3 3 3; 4 5 6];
 B=A(1:2,2);

Warto zauważyć, że jeśli wywołamy polecenie x=A(:,3) to jest to równoważne x=A(1:2,3), czyli $ x $ staje się trzecią kolumną macierzy $ A $.

Macierze możemy tworzyć z innych macierzy blokowo:

B=[A,A];
C=[A;A];

Jeśli $ A $ miała wymiar $ M\times N $, to
macierz $ B $ ma wymiar $ M\times 2N $,
a jej odpowiednie podmacierze złożone z pierwszych $ N $ kolumn i ostatnich $ N $ kolumn są równe macierzy $ A $.
Macierz $ C $ analogicznie ma wymiar $ 2M\times N $, a jej podmacierze złożone z pierwszych i ostatnich $ M $ wierszy są równe $ A $. Przy takich definicjach musimy operować macierzami o
odpowiednich wymiarach.

Popatrzmy na wyniki:

octave:6> A=[3 3 3; 4 5 6]
A =

   3   3   3
   4   5   6

octave:7> B=A(1:2,2)
B =

   3
   5

octave:8> B=[A,A]
B =

   3   3   3   3   3   3
   4   5   6   4   5   6

octave:9> C=[A;A]
C =

   3   3   3
   4   5   6
   3   3   3
   4   5   6

octave:10>

Istnieją funkcje definiujących określone macierze - wymienimy teraz podstawowe:

  • zeros(M,N) - tworzy macierz zerową o $ M $ wierszach i $ N $ kolumnach.
  • ones(M,N) - tworzy macierz o wszystkich wartościach równych jeden o $ M $ wierszach i $ N $ kolumnach
  • eye(M,N) - tworzy macierz identycznościową o $ M $ wierszach i $ N $ kolumnach, tzn. o wartościach: jeden - na głównej przekątnej i zero - na pozostałych pozycjach
  • rand(M,N) - tworzy macierz o $ M $ wierszach i $ N $ kolumnach o losowych wartościach według rozkładu jednostajnego
    na odcinku $ [0,1] $
  • hilb(N) - tworzy macierz Hilberta $ N\times N $
  • vander(x) - tworzy macierz Vandermonde'a dla punktów $ \{x_k\} $ z wektora $ x=[x(1),\ldots,x(N)] $
  • diag() - tworzy macierz diagonalną, ale również pozwala otrzymać odpowiednie diagonale macierzy przy odpowiednim wywołaniu

Większość funkcji elementarnych typu $ \sin() $ może być wywoływana na macierzach. Wtedy:
w octave'ie a=sin(A) zwraca macierz tego samego wymiaru co $ A $,
ale o elementach $ (\sin(A(i,j))_{i,j} $. Mówimy, że implementacja takiej funkcji jest wektorowa. Popatrzmy na przykład:

octave:6> A=0.5*pi*[0:6;-3:3]'
A =

   0.00000  -4.71239
   1.57080  -3.14159
   3.14159  -1.57080
   4.71239   0.00000
   6.28319   1.57080
   7.85398   3.14159
   9.42478   4.71239

octave:7> a=sin(A)
a =

   0.00000   1.00000
   1.00000  -0.00000
   0.00000  -1.00000
  -1.00000   0.00000
  -0.00000   1.00000
   1.00000   0.00000
   0.00000  -1.00000

octave:8>

Macierze możemy dodawać i mnożyć, o ile zgadzają się wymiary:

octave:15> A=2+eye(2,2)
A =

   3   2
   2   3

octave:16> B=ones(2,2)+2*A
B =

   7   5
   5   7

octave:17> C=A*B
C =

   31   29
   29   31

octave:18>

Operacje arytmetyczne możemy też wywołać w wersji wektorowej, tzn.
jeśli np. dwie macierze $ A,B $ mają ten sam wymiar, to wywołanie
C=A.*B zwróci macierz $ C $ tego samego wymiaru taką, że
$ C(k,l)=A(k,l)*B(k,l) $. To samo dotyczy innych działań:

  • dzielenia C=A.\B
  • odejmowania C=A.-B
  • dodawania C=A.+B
  • podnoszenia do potęgi C=A.^B

W przypadku dodawania i odejmowania operatory bez kropek działają tak samo.

Macierz rzeczywistą możemy transponować operatorem C=A'.
W przypadku macierzy zespolonej operator ten zwróci sprzężenie hermitowskie, a operator z kropką A.' zwróci zwykłe transponowanie. Najlepiej zobaczyć to na przykładzie:

octave:37> A=[1+i -3+4i;-8+2i 5+7i]
A =

   1 + 1i  -3 + 4i
  -8 + 2i   5 + 7i

octave:38> A'
ans =

   1 - 1i  -8 - 2i
  -3 - 4i   5 - 7i

octave:39> A.'
ans =

   1 + 1i  -8 + 2i
  -3 + 4i   5 + 7i

octave:40>

Istnieje wiele funkcji pozwalających na manipulowanie macierzami np.:

  • fliplr() - odwraca kolejność kolumn
  • flipud() - odwraca kolejność wierszy
  • reshape() - tworzy nową macierz o określonych wymiarach z elementów macierzy wyjściowej
  • rot90() - obraca macierz o wielokrotność $ 90 $ stopni w kierunku przeciwnym obrotowi wskazówek zegara

Wektory traktowane są jako macierze o jednej kolumnie, albo jednym wierszu, jakkolwiek istnieją specjalne funkcje działające tylko na wektorach, np.
funkcja length(x) zwraca długość wektora $ x $ pionowego lub poziomego.

Polecenie [m,n]=size(A) zwraca wymiar macierzy jako dwie wartości w wektorze.

octave:40> x=[1 2 3 4]
x =

   1   2   3   4

octave:41> length(x)
ans =  4
octave:42> size(x)
ans =

   1   4

octave:43>

Octave jako język programowania - skryptu i m-pliki

Octave może być również traktowany jako język programowania.

Istnieje możliwość pisania skryptów, a następnie ich wywoływania
z linii komend octave'a.

Skrypty to po prostu pliki tekstowe o rozszerzeniu \textsl{.m}, tzn. np. skryptem może być plik o nazwie \textsl{nazwaskryptu.m}, zawierające w kolejnych liniach komendy octave'a
utworzone przez dowolny edytor tekstowy, np. vi w unixie, czy \lstinline{wordpad} w systemie windows.

Skrypt wywołujemy podając nazwę bez rozszerzenia, czyli skrypt zawarty w pliku tekstowym \textsl{skrypt.m} wywołujemy z linii komend komendą:

skrypt

Po wywołaniu np. następującego skryptu:

#skrypt.m
#Prosty skrypt - ta linia to komentarz
 a=2+3

octave wykona wszystkie polecenia po kolei - o ile nie będzie w nich jakiegoś błędu. W tym przypadku zwróci a=5. Tak więc octave jest interpreterem - nie ma etapu kompilacji skryptu przed wykonaniem.

Pisząc skrypty w octavie najlepiej mieć na pulpicie komputera co najmniej dwa okna - jedno z sesją octave'a, a drugie z edytorem, w którym
edytujemy plik tekstowy ze skryptem.

Tu warto dodać, że komentarze w octavie wstawiamy po znakach # lub %.

Octave jako język programowania dostarcza również większość konstrukcji imperatywnych, m.in.

  • instrukcja warunkowa if - przykład:
    #instrukcja warunkowa
    if(x>0)
     x=1;
    else
     x=0;
    endif;
  • instrukcja switch wielokrotnego wyboru zastępująca
    if jeśli $ k $ przyjmuje określoną ilość wartości

    switch (taknie)
       case {"T" "t"}
            x = 1;
        case {"N" "n" }
             x = 0;
         otherwise
             error ("tylko t lub n");
    endswitch
  • pętla while( )
    #sumowanie
    k=0;
    a=0;
    while(k<=10)
    do
     a+=k;
     k++;
    endwhile;!
  • pętla z warunkiem zatrzymania na końcu:
    do  intrukcje  until( );

    k=0;
    do
     k++;
    until(k>10)
  • pętla for:
    a=0;
    for k=1:10,
     a+=k;
    endfor
  • funkcje octave'a function []=f()
    function [s]=f(a,b)
     s=a+b;
    endfunction

    Więcej informacji o funkcjach podamy w następnym rozdziale.

W tym miejscu omówimy operacje logiczne octave'a.
Ogólnie ujmując - wartość różna od zera traktowana jest jako prawda, $ 0 $ to fałsz.
Operatory logiczne octave'a są takie same jak w języku C tzn:

  •  a && b - operator logiczny dwuargumentowy $ a \;I \;b $ zwraca wartość jeden (prawda), jak
    $ a $ i $ b $ są różne od zera, w przeciwnym przypadku operator zwraca zero (fałsz)
  • a||b - operator logiczny dwuargumentowy $ a\; LUB\; b $ zwraca jeden, jeśli któreś z
    $ a $ lub $ b $ jest równe jeden, w przeciwnym przypadku zwraca zero
  • \lstinline+!a+ - operator logiczny jednoargumentowy negacji $  NIEPRAWDA \; a $, jeśli $ a $ jest różne od zera - to operator zwraca zero, w przeciwnym przypadku operator zwraca jeden
  • a==b - operator dwuargumentowy równości, jeśli $ a $ ma tę samą wartość co $ b $, to operator zwraca jeden, w przeciwnym razie operator zwraca zero
  • \lstinline+a!=b+ - operator dwuargumentowy bycia różnym, jeśli $ a $ jest różne od $ b $ to zwraca jeden, w przeciwnym razie operator zwraca zero
  •  a < b - operator dwuargumentowy porównania, zwraca jeden gdy $ a $ jest mniejsze od $ b $, w przeciwnym razie operator zwraca zero. Pozostałe operatory porównania ( a > b,  a<=b,  a>= b) są zdefiniowane analogicznie

Funkcje

(#)

Funkcje w octave'ie możemy definiować z linii komend, czy jako część kodu w skrypcie:

function [w,s]=f(a,b)
#pomoc do danej funkcji
 s=a+b;
 w=a-b;
endfunction

Wywołujemy je z linii komend (lub w linii skryptu, czy w innej funkcji)

w=f(1,2)
[w,s]=f(1,2)

Pierwsze wywołanie spowoduje, że funkcja zwróci tylko pierwszą wartość tj. w.

Niektóre ostatnie parametry funkcji mogą przyjmować domyślną wartość tzn.

function [w,s]=f(a,b=0)
#pomoc do danej funkcji
 s=a+b;
 w=a-b;
endfunction

Wtedy wywołanie funkcji

[w,s]=f(1)

zwróci wartości funkcji dla parametru $ b=0 $.

Funkcje można przekazywać jako parametry do innych funkcji poprzez przekazanie nazwy funkcji jako napisu - ciągu znaków lub jako tzw. uchwyty do funkcji - tu przykład przekazania nazwy funkcji $ sin() $ do funkcji quad() obliczającej całki:

> f=@sin;
> quad("sin",0,1)
ans =  0.45970
> quad(f,0,1)
ans =  0.45970
> quad(@sin,0,1)
ans =  0.45970

Oba sposoby są równoważne.

Proste funkcje można definiować jako tzw. funkcje anonimowe:

> f=@(x) + sin(x); #f(x)=x+sin(x)
> f(1)
ans =  0.84147

czyli definiujemy od razu uchwyt do funkcji.

W octavie wprowadzono za matlabem tzw. m-pliki, czy - inaczej - pliki funkcyjne. Są to pliki z rozszerzeniem \textsl{m} np. \textsl{f.m}, zawierające definicje funkcji o tej samej nazwie co plik (bez rozszerzenia).
M-plik o nazwie \textsl{test.m} powinien więc zawierać definicję
funkcji test(), czyli w pierwszej linii powinien znajdować się
nagłówek funkcji.

Poniżej widzimy ogólną strukturę funkcji:

function [zwracane wartosci]=test(parametry)
##pomoc do funkcji
#   komendy

   zwracane wartosci powinny byc zdefiniowane

endfunction

Wywołanie polecenia help test powinno wyświetlić pomoc do funkcji tzn. wykomentowane linie spod nagłówka funkcji.

Funkcje można też definiować w skryptach, które mają te same rozszerzenie \textsl{m}.

Ważne, aby w pierwszej niepustej linii skryptu nie było słowa kluczowego function, tzn. aby przed definiowaniem funkcji w skrycie
podać jakąkolwiek komendę. W przeciwnym przypadku octave potraktuje skrypt jako m-plik.

W m-plikach, pod definicją głównej funkcji o tej samej nazwie co m-plik, mogą się znajdować tzw. pod-funkcje
(ang. subfunctions), które mogą być wywołane przez tę główną funkcję z m-pliku i inne pod-funkcje tylko w tym m-pliku.
Tutaj podajemy przykład takiego m-pliku z funkcją główną obliczającą przybliżenie normy typu $ L^2 $ i pod-funkcjami:

#m-plik czyli plik tekstowy normaL2.m

function [nL2]=normaL2(FCN,a,b)
#norma L^2(a,b) z funkcji o uchwycie FCN
# na [a,b]
   sf=uchwytkwad(FCN);
   nL2=sqrt(quad(sf,a,b));
endfunction

function [SFCN]=uchwytkwad(f)
#tworzy uchwyt do funkcji g(x)=f(x)*f(x)
 SFCN=@(x) f(x)*f(x);
endfunction

#koniec m-pliku normaL2.m

Oczywiście pod-funkcja jest tutaj wprowadzona sztucznie, nie jest ona konieczna.

W funkcji octave'a można używać automatycznych zmiennych nargin i nargout, które
przy wywołaniu danej funkcji przyjmują wartości:

  • nargin -
    ilości parametrów, z jaką funkcja została wywołana.
  • nargout -
    ilości zwracanych wartości, z jaką funkcja została wywołana.

Tu podajemy przykład wykorzystania tych zmiennych w konkretnej funkcji:

function [y1,y2]=test(a,b)
# test nargin nargout
  if(nargin<2)
   b=a;
  endif
  y1=a+b;
  if(nargout>1)
    y2=a-b;
  endif
endfunction

Sprawdźmy, jak taka funkcja działa:

octave:46> test(1)
ans =  2
octave:47> test(1,2)
ans =  3
octave:48> [y1]=test(1,2)
y1 =  3
octave:49> [y1,y2]=test(1,2)
y1 =  3
y2 = -1
octave:50>

Podstawowa grafika

Octave posługuje się zazwyczaj zewnętrznym narzędziem do grafiki w linuxie.
Najczęściej jest to program gnuplot, ewentualnie grace lub inny.
W przypadku systemu windows w binarnej wersji octave'a zawarta jest też grafika.
mo-wyksin.jpg
Podstawową funkcją umożliwiającą rysowanie wykresów jest plot().

Najprostsze wywołanie: plot(x) dla wektora $ x $ otworzy nowe okno i w nim
narysuje punkty $ (k,x(k)) $. Punkty mogą pozostać połączone prostymi w zależności od ustawienia domyślnych opcji octave'a.

Innym wywołaniem plot() pozwalającym np. na rysowanie wykresów jest plot(x,y) dla wektorów $ x,y $ tej samej długości. W tym przypadku octave narysuje punkty $ (x(k),y(k)) $.

Tak więc, jeśli stworzymy w wektorze $ x $ siatkę stupunktową równomierną na $ [0,3] $ - wygodnym poleceniem jest linspace(,a,b,N) - a potem policzymy wartości np. funkcji $ y(k)=\sin(x(k)) $, to polecenie
plot(x,y) narysuje wykres $ \sin() $ na tym odcinku, por. rysunek  [link].

x=linspace(0,3);
y=sin(x); #sin jest wektorowa funkcja
plot(x,y)

mo-wyksc.jpg
Istnieje możliwość umieszczania kilku wykresów na jednym obrazku np.

x=linspace(0,3);
y=sin(x); #sin jest wektorowa funkcja
z=cos(x);
plot(x,y,x,z)

Tu plot(x,y,x,z) narysuje punkty $ (x(k),y(k) $ i $ (x(k),z(k) $, czyli wykresy $ \sin() $ i $ \cos() $ na jednym rysunku.

mo-wykstyle.jpg
Inną możliwością jest zastosowanie odpowiednich funkcji dotyczących grafiki; np.
hold on powoduje, że kolejne wywołania plot sprawia rysowanie wykresów na jednym obrazku.
Polecenie hold off likwiduje tę własność.

Możemy dobrać kolor wykresu, styl wykresu, dodać opis, np. wykres sinusa może być w kolorze czerwonym narysowany plusami, a cosinusa w kolorze niebieskim narysowany gwiazdkami

x=linspace(0,3,40);
y=sin(x); #sin jest wektorowa funkcja
z=cos(x);
plot(x,y,"+r;sin;",x,z,"*b;cos;")
#+ - wykres plusami r - red (czerwony)
#* wykres gwiazdkami, b -blue (niebieski)

Podamy teraz kilka przydatnych funkcji związanych z grafiką:

  • xlabel - etykieta osi poziomej
  • ylabel - etykieta osi pionowej
  • title - tytuł rysunku
  • figure - kontroluje okienka z
    wykresami - możemy ich mieć kilka, czyli

    x=linspace(0,3,40);
    figure(1);
    y=sin(x);
    plot(x,y)
    figure(2)
    z=cos(x);
    plot(x,z)

    stworzy dwa okienka z wykresem $ \sin() $ w pierwszym i $ \cos() $ w drugim.

  • semilogx - wykres w skali półlogarytmicznej - względem zmiennej poziomej
  • semilogy - wykres w skali półlogarytmicznej - względem zmiennej pionowej
  • loglog - wykres w skali logarytmicznej
  • bar - wykres słupkowy.
  • stairs - wykres schodkowy.

Zachęcam do zapoznania się z tymi funkcjami używając pomocy octave'a.

Istnieje też możliwość używania prostej grafiki trójwymiarowej, ale
w przypadku matematyki obliczeniowej nie będziemy tego potrzebowali.

Przykładowy skrypt

Poniżej zamieszczamy prosty skrypt, w którym
przedstawiamy przykłady użycia większości operacji, zmiennych omawianych w tym rozdziale.

\lstinputlisting{./mobasic.m}

Przykładowa sesja matlaba

W tym rozdziale przedstawimy przykładową sesję matlaba.
mo-matlab.jpg

mo-matlab1.jpg

W przeciwieństwie do octave'a, matlab posiada środowisko graficzne.
Ogólnie matlaba wywołujemy klikając na odpowiednią ikonę, czy wybierając program matlab w odpowiednim menu systemu.
W linuxie możemy też wywołać ten program z linii komend komendą:
\textsl{matlab}. Wówczas powinno otworzyć się środowisko matlaba, por. rysunek  [link].

Jak widać, w środku mieści się \textsl{Command Window}, czyli okno komend, w którym możemy wpisywać komendy matlaba, z lewej strony widać podgląd bieżącego katalogu pod nazwą \textsl{Current Folder}, z prawej strony, z góry pod nazwą \textsl{Workspace}, widać okno ze spisem zdefiniowanych zmiennych oraz z dołu, z lewej strony - małe okno historii komend, por. rysunek  [link].

Ustawienia można zmieniać, np. można część okien ukryć, czy odłączyć.
Po szczegóły odsyłam do pomocy do matlaba.

Wpiszmy kilka komend takich samych jak w octave'ie w okno komend:

A=[1:3;-2 4 7]
B=[A; A]

W oknie komend widać, że wynik jest taki sam jak w octave'ie. Natomiast w oknie
\textsl{Workspace} widzimy, że w sesji są zdefiniowane dwie zmienne A,B, a w oknie historii widzimy dwie nasze komendy, por rysunek  [link].

W środowisku matlaba najważniejsze jest okno komend, które pełni tę samą rolę co terminal z wywołanym octavem. Pozostałe okna pełnią rolę pomocniczą ale
niewątpliwie
środowisko matlaba ułatwia pracę.
Dalej nie będziemy się zajmować szczegółowo matlabem, ale zdecydowana większość kodów z tego skryptu działa zarówno w octave'ie, jak i w matlabie.

Arytmetyka zmiennopozycyjna

Wszystkie obliczenia w octavie są wykonywane w arytmetyce zmiennopozycyjnej (inaczej - arytmetyce fl) podwójnej precyzji (double) - choć w najnowszych wersjach octave'a istnieje możliwość używania zmiennych typu single, czyli zmiennych w pojedynczej precyzji arytmetyki zmiennopozycyjnej.

Dokładność względna arytmetyki podwójnej precyzji wynosi ok. $ 10^{-16} $, a dokładność arytmetyki pojedynczej precyzji wynosi ok. $ 10^{-7} $.

W octavie jest kilka funkcji związanych bezpośrednio z własnościami arytmetyki fl, mianowicie:

  • eps - zwraca epsilon maszynowy arytmetyki, tzn. najmniejszą liczbę $ \epsilon $ taką, że $ fl(\epsilon+1)>1 $. Tu $ fl() $ oznacza wynik działania w arytmetyce zmiennopozycyjnej. Octave domyślnie działa na liczbach zmiennopozycyjnych w podwójnej precyzji arytmetyki więc oczywiście eps zwróci nam epsilon maszynowy dla arytmetyki podwójnej precyzji. Z kolei wywołanie eps(a) - dla $ a $ liczby zwróci odległość od $ a $ do najbliższej większej od $ a $ liczby w arytmetyce zmiennopozycyjnej
  • single(x) - funkcja konwertująca zmienną np. typu double, czyli w podwójnej precyzji arytmetyki, do zmiennej typu single, czyli w pojedynczej precyzji arytmetyki zmiennopozycyjnej
  • double(x) - funkcja konwertująca zmienną do typu zmiennej w podwójne precyzji arytmetyki, czyli odwrotna do funkcji single(x). Podwójna precyzja jest w octavie domyślna, więc wydaje się że, głównym zastosowaniem tej funkcji jest konwersja zmiennej typu single z powrotem do typu double

Zastosujmy funkcję eps - sprawdźmy, czy rzeczywiście $ eps+1 $ w fl jest większe od jeden, a np. $ eps/2+1 $ już nie - tu fragment kodu octave'a:

 if((eps+1)>1)
  printf("%g+1>1 w fl\n",eps);
 else
  printf("%g nie jest epsilonem maszynowym\n",eps);
 endif

 if((1+eps/2)==1)
  printf("fl(%g+1)=1\n",eps/2);
 endif

Dwa sposoby obliczania <img class= w arytmetyce zmiennopozycyjnej" />
Przy pomocy eps możemy znaleźć najmniejszą liczbę dodatnią w arytmetyce podwójnej precyzji:

octave:33> x=eps(0)
x =  4.9407e-324
octave:34>

Sprawdźmy:

octave:34> x/2>0
ans = 0
octave:35>

Sprawdźmy, jak działają funkcje związane ze zmianą precyzji arytmetyki single
i double - tu fragment sesji octave'a:

 b=single(eps)
 if((single(eps)+1)==1)
 printf("fl(single(eps)+1)==1\n");
 endif

 if((double(b)+1)>1)
  printf("fl(double(single(eps))+1)>1.\n");
 endif
 #a teraz epsilon maszynowy precyzja pojedyncza
 eps(single(1))

Otrzymaliśmy, że epsilon maszynowy w podwójnej precyzji wynosi:
$ 2.2204e-16 $, a w pojedynczej precyzji arytmetyki:
$ 1.1921e-07  $.

Teraz znajdziemy najmniejszą liczbę dodatnią w arytmetyce pojedynczej precyzji:

octave:52> eps(single(0))
ans =  1.4013e-45
octave:53>

Redukcja cyfr przy odejmowaniu

Inną własnością arytmetyki fl, która może powodować problemy, jest tzw. redukcja cyfr przy odejmowaniu liczb tego samego znaku.

Przykład

Przetestujmy to na dwóch równoważnych wzorach na
funkcję

$$f(x)=\cos(x)-1=1-2*\sin^2(0.5*x)-1=-2*\sin^2(0.5*x).$$

Przetestujmy oba wzory dla $ x $ bliskich zeru:

f1=@(x) cos(x)-1;
f2=@(x) -2*sin(0.5*x).*sin(0.5*x);
a=1e-7;
x=linspace(-a,a,1000);
plot(x,f1(x),";wzor 1;",x,f2(x),";wzor 2;");

Wykres na rysunku  [link] pokazuje, że wzór drugi jest lepszy ze względu na arytmetykę fl.

Możemy oczywiście policzyć też błąd pomiędzy wynikami:

 er=f1(1e-7)-f2(1e-7)

Różnica $ 3.9964e-18 $ jest
pozornie mała, bo rzędu $ 10^{-18} $. Jednak
$ f2(10^{-7}) $ jest rzędu $ 10^{-14} $, więc
dokładność względna jest rzędu $ 10^{-4} $.
Jest to bardzo mała dokładność, jak na arytmetykę podwójnej precyzji, w której błąd względny obliczeń jest na poziomie $ 10^{-16} $.

Możemy policzyć względny błąd względem $ f2(10^{-7}) $:

 abs((f1(1e-7)-f2(1e-7))/f2(1e-7))

i otrzymamy $ 7.9928e-04 $.

Obliczanie funkcji z rozwinięcia w szereg

Inny przykład na problemy redukcją cyfr, to obliczanie wartości funkcji
korzystające z rozwinięcia w szereg o wyrazach o różnych znakach.

Rozpatrzmy funkcję $ \sin(x) $, która ma następujące rozwinięcie w globalnie zbieżny
szereg:

$$<br />
 \sin(x)=x - \frac{x^3}{3!}+\frac{x^5}{5!} - \frac{x^7}{7!}+\ldots.<br />
$$

Wykres funkcji sinus  obliczanej z rozwinięcia w szereg na <img class=" />
Zaimplementujmy w octavie funkcję obliczającą przybliżenie funkcji matematycznej $ \sin(x) $, korzystające z odpowiedniego obcięcia tego rozwinięcia dla dowolnego $ x $ rzeczywistego. Implementacja tej funkcji będzie wektorowa, tzn. wywołanie jej z macierzą $ X $ jako argumentem zwróci macierz przybliżeń wartości $ \sin $ na odpowiednich elementach macierzy $ X $.

function y=naszsin(x,N=300)
#obliczamy wartosc sin(x)
#korzystajac z rozwiniecia w szereg :
#x-x^3/3!+x^5/5!-....(-1)^(k+1)*x^{2k+1}/(2k+1)!+...
#N - ilosc wyrazow w szeregu - domyslnie 300
 kx=-x.*x;
 y=x;
 for k=1:N,
   x=x.*kx./((2*k)*(2*k+1));
   y+=x;
 endfor
endfunction

Przetestujmy naszą funkcję na $ [-\pi,\pi] $.

w=linspace(-pi,pi);
plot(w,sin(w),";sin;",w,naszsin(w),";naszsin;")

Na wykresie nie widać różnicy.

Wykres błędu w skali półlogarytmicznej obliczania funkcji sinus z rozwinięcia w szereg
Przetestujmy naszą funkcję dla większych argumentów: $ [11*\pi,13*\pi] $:

w=linspace(11*pi,13*pi);
plot(w,sin(w),";sin;",w,naszsin(w),";naszsin;")

Widać z rysunku  [link], że rozwijanie w szereg nie zawsze ma sens z powodu
własności arytmetyki zmiennopozycyjnej.

Można postawić pytanie, czy w rozwinięciu nie było za mało wyrazów szeregu?

w=linspace(11*pi,13*pi);
plot(w,sin(w),";sin;",w,naszsin(w,100000),";naszsin;")

Proszę sprawdzić, że zwiększenie ilości wyrazów szeregu nic nie zmieniło.

Popatrzmy, jak rośnie błąd na skali logarytmicznej, por. rysunek  [link]:

w=linspace(eps,13*pi,1000);
semilogy(w,abs(sin(w)-naszsin(w)))

Błąd w skali logarytmicznej rośnie liniowo. To oznacza, że
błąd rośnie wykładniczo.
Wykres błędu przy obliczaniu pochodnej ilorazem różnicowym <img class= dla $ x=1 $ w zależności od $ h $." />

Obliczanie przybliżonej wartości pochodnej

Kolejnym przykładem sytuacji, w której mogą pojawić się problemy, to obliczanie pochodnej za pomocą ilorazów różnicowych, czyli wprost z definicji.
Wykres błędu przy obliczaniu pochodnej ilorazem różnicowym <img class= dla $ x=1 $ w zależności od $ h $ w skali półlogarytmicznej" />
Przetestujmy przybliżone obliczanie pochodnej z definicji, czyli ze wzoru:

$$<br />
  f'(x)\approx \frac{f(x+h)-f(x)}{h}.<br />
$$

Wydawałoby się, że im $ h $ jest mniejsze, tym lepsze przybliżenie otrzymamy.
Ale np. dla $ f(x)=\sin(x) $ z $ x=1 $ i $ h<eps $, np. dla $ h=eps/2 $, otrzymamy w arytmetyce fl:

$$<br />
 fl(x+h)=fl(1+h)=1.<br />
$$

Sprawdźmy:

> h=eps/2
h =  1.1102e-16
> sin(1+h)-sin(1)
ans = 0
#czy na pewno
>sin(1+h)==sin(1)
ans=1

Czyli rzeczywiście $ h $ nie może być dowolnie małe.

Spróbujmy znaleźć dla tego przykładu optymalne $ h $ postaci $ h=2^{-n} $ eksperymentalnie.

f=@sin;
h=1;x=1;
fx=f(x); dfx=cos(x);
er=zeros(53,1);
eropt=er(1)=abs((f(x+h)-fx)/h-dfx);
hopt=h;n=nopt=0;
while(h>=eps)
  h/=2;
  n++;
  er(n)=abs((f(x+h)-fx)/h - dfx);
 if(er(n)<eropt)
    eropt=er(n);
    nopt=n;
    hopt=h;
 endif
endwhile
printf("optymalne h=%g=2^(-%d)\n",hopt,nopt);
plot(er);

Na wykresie nie widać dokładnie, por. rysunek  [link], gdzie jest minimum. Na wykresie w skali pół-logarytmicznej minimum jest widoczne, por. rysunek  [link]:

 semilogy(er)

Widzimy, że optymalne $ h $ w tym przypadku jest rzędu $ 2^{-27} $, czyli ok $ 10^{-9} $.

Można powtórzyć te same obliczenia dla przybliżania pochodnej za pomocą tzw. ilorazu centralnego:

$$<br />
  f'(x)\approx \frac{f(x+h)-f(x-h)}{2h}.<br />
$$

W kolejnych rozdziałach będziemy omawiali wpływ
arytmetyki zmiennopozycyjnej na wyniki obliczeń dla konkretnych zadań.

Interpolacja wielomianowa i splajnowa

Wielomiany w octavie

W octavie istnieje cała gamma funkcji związanych z wielomianami:

  • polyval( ) - funkcja pozwalająca obliczać wartość wielomianu
    zadanego w bazie potęgowej dla danej wartości $ x $, czy całej macierzy wartości
  • polyfit( ) - funkcja znajdująca współczynniki wielomianu zadanego stopnia najlepiej dopasowanego do zadanej tabelki punktów.
  • polyinteg( ) - zwraca współczynniki wielomianu będącego całką nieoznaczoną z danego wielomianu
  • polyderiv( ) - zwraca współczynniki wielomianu będącego pochodną z danego wielomianu
  • roots( ) - zwraca zera wielomianu
  • conv(a,b) - zwraca współczynniki wielomianu będącego iloczynem wielomianów o współczynnikach odpowiednio z wektorów $ a,b $

Będziemy korzystali przede wszystkim z dwóch pierwszych funkcji.

We wszystkich tych funkcjach przyjmowane jest, że rozpatrujemy współczynniki wielomianu w bazie potęgowej, ale w następującej kolejności:

$$<br />
  (x^n,x^{n-1},\ldots,1).<br />
$$

Współczynniki indeksowane są od jedynki, tzn. wielomian stopnia nie większego od $ n $:

$$w(x)=\sum_{k=1}^{n+1} a_kx^{n+1-k}$$

jest reprezentowany przez wektor współczynników: $ (a_1,\ldots,a_{n+1}) $.

Jeśli chcemy policzyć wartość $ w $ w danym punkcie $ x $, czy ewentualnie dla tablicy punktów umieszczonych w macierzy $ X $,
wywołujemy polyval(a,x) lub polyval(a,X).

Tak więc, aby narysować wykres funkcji $ x^3-2x+1 $ na $ [-3,4] $ można wykonać następującą sekwencję komend :

 a=[1,0,-2,1];  #definiujemy wsp.
 x=linspace(-3,4);#siatka
 wx=polyval(a,x); #wartosc w na siatce
 plot(x,wx); #wykres

Interpolacja wielomianowa

Interpolacja wielomianowa polega na tym, że szukamy funkcji $ p $ z pewnej przestrzeni wielomianów (zazwyczaj wielomianów stopnia nie większego od $ n $), które w tych punktach spełniają odpowiednie warunki interpolacyjne, tzn. $ p $ oraz jej pochodne mają w tych punktach zadane wartości.

Interpolacja Lagrange'a

mo-ILsimplest.jpg
W interpolacji wielomianowej Lagrange'a dla zadanych $ n+1 $ różnych punktów $ x_k $ i wartości $ y_k=f(x_k) $ szukamy wielomianu $ p(x) $ stopnia co najwyżej $ n $ takiego, że spełnione są następujące warunki interpolacyjne:

$$<br />
 f(x_k)=p(x_k)=y_k \qquad k=1,\ldots,n+1.<br />
$$

Funkcja polyfit(x,y,n) znajduje
współczynniki wielomianu $ p $ w bazie potęgowej dla wektorów $ x,y $ długości $ n+1 $. Ważne aby wartości w $ x $ były różne.

Przetestujmy tę funkcję dla następujących danych:
$ x=(-1,0,2,4) $, $ y=(1,1,2,-1) $:

 x=[-1,0,2,4];
 y=[1,1,2,-1];
 a=polyfit(x,y,3);
 s=linspace(-2,5);
 plot(s,polyval(a,s),x,y,"r+;war interp;")

Widać na rysunku  [link], że wielomian spełnia warunki interpolacyjne.
Możemy to sprawdzić też obliczeniowo:

max(abs(y-polyval(a,x)))

Otrzymaliśmy wynik: błąd w węzłach jest równy prawie zero.

Przetestujmy polyfit(x,y,n) dla dużych $ n $ dla węzłów równoodległych na $ [0,1] $ i losowe wartości $ y $ z $ [-1,1] $, następnie policzmy błąd w węzłach w zależności od $ N $:

 N=100;
 er=zeros(N+1,1);
 for k= 0:N,
   x=linspace(0,1,k+1);
   y=2*(rand(size(x))-0.5);
   a=polyfit(x,y,k);
   er(k+1)=max(abs(y-polyval(a,x)));
 endfor
 max(er)

Widzimy, że możemy otrzymać niepoprawny wynik dla dużych $ N $.
Wynika to ze złych własności algorytmu stosowanego przez octave'a ze względu na zaburzenia powodowane przez niedokładne obliczenia w arytmetyce zmiennopozycyjnej.

Zbadajmy błąd pomiędzy funkcją, a jej wielomianem interpolacyjnym.
Na początek rozpatrzmy analityczną funkcję $ f(x)=\sin(x) $ i węzły równoodległe, oraz normę supremum na odcinku $ [-4,6] $. Normę supremum:

$$\|g\|_{\infty,[a,b]}:=\sum_{t\in[a,b]}|g(t)|$$

przybliżymy poprzez dyskretną normę na siatce $ N $ równomiernie rozłożonych punktów na tym odcinku $ [a,b] $, czyli

$$<br />
\|f\|_{h,\infty,[a,b]}=\max_k |f(a+k*h)|<br />
$$

z $ h=(b-a)/N $:

 x=[-1,0,2,4];
 y=[1,1,2,-1];
 a=polyfit(x,y,3);
 s=linspace(-2,5);
plot(s,polyval(a,s),x,y,"r+;war interp;")

mo-IHsin3x.jpg

Interpolacja Hermite'a

Rozpatrzmy teraz zadanie interpolacji Hermite'a.
W interpolacji wielomianowej Hermite'a dla zadanych $ n+1 $ różnych punktów $ x_k $ i naturalnych krotności węzłów $ p_k $ szukamy wielomianu $ w(x) $ stopnia co najwyżej $ N $ dla $ N=\sum_{k=0}^np_k-1 $ takiego, że spełnione są następujące warunki interpolacyjne:

$$<br />
 w^{(j)}(x_k)=y_{k,j} \qquad k=1,\ldots,n+1, \qquad j=0,\ldots,p_k-1.<br />
$$

Tutaj $ y_{k,j} $ to $ N+1 $ zadanych wartości.

W octavie nie ma funkcji realizującej interpolacje Hermite'a.

Ale czy w szczególnych przypadkach nie możemy łatwo rozwiązać tego problemu korzystając z odpowiednich funkcji octave'a?

Rozpatrzmy przypadek węzłów tej samej krotności; np. $ n $ różnych dwukrotnych węzłów.

Chcemy znaleźć współczynniki wielomianu $ \sum_{k=1}^N a_kx^k $ dla $ N=2(n+1) $ takie, że

$$<br />
w(x_k)=f(x_k), \qquad w'(x_k)=f'(x_k).<br />
$$

Czyli rozwiązanie spełnia następujący układ równań liniowych:
\begin{eqnarray*}
\sum_{k=0}^N a_kx_j^k=f(x_k) && j=0,\ldots,n \\
\sum_{k=1}^N ka_kx_j^{k-1}=f'(x_k) && j=0,\ldots,n.
\end{eqnarray*}

Utwórzmy macierz tego układu z pomocą funkcji vander(x), która tworzy macierz Vandermonde'a dla węzłów podanych w wektorze $ x $:

 function [a,C,f]=interpolyH(x,y,dy)
#tworzy wspolczynniki wielomianu Hermite'a dla
#wezlow dwukrotnych z x
#y= wartosci w wezlach x (wektor pionowy)
#wartosci pochodnej w wezlach x (wektor pionowy)
#output: wspolczynniki wielomainu w bazie potegowej
# (Zgodne z polyval())
#C- opcja - macierz ukladu
n=length(x);
N=2*n-1;
A=vander(x,N+1);
D=diag(N:-1:0);
B=zeros(size(A));
B(:,1:N)=A(:,2:N+1);
B=B*D;
C=[A;B];
f=[y;dy];
a=C\[y;dy];
endfunction

mo-IH3w.jpg
Na początku przetestujmy tę funkcję
dla węzłów $ -1,1 $ dwukrotnych i funkcji $ \sin(3*x) $.

x=[-1;1];
f=@(x) sin(3*x);
df=@(x) 3*cos(3*x);
y=f(x);
dy=df(x);
a=interpolyH(x,y,dy);
z=linspace(-1.1,1.1);
plot(z,polyval(a,z),";wiel H;",z,f(z),";f;",x,y,"r+");

Z wykresu funkcji i wielomianu widzimy, że
wielomian przecina wykres i jest styczny w punktach $ -1,1 $, por. rysunek  [link].

Zobaczmy co się stanie w przypadku trzech węzłów.
Testujemy funkcję dla węzłów $ -2,0,2 $ dwukrotnych i $ \sin(3*x) $.

x=[-1;0;1];
f=@(x) sin(3*x);
df=@(x) 3*cos(3*x);
y=f(x);
dy=df(x);
a=interpolyH(x,y,dy);
z=linspace(-1.2,1.2);
plot(z,polyval(a,z),";wiel H;",z,f(z),";f;",x,y,"r+");

Z rysunku  [link] widać, że funkcja działa poprawnie również w tym przypadku.
%%
% SECTION SPLAJNY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

Interpolacja zespolona. Algorytm FFT

(#)
%%
% SECTION FFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W tym rozdziale omówimy krótko interpolację zespoloną, ale tylko w przypadku
określonych węzłów zespolonych równomiernie rozmieszczonych na sferze jednostkowej.

Chcemy znaleźć współczynniki zespolone $ a_k\in \mathbb{C} $, $ k=0,\ldots,N $ takie, że

$$<br />
   \sum_{k=0}^N a_k z_j^k=b_k \quad j=0,\ldots,N<br />
$$

dla danych zespolonych $ b_k $ i $ z_j=\exp\left(\frac{2\pi*i*j}{N+1}\right) $ dla $ j=0,\ldots,N $, czyli pierwiastków z jedynki stopnia $ N+1 $.

Równoważnie możemy to zadanie sformułować jako zadanie znalezienia $ a_k $ takich, że

$$<br />
   \sum_{k=0}^N a_k\exp\left(\frac{2\pi*i*j*k}{N+1}\right)=b_k \quad j=0,\ldots,N.<br />
$$

Okazuje się, że rozwiązanie możemy wyrazić poprzez operator dyskretnej transformaty Fouriera, czy równoważnie jako mnożenie przez macierz $ \mathcal{F}_{N+1}=\frac{1}{N+1}(\omega^{k*l}_{N+1})_{k,l=0}^N $:

$$<br />
   \vec{a}=\mathcal{F}_{N+1} \vec{b}<br />
$$

gdzie $ \vec{a}=(a_k)_{k=0}^N,\vec{b}=(b_k)_{k=0}^N $, a $ \omega_{N+1}=\exp\left(\frac{-2\pi*i}{N+1}\right) $.

Mnożenie przez macierz $ \mathcal{F}_{N+1} $ można szybko wykonać z wykorzystaniem algorytmu szybkiej transformacji Fouriera, czyli FFT (Fast Fourier Transform). Odpowiednia wersja algorytmu FFT służy też szybkiemu mnożeniu przez macierz odwrotną do $ \mathcal{F}_{N+1} $: $ \mathcal{F}_{N+1}^{-1}=(N+1)\overline{\mathcal{F}}_{N+1}=(\overline{\omega}^{k*l}_{N+1})_{k,l=0}^N $.

Obliczając wartość wielomianu $ \sum_k a_k z^k $ w punktach $ z_j $, czy równoważnie wartości wielomianu trygonometrycznego $  \sum_{k=0}^N a_k\exp\left(i*k*x\right) $ w punktach $ x_j=\frac{2\pi*j}{N+1} $, to przyjmując
$ b_j=\sum_k a_k z_j^k $ otrzymujemy:

$$<br />
   \vec{b}=\mathcal{F}_{N+1}^{-1}\vec{a}.<br />
$$

Warto dodać, że często macierz DFT definiuje się bez czynnika $ \frac{1}{N+1} $, oczywiście wtedy macierz odwrotna też musi być odpowiednio przeskalowana.

Funkcja

 fft()

służy w octave'ie mnożeniu przez macierz $ \mathcal{F}_{N+1}. $
Jej najprostsze wywołanie to

 a=fft(b)

Sprawdźmy, czy funkcja fft() oblicza wartość DFT zgodnie z naszą definicją. Policzmy $ \mathcal{F}_4(4,0,0,0)^T $, powinniśmy otrzymać
wektor samych jedynek:

 a=fft([4;0;0;0])

a otrzymaliśmy:

a =

   4 + 0i
   4 + 0i
   4 + 0i
   4 - 0i

To oznacza, że funkcja fft() oblicza wartość mnożenia przez
$ (N+1)*\mathcal{F}_{N+1} $.

Z kolei funkcja ifft() powinna obliczać mnożenie przez macierz odwrotną do $ (N+1)*\mathcal{F}_{N+1} $.

Sprawdźmy, jak to działa:

 x=[1, -3, 4,5];
 a=fft(x);
 xx=ifft(a);
 x-xx
 norm(x-xx,2)

Powtórzmy to dla większego $ N $:

 x=rand(100);
 a=fft(x);
 xx=ifft(a);
 norm(x-xx,2)

Stwórzmy macierz $ (N+1)*\mathcal{F}_{N+1} $ i porównajmy szybkość mnożenia przez tę macierz wykonaną za pomocą standardowego operatora octave'a, czyli operatora *, z szybkością działania funkcji fft().

Stwórzmy najpierw funkcję tworzącą macierz $ (N+1)*\mathcal{F}_{N+1} $:

 function F=DFTmac(N=3)
  om=exp((-2*pi*i*(0:N))/(N+1));
  F=zeros(N+1,N+1);
  for k=0:N,
   F(k+1,:)=om.^k;
  endfor
 endfunction

Przetestujmy ją na początek dla $ N=3 $:

 F=DFTmac(3)
 F*F'
 a=rand(4);
 norm(fft(a)-F*a,2)

W tym przypadku wyniki się pokrywają.

Zgodnie z teorią, por. np. [Jankowska:1981:MN], koszt obliczenia DFT przy zastosowaniu algorytmu FFT to
$ O(n\log_2(n)) $, a standardowe mnożenia przez macierz $ \mathcal{F}_n $ kosztuje $ O(n^2) $.

A teraz testujemy szybkość:

 n=4*512
 F=DFTmac(n-1);
 x=rand(n);
 tic;y=F*x;t1=toc
 tic;yy=fft(x);t2=toc
 t1/t2

Na moim komputerze FFT było ponad $ 63 $ razy szybsze dla $ n=2048 $.

Policzmy jeszcze wykres realnego kosztu względem $ n $:

 n=t1=t2=zeros(100,1);
 n(1)=100;
 for k=1:100,
  F=DFTmac(n(k)-1);
  x=rand(n(k));
  tic;y=F*x;t1(k)=toc;
  tic;yy=fft(x); t2(k)=toc;
  n(k+1)=n(k)+10;
 endfor
 n=n(1:100);
 plot(n,t1,";mnozenie przez macierz;",n,t2,";fft;")

Na wynik chwilę musimy poczekać, ale potwierdza on teorię, że FFT jest wyraźnie szybszym algorytmem.

Interpolacja splajnowa

(#)
W przypadku interpolacji splajnowej rozpatrujemy dane węzły:

$$a=x_0,\ldots,x_N=b$$

na odcinku $ [a,b] $, oraz
funkcje splajnowe (splajny), czyli funkcje, które są odpowiedniej klasy gładkości (czyli są w $ C^k([a,b]) $, oraz obcięte do dowolnego pododcinka $ [x_k,x_{k+1}] $ dla $ k=0,\ldots,N-1 $, są wielomianami co najwyżej ustalonego stopnia.

W przypadku splajnów kubicznych rozważamy funkcje, które są klasy $ C^2 $ na $ [a,b] $, oraz są wielomianami kubicznymi na każdym pododcinku.

mo-SplineKub-simplest.jpg

mo-SplineKub-simplest2.jpg

Zadanie interpolacji splajnowej kubicznej polega na znalezieniu splajnu
kubicznego $ s $ takiego, że

$$<br />
 s(x_k)=y_k \quad k=0,\ldots,N<br />
$$

dla zadanych $ y_k $,
oraz spełniającego odpowiednie warunki brzegowe. Np. w przypadku splajnu hermitowskiego
$ s $ musi spełnić warunki brzegowe hermitowskie:

$$<br />
 s'(a)=dy_a \qquad  s'(b)=dy_b<br />
$$

dla zadanych dodatkowych dwóch wartości $ dy_a,dy_b $.
W przypadku splajnu naturalnego warunki brzegowe to

$$<br />
  s''(a)=s''(b)=0.<br />
$$

Rozpatruje się również splajny typu not-a-knot. W tym przypadku,
zamiast warunków brzegowych, dodaje się sztucznie dwa warunki ciągłości trzeciej pochodnej w węzłach $ x_1 $ i $ x_{N-1} $.

Z kolei splajny okresowe spełniają następujące warunki brzegowe okresowe:

$$<br />
  s^{(j)}(a)=s^{(j)}(b) \quad j=0,1,2.<br />
$$

W każdym z tych czterech przypadków splajn interpolacyjny jest wyznaczony jednoznacznie, por. np.
[Kincaid:2006:AN].

W octave'ie istnieje kilka funkcji związanych z interpolacją splajnami kubicznymi, czy ogólnie - z funkcjami wielomianowymi na pododcinkach:

  • spline() - funkcja służąca
    znalezieniu splajnu interpolacyjnego hermitowskiego, czy
    typu not-a-knot
  • ppval() - funkcja służąca obliczeniu wartości splajnu zadanego w formacie octave'a
  • mkpp() - tworzy funkcję wielomianową na pod-odcinkach (szczegóły help mkpp())
  • unmkpp() - ze struktury splajnu w octave'ie zwraca współczynniki wielomianów na pod-odcinkach (szczegóły help unmkpp())

Podstawową funkcją służącą
znalezieniu splajnu interpolacyjnego hermitowskiego, czy
typu not-a-knot jest funkcja spline().
Jej najprostsze wywołanie to

p=spline(x,y)

gdzie $ x $ to wektor wymiaru $ N $ z węzłami interpolacji splajnowej, a wektor $ y $ ma tę samą długość co $ x $ i zawiera wartości jakie ma przyjąć splajn w tych węzłach.
Druga możliwość to - przy takim samym wektorze węzłów $ x $, podanie wektora $ y $ o długości $ N+2 $. Wtedy pierwsza i ostatnia wartość wektora $ y $ to wartości pochodnych splajnu w końcowych węzłach $ x $. Pozostałe wartości $ y $ tzn. $ y_k $ dla $ k=2,\ldots,N+1 $ zawierają wartości splajnu w węzłach z $ x $, czyli są takie same jak w pierwszym przypadku.
Funkcja zwraca strukturę typu $ pp $, czyli w odpowiednim formacie octave'a splajnu kubicznego typu not-a-knot w pierwszym przypadku, a w drugim - splajnu hermitowskiego. Następnie, korzystając z funkcji ppval() można obliczyć wartość tego splajnu w punkcie, czy tablicy punktów.

Policzmy splajn który w węzłach $ -2,1,0,1,2 $ przyjmie wartości $ 1,0,0,0,1 $ i narysujmy jego wykres, por. rysunek  [link]:

x=-2:2;
y=[1,0,0,0,1];
pp=spline(x,y);
z=linspace(-2,2);
plot(z,ppval(pp,z),";splajn n-a-k;",x,y,"r+");

Korzystając z tej samej funkcji, stwórzmy splajn hermitowski przyjmując, że jego pochodne w końcach wynoszą zero. Następnie narysujmy wykresy obu splajnów na odcinku $ [-2,2] $, por. rysunek  [link], oraz blisko lewego końca, por. rysunek  [link]:

x=-2:2;
y=[1,0,0,0,1];
pp=spline(x,y);
yh=[0,y,0];
ph=spline(x,yh);
z=linspace(-2,2);
plot(z,ppval(ph,z),";splajn hermitowski;",...
z,ppval(pp,z),";splajn n-a-k;",x,y,"r+");
pause(2);
z=linspace(-2,-2+0.1);
plot(z,ppval(ph,z),";splajn hermitowski;",...
z,ppval(pp,z),";splajn n-a-k;",-2,1.3);

Widać, że splajny są różne, oraz że splajn hermitowski ma pochodną równą zero w lewym i prawym końcu.
mo-SplineKub-simplest3.jpg

Czy w octave'ie można wyznaczyć w prosty sposób inne splajny, np. naturalny lub okresowy?

Okazuje się, że tak, ale trzeba wykorzystać funkcję z rozszerzenia octave'a, czyli octave-forge'a: csape().

Jej wywołanie to

pp=csape(x,y,cond,valz)

gdzie $ x,y $ to wektory długości $ N $ z węzłami i wartościami splajnu w węzłach, a cond przyjmuje wartości:

  • 'variational' w przypadku splajnu interpolacyjnego kubicznego naturalnego
  • 'complete' w przypadku splajnu interpolacyjnego kubicznego hermitowskiego,
    wartości pochodnych w końcach są podane w parametrze valc
  • 'not-a-knot' w przypadku splajnu interpolacyjnego kubicznego typu \textsl{not-a-knot}
  • 'periodic'
    w przypadku splajnu interpolacyjnego kubicznego okresowego
  • 'second'
    w przypadku splajnu interpolacyjnego kubicznego z ustalonymi wartościami drugiej pochodnej w końcach, zgodnymi z tym, co podano w parametrze valc

Funkcja zwraca strukturę typu $ pp $, czyli dane splajnu kubicznego w formacie octave'a.

Porównajmy wyniki tej funkcji z wynikiem funkcji spline() dla danych z naszego prostego przykładu i splajnu typu not-a-knot. Policzymy dyskretną normę maksimum (na siatce równomiernej 300 punktów na $ [-2,2] $ z wyników obu funkcji:

x=-2:2;
y=[1,0,0,0,1];
pp=spline(x,y);
pp1=csape(x,y,'not-a-knot');
z=linspace(-2,2,300);
norm(ppval(pp,z)-ppval(pp1,z),'inf')

Otrzymaliśmy zero.

Na jednym wykresie narysujmy wykresy trzech splajnów interpolacyjnych dla $ N=6 $ z tymi samymi warunkami interpolacyjnymi, ale z różnymi warunkami brzegowymi: hermitowskim, not-a-knot i naturalnym.

N=4;
x=-2:2;
y=[1,0,0,0,1];
pp=spline(x,y);
ppn=csape(x,y,'variational');
yh=[0,y,0];
pph=spline(x,yh);
z=linspace(-2,2);
plot(z,ppval(pph,z),";spl. hermit.;",...
z,ppval(pp,z),";spl. n-a-k;",...
z,ppval(ppn,z),";spl. natur.;",x,y,"r+");

Wszystkie trzy splajny są różne w tym przypadku, por.
rysunek  [link].
mo-SplineKub-diff4.jpg

mo-SplineKub-noteqknots6.jpg

Oczywiście obie funkcje umożliwiają dowolny wybór węzłów, np. weźmy węzły
$ [-2,-1,2,3,4] $ z tymi samymi wartościami i stwórzmy splajny obu typów, por. rysunek  [link]:

x=[-3,2,3,4];
y=[1,0,0,0,1];
pp=spline(x,y);
yh=[0,y,0];
ph=spline(x,yh);
z=linspace(-3,4);
plot(z,ppval(ph,z),";splajn hermitowski;",...
z,ppval(pp,z),";splajn n-a-k;",x,y,"r+");

mo-SplineKub-err5.jpg

Popatrzmy na błędy aproksymacji w normie supremum.
Ustalmy, że interpolujemy splajnami kubicznymi znaną funkcję gładką $ f(x) $ na czterech węzłach równoodległych na odcinku $ [a,b] $.
Jeśli $ pp $ to struktura typu $ pp $ opisująca splajn interpolujący $ f $, to dyskretną normę maksimum różnicy między $ f $, a splajnem $ s $, np. na siatce o tysiącu punktach, możemy obliczyć komendą:

z=linspace(a,b,1000);
y=f(z);
s=ppval(pp,z);
errmax=norm(s-y,'inf')

Założyliśmy, że funkcja $ f $ jest zaimplementowana wektorowo, tzn. że wywołanie w octave'ie f(z) dla $ z $ wektora zwróci wektor $ (f(z(k))) $. W przeciwnym razie należałoby użyć pętli do wyznaczenia wektora $ y $:

z=linspace(a,b,1000);
for k=1:length(z),
 y(k)=f(z(k));
endfor
s=ppval(pp,z);
errmax=norm(s-y,'inf')

W poniższym kodzie obliczymy przybliżoną normę maksimum na odcinku $ [-1,2] $ pomiędzy
$ f(x)=\sin(x) $, a jej dwoma splajnami interpolacyjnymi kubicznymi: naturalnym i typu not-a-knot na czterech węzłach równoodległych:

a=-1;
b=2;
f=@sin;
N=4;
x=linspace(a,b,4);
y=f(x);
pp=spline(x,y);
ppn=csape(x,y,'variational');
z=linspace(a,b,1000);
fz=f(z);
pz=ppval(pp,z);
pn=ppval(ppn,z);
err=norm(fz-pz,'inf')
errn=norm(fz-pn,'inf')

Błąd dla splajnu naturalnego był większy:
errn =  0.048190, niż dla splajnu typu not-a -knot:
err =  0.025120, co widać też na rysunku  [link].

Jako zadanie pozostawiamy policzenie błędów dla innych funkcji, węzłów i typów splajnów kubicznych interpolacyjnych.

Obliczanie całek. Kwadratury

W tym rozdziale zajmiemy się zadaniem obliczenia przybliżonego całek postaci:

$$<br />
 \int_a^b f(t) \:dt,<br />
$$

dla funkcji $ f $,
czy ogólniej:

$$<br />
  \int_a^b f(t)\rho(t) \:dt,<br />
$$

dla $ \rho $ danej wagi.

Funkcja octave'a quad()

Do obliczania przybliżonego całek jednowymiarowych tzn. typu

$$c=\int_a^b f(t)\: d t $$

służy funkcja octave'a:

 quad()

Jej najprostsze wywołanie to:

 c=quad(FCN,a,b),

gdzie $ a,b $ to początek i koniec odcinka, a $ FCN $ to wskaźnik (uchwyt) do funkcji
octave'a postaci:

 function y=f(x)

  #komendy octave'a

  y=.....;

endfunction

Jeśli funkcja, której całkę chcemy obliczyć jest już wcześniej zdefiniowana, to uchwyt do niej zwraca operator octave'a: @, np. uchwyt do funkcji $ \sin $ zwróci komenda: @sin.

Na początek kilka prostych przykładów całek, które sami potrafimy obliczyć analitycznie:

  • całka z $ \sin(x) $ na $ [-1,1] $
     quad(@sin,-1,1)

    Otrzymaliśmy zero - jak tego oczekiwaliśmy.

  • $ \int_0^\pi \sin(x) \: d x=2 $
     quad(@sin,0,pi)

    Otrzymaliśmy wynik zgodny z oczekiwaniami.

  • całka z funkcji mało regularnej $ \int_0^1 \sqrt{t} \: d t $:
     c=quad(@sqrt,0,1)

    Sprawdźmy c-2/3, otrzymaliśmy $ ans =  2.2204e-16 $, czyli zero na poziomie błędu zaokrągleń.

mo-quad-eps1.jpg
Czy funkcja quad() zawsze daje dobre wyniki?
Rozpatrzmy prosty przykład
całki z funkcji z parametrem

$$f_\epsilon(t)=\epsilon^{-1}f1(\frac{t}{\epsilon})$$

na $ [-1,1] $, por. rysunek  [link], dla

$$<br />
  f1(t)=\left\{<br />
      \begin{array}{lr}<br />
       1-|t|  & x\in(-1,1)\\<br />
       0 & \mathrm{w \; przeciwnym\; przypadku}<br />
      \end{array}<br />
\right.<br />
$$

mo-quad-eps2.jpg

Całka powinna być zawsze równa jeden.
Zdefiniujemy najpierw funkcję $ f1 $:

 function y=f1(x)
  y=1-abs(x);
  y=y.*(y>0);
 endfunction
 c=quad(@f1 ,-1,1)

Otrzymaliśmy jeden. Powtórzmy obliczenia dla $ \epsilon=10^{-k} $ dla $ k=1,3,5,7 $:

 M=6;
 c=epsi=zeros(M+1,1);
 epsi(1)=1;
 for k=1:M,
   c(k)=quad(@(x) epsi(k)*f1(epsi(k)*x) ,-1,1);
   epsi(k+1)=epsi(k)*10;
 endfor

Od pewnego momentu obliczone przybliżenia całek są równe zero zamiast jeden.

Wychwyćmy ten moment dokładniej.
Powtórzmy obliczenia dla $ \epsilon $ z $ [1e-3,1e-2] $.

 epsi=linspace(100,1000,300);
 M=length(epsi);
 c=zeros(M,1);
 for k=1:M,
   c(k)=quad(@(x) epsi(k)*f1(epsi(k)*x) ,-1,1);
 endfor
 plot(epsi,c)

Jest to dość gwałtowna zmiana, por. rysunek  [link].

Dlaczego tak jest? Oczywiście funkcja octave'a quad(() może obliczyć przybliżenie całki korzystając z co najwyżej skończonej ilości wartości funkcji na odcinku całkowania. Nośnik funkcji $ f_\epsilon $ jest bardzo
mały dla $ \epsilon\approx 1e-3 $ względem odcinka całkowania $ [-1,1] $, stąd być może wszystkie obliczone wartości były równe zero.

Można się domyśleć, że funkcje octave'a działają na zasadzie czarnej skrzynki: podajemy parametry do funkcji, a dana funkcja, w tym przypadku quad() powinna zwrócić przybliżenie całki.

Spróbujmy obliczyć tę całkę
dzieląc odcinek $ [-1,1] $ na małe pododcinki i całkując po nich tzn.

$$<br />
 \int_a^b f(t)\: d t=\sum_{k=0}^{N-1} \int_{x_k}^{x_k+h} f(t)\: d t \qquad h=(b-a)/N;\quad x_k=a+k*h ,<br />
$$

a każdą całkę $ \int_{x_k}^{x_k+h} f(t)\: d t $ możemy obliczyć przy pomocy funkcji
octave'a quad().

 epsi=1000;
 N=100;
 h=2/N;
 xk=-1;
 g=@(x) epsi*f1(epsi*x);
 c=0;
 for k=0:(N-1),
   c+=quad(g ,xk,xk+h);
   xk+=h;
 endfor
 abs(c-1)

Otrzymaliśmy błąd $ abs(c-1) $ równy ans =  6.8505e-12,
podczas gdy
abs(quad(g ,-1,1)-1) zwraca błąd równy jeden.

Oczywiście znając nośnik funkcji możemy policzyć całkę po jej nośniku.
Otrzymamy:

 epsi=1000;
 g=@(x) epsi*f1(epsi*x);
 c=quad(g ,-1e-3,1e-3)
 abs(c-1)

Błąd jest równy ans =  1.1102e-16, czyli jest na poziomie błędu zaokrągleń.

Również dla funkcji o dużej zmienności, np. silnie oscylujących, funkcja
octave'a quad(() może zwrócić zły wynik.

Spróbujmy scałkować np. $ \sin(999*x) $ na odcinku $ [0,\pi] $ i porównajmy z dokładnym wynikiem $ \int_0^\pi \sin(999*x)d t=1e-3*(1-\cos(999*\pi)) $.

 a=999;
 g=@(x) sin(a*x);
 c=quad(g ,0,pi)
 abs(c-(1/a)*(1-cos(a*pi)))

Tu funkcja quad() wydrukowała na ekranie ostrzeżenie:

 ABNORMAL RETURN FROM DQAGP

Warto wspomnieć, że funkcja quad() oblicza również całki po odcinkach nieograniczonych. Za wartości $ a,b $ w wywołaniu quad(f,a,b) możemy przyjąć
inf lub odpowiednio -inf, np. możemy scałkować funkcję $ \exp(\frac{x^2}{2}) $ po całej prostej rzeczywistej, czy jakiejś półprostej:

 g=@(x) exp(-0.5*x*x);
 c1=quad(g ,-inf,0)
 c2=quad(g,0,inf)
 c3=quad(g,-inf,inf)
 c3-c1-c2

Czy wynik jest zgodny z oczekiwaniami?

Kwadratury złożone

Oczywiście w octave'ie możemy zaimplementować
najprostsze kwadratury złożone, np. najprostszą kwadraturę złożoną:

$$<br />
  \int_a^b f(t) \: d t P_N(f)\approx \frac{b-a}{N}\sum_{k=1}^N f(a+k*h)<br />
$$

Implementacja w octave'ie tej kwadratury to:

function c=prostakw(f,a,b,N=300)
 #kwadratura prostokatow prawostronna
 #Input: a,b, konce przedzialu calkowania
 #N- ilosc punktow - domyslnie 100
 #f- wskaznik do funkcji -
 #     zwracajacej wektor wartosci dla wektora danych
 #Output: c - przyblizenie calki
 h=(b-a)/N;
 x=(1:N)*h;
 y=f(x);
 c=h*sum(y);
endfunction

Przetestujmy tę kwadraturę. Na początku na jednomianach:

prostakw(@(x) x,0,1)-0.5
prostakw(@(x) x.*x,0,1)-1/3

Wynik jest poprawny. Jeżeli zwiększymy $ N $, to być może uzyskamy lepszy wynik:

prostakw(@(x) x,0,1,1000)-0.5
prostakw(@(x) x.*x,0,1,1000)-1/3

Błąd jest na poziomie $ 10^{-4} $.

Dla funkcji klasy $ C^1([a,b]) $
błąd można oszacować przez

$$<br />
| \int_a^b f(t) \: d t - P_N f|\leq \|f'\|_{\infty,[a,b]}h=O(N^{-1}),<br />
$$

przy czym dla funkcji $ f(x)=x $ w tym oszacowaniu widzimy równość.

Można pokazać, że dla dostatecznie gładkiej funkcji zachodzi:

$$<br />
  \int_a^b f(t) \: d t - P_N f=-f'(\xi)h<br />
$$

dla pewnego punktu $ \xi\in(a,b) $, a dla bardziej regularnych funkcji

$$<br />
  \int_a^b f(t) \: d t - P_N f=Ch + O(h^2)<br />
$$

dla pewnej ujemnej stałej $ C $.

\[<br />
\begin{tabular}{|c|c|c|c|}<br />
\hline<br />
N & h & $E_N$ & $E_N/E_{N/2}$\\<br />
\hline<br />
10 & 1.00e-01 & 5.17e-02 & 0.00000 \\<br />
 20 & 5.00e-02 & 2.54e-02 & 2.03279 \\<br />
 40 & 2.50e-02 & 1.26e-02 & 2.01653 \\<br />
 80 & 1.25e-02 & 6.28e-03 & 2.00830 \\<br />
 160 & 6.25e-03 & 3.13e-03 & 2.00416 \\<br />
 320 & 3.13e-03 & 1.56e-03 & 2.00208 \\<br />
 640 & 1.56e-03 & 7.82e-04 & 2.00104 \\<br />
 1280 & 7.81e-04 & 3.91e-04 & 2.00052 \\<br />
\hline<br />
\end{tabular}<br />
\]

Badanie rzędu kwadratury $ P_N f $ dla całki z funkcji $ f(x)=x^2 $ na $ [0,1] $. Błąd $ E_N=|\int_0^1f - P_Nf| $. (#)

Przetestujmy teraz tę kwadraturę dla ustalonej analitycznej funkcji $ f $ i podwajanych
wartości $ N $, tzn. obliczymy $ P_Nf $ dla $ N=N_0,2*N_0,\ldots $ i następnie
policzymy:

$$<br />
   \frac{E_k}{E_{k+1}}=\frac{C*h+O(h^2)}{C*h/2+O(h^2)} \approx 2,<br />
$$

dla $ E_k=| \int_a^b f(t) \: d t - P_{2^kN_0} f| $.

\[<br />
\begin{tabular}{|c|c|c|c|}<br />
\hline<br />
N & h & $E_N$ & $E_N/E_{N/2}$\\<br />
\hline<br />
10 & 3.14e-01 & 1.65e-02 &  \\<br />
 20 & 1.57e-01 & 4.11e-03 & 4.00495 \\<br />
 40 & 7.85e-02 & 1.03e-03 & 4.00123 \\<br />
 80 & 3.93e-02 & 2.57e-04 & 4.00031 \\<br />
 160 & 1.96e-02 & 6.43e-05 & 4.00008 \\<br />
 320 & 9.82e-03 & 1.61e-05 & 4.00002 \\<br />
 640 & 4.91e-03 & 4.02e-06 & 4.00000 \\<br />
 1280 & 2.45e-03 & 1.00e-06 & 4.00000 \\<br />
\hline<br />
\end{tabular}<br />
\]

Badanie rzędu kwadratury $ P_N f $ dla całki z funkcji $ f(x)=\sin(x) $ na $ [0,\pi] $. Błąd $ E_N=|\int_0^\pi f - P_Nf| $. (#)

 N=10; M=8;
 #c=2; #znana wartosc calki
 e=0;
 for k=1:M,
  kw=prostakw(f,a,b,N);
  ep=e;
  e=abs(c-kw);
  printf("[%d] e=%g ep/e=%6.5f  \n",N,e,ep/e);
  N*=2;
 endfor

Wyniki dla funkcji $ f(x)=x^2 $ i odcinka $ [0,1] $ są zawarte w Tabeli  [link], a dla funkcji $ \sin(x) $ i
całki z odcinka $ [0,\pi] $ w Tabeli  [link].
Dlaczego dla $ \sin(x) $ wyniki wskazują na błąd $ E_N\approx O(h^2) $?

\[<br />
\begin{tabular}{|c|c|c|c|}<br />
\hline<br />
N & h & $E_N$ & $E_N/E_{N/2}$\\<br />
\hline<br />
10 & 1.00e-01 & 4.17e-02 &  \\<br />
 20 & 5.00e-02 & 2.09e-02 & 1.99085 \\<br />
 40 & 2.50e-02 & 1.05e-02 & 1.99544 \\<br />
 80 & 1.25e-02 & 5.25e-03 & 1.99772 \\<br />
 160 & 6.25e-03 & 2.63e-03 & 1.99886 \\<br />
 320 & 3.13e-03 & 1.31e-03 & 1.99943 \\<br />
 640 & 1.56e-03 & 6.57e-04 & 1.99972 \\<br />
 1280 & 7.81e-04 & 3.29e-04 & 1.99986 \\<br />
\hline<br />
\end{tabular}<br />
\]

Badanie rzędu kwadratury $ P_N f $ dla całki z funkcji $ f(x)=\sin(x) $ na $ [0,1] $. Błąd $ E_N=|\int_0^1 f - P_Nf| $. (#)

Zmieńmy odcinek na $ [0,1] $ i wyniki będą zgodne z oczekiwanymi, por. Tabela  [link].

Normy całkowe

Zastanówmy się, jak w sposób przybliżony wyznaczyć normę typu $ L^p $ na odcinku $ [a,b] $ dla zadanej funkcji $ f(x) $ zdefiniowanej w octave'ie jako

 function y=f(x)

Spróbujmy zdefiniować funkcję octave'a, która korzystając z quad()
wyznaczy przybliżenie normy $ \|f\|_{L^p(a,b)}=\left(\int_a^b|f(t)|^p\: d t<br />
\right)^{1/p} $.
Jako parametry tej funkcji chcemy podawać wskaźnik do funkcji typu y=f(x), $ p $, oraz końce odcinka $ [a,b] $.
Żeby wywołać funkcję quad() musimy podać funkcję octave'a jednoargumentową $ x\mapsto |f(x)|^p $ jako pierwszy argument.
Trzeba ją odpowiednio zdefiniować wewnątrz naszej funkcji. Najwygodniejszym sposobem jest użycie mechanizmu funkcji anonimowej, por Rozdział  [link]:

fp=@(x) (abs(f(x))^p;

oczywiście $ p $ musi być wcześniej zdefiniowane.
Tak więc można by napisać tę funkcję następująco:

function n=Lpnorm(f,a,b,p=2)
 fp=@(x) abs(f(x))^p;
 n=quad(fp,a,b);
 n=n^(1/p);
endfunction

Sprawdźmy dla $ p=2 $ i prostych funkcji
np.:

Lpnorm(@(x) 1,0,1)-1
Lpnorm(@(x) x,0,1)-1/sqrt(3)
Lpnorm(@cos,0,pi)-sqrt(pi/2)

Wyniki są poprawne.

Mając zaimplementowaną taką funkcję, możemy np. przeprowadzić badanie rzędów zbieżności interpolacji splajnowej w normie $ L^2 $, zamiast w normie maksimum, tak jak w Rozdziale  [link].

Rozwiązywanie równań nieliniowych

W tym rozdziale zajmiemy się metodami rozwiązywania
równań nieliniowych skalarnych.

Interesuje nas znalezienie zera nieliniowej funkcji $ f:[a,b]\rightarrow \mathbb{R} $:

$$<br />
 f(x)=0.<br />
$$

Przetestujemy kilka metod.

Funkcja octave'a fzero()

Na początek zapoznamy się z funkcją octave'a służącą do znajdowania zer równania nieliniowego, czyli fzero().

Funkcja posiada dwa wywołania w pierwszym:

 x=fzero(f,ap)

$ f $ jest wskaźnikiem (uchwytem) do funkcji

 function y=f(x)

a $ app=(a,b) $ jest wektorem z dwoma liczbami lokalizującymi zero funkcji, tzn. że to zero leży pomiędzy $ a $ i $ b $.

Przetestujmy tę funkcję z tym wywołaniem na kilku prostych przykładach:

fzero(@(x) x^5,[-1,2])
fzero(@(x) x*x-2,[1,2])-sqrt(2)
fzero(@sin,[-1,2])
fzero(@sin,[3,4])-pi

Funkcja działa. Przetestujmy ją w sytuacji kiedy zer może być kilka np.:

fzero("sin",[-4,5])

Funkcja znalazła pierwiastek $ -\pi $, ale może nie zadziałać:

fzero("sin",[-4,1])

Zwracając:

error: fzero: not a valid initial bracketing

Nietrudno się domyśleć, że dla tej funkcji właściwy przedział to taki, że funkcja w końcach przyjmuje przeciwne znaki.

Kolejnym możliwym wywołaniem tej funkcji jest podanie $ ap $ jako liczby -przybliżenia zera:

fzero(@(x) x*x-2,2)-sqrt(2)
fzero(@sin,1)
fzero(@sin,3)-pi

To wywołanie też działa.
Inną funkcją octave'a służącą rozwiązywaniu równań nieliniowych jest
fsolve(). Jej wywołanie jest takie samo, jak druga wersja wywołania funkcji fzero(), tzn. podajemy wskaźnik (uchwyt) do funkcji i liczbę będącą przybliżeniem pierwiastka, np.:

fsolve(@(x) x*x-4,2)

Przetestujmy na kilku innych przykładach:

fsolve(@(x) x*x*x-27,5)-3
fsolve(@sin,1)
fsolve(@sin,3)-pi

Wyniki są poprawne.

Funkcja fsolve() jest bardziej ogólna od fzero(), może służyć znajdowania zer układów równań nieliniowych, ale to wykracza poza standardowwy zakres
wykładu z Matematyki Obliczeniowej.

Metoda bisekcji

Najprostszą metodą znajdowania zera jest metoda bisekcji.
Dla funkcji ciągłej $ f $ - jeśli $ f(a)*f(b)<0 $ - to istnieje zero pomiędzy punktami $ a,b $.

Idea metody bisekcji (połowienia odcinka), to
przyjęcie za przybliżenie zera $ x=(a+b)/2 $ i sprawdzenie znaku $ f $ w $ x $. Albo $ f(x) $ jest równe zero - wtedy otrzymaliśmy, że $ x $ jest szukanym zerem. W przeciwnym razie zastępujemy jeden z końców odcinka przez $ x $ - ten w którym znak wartości $ f $ w tym końcu jest ten sam co dla $ f(x) $. Dalej postępujemy analogicznie z nowym odcinkiem o długości równej połowie poprzedniego.

Poniżej widzimy prostą implementację metody bisekcji:

function x=bisekcja(f,a,b,eps=1e-4,itmax=200)
#[x]=bisekcja(f,a,b,eps=1e-4,itmax=200)
#f- wskaznik do funkcji
#a,b poczatek i koniec odcinka
#eps - tolerancja met zatrzymuje sie o ile (b-a)<eps
#itmax - max ilosc iteracji
#Output:
#x- przyblizone zero
#Zakladamy, ze uzytkownik poda dwa punkty a<b
# w ktorych f(a)f(b)<0
ya=f(a);
yb=f(b);
it=0;
x=(a+b)/2;
while( (b-a>eps) && (it<itmax))
y=f(x);
if(sign(y)==sign(ya))
 ya=y;
 a=x;
else
 yb=y;
 b=x;
end
 x=(a+b)/2;
endwhile
endfunction

Przetestujmy dla kilku funkcji:

bisekcja(@(x) x*x-2,1,2)-sqrt(2)

Wynik jest poprawny.

Jeśli $ |b-a| $ jest duże to zbieżność jest dość wolna:

tic;bisekcja(@(x) x*x-2,1,1e30)-sqrt(2);toc

Czas potrzebny do obliczenia zera wydaje się mały, ale wyobraźmy sobie, że musimy rozwiązać kilka milionów takich równań.
Spróbujmy znaleźć zero dla funkcji bardziej skomplikowanej $ g(x)=\int_0^x \exp(-0.5*t^2)\:d t=0 $:

f=@(x) exp(-x*x/2);
g=@(x) quad(f,0,x)-1;
x=bisekcja(g,1,2)
g(x)

Metoda Newtona

Kolejną bardzo ważną metodą jest metoda Newtona:

$$<br />
  x_{n+1}=x_n - f(x_n)/f'(x_n) \qquad n\geq 0<br />
$$

poprawnie określona, o ile $ f'(x_n)\not = 0 $.

Zaimplementujmy w octave'ie metodę Newtona. Poniżej podajemy kod funkcji octave'a
szukania zer za pomocą metody Newtona:

function [x,fval,info]=newton(f,df,x0,...
                     epsr=1e-7,epsa=1e-9,maxit=30)
#Metoda Newtona skalarna
##Input:
#f - function handle do funkcji postaci y=f(x),
#ktorej zera szukamy przy czym zwracana wartosc
#y to wartosc f(x)
#df - function handle do pochodnej funkcji postaci
#y=df(x) ktorej zera szukamy
#x - startowe przyblizenie
##Pozostale parametry opcjonalne
#epsr = wzgledna tolerancja
#epsa - bezwgledna tolerancja
#maxit - maksymalna ilosc iteracji
#Output:
#y - przyblizenie zera
#fval - wartosc f(x)
#info - wynik 0 - zbieznosc;
#             1   brak zbieznosci -
#                  przekroczyl max ilosc iteracji
#             2 - df(x) - zero nie mozna policzyc
#                        kolejnego przyblizenia
# iteracje zatrzymuja sie jesli
# zachodzi jedna z dwu
#   1/|f(x)|>epsa + epsr*|f(x0)|
#   lub
#   2/ilosc iteracji = maxit
## Przyklad:
# function y=f(x)
#  y=x*x-2;
# endfunction
## Wywolujemy:
# y=newton(@f,2)
# [y,fy,info]=newton(@f,4,1e-2,1e-3) itd

fx=f(x0);
fval=fx(1);
it=0;
info=0;
eps=0.5*(epsr*abs(fval)+epsa);

x=x0;
while((abs(fval)>eps ) && (it<maxit))
#liczymy pochodna
 dfx=df(x);
 if(abs(dfx)>1e-12)
  x=x - fval/dfx;
 else
  info=2;
  printf("Pochodna df w punkcie x_%d=%g zero\n",it,x);
  return;
 endif
 it++;
 fx=f(x);
 fval=fx(1);
endwhile
if((abs(fval)>eps ))
 info=1; #brak zbieznosci
 printf("Brak zbieznosci! Po ...
%d iteracjach  |f(%g)|=%g>tol=%g!\n",it-1,x,fval,eps);
endif
endfunction

Zbadajmy zbieżność metody Newtona dla modelowej funkcji $ f(x)=x^2-4 $, której pierwiastkami są $ -2,2 $.

newton(@(x) x*x-4,@(x) 2*x,4)-2

Teraz wywołamy tę funkcję ze wszystkimi zwracanymi wartościami:

[x,fval,info]=newton(@(x) x*x-4,@(x) 2*x,4)

Przetestujmy, czy metoda zbiega dla $ x_0=100 $:

[x,fval,info]=newton(@(x) x*x-4,@(x) 2*x,100)

czy dla $ x_0 $ jeszcze bardziej oddalonych od pierwiastków:

[x,fval,info]=newton(@(x) x*x-4,@(x) 2*x,1000)
[x,fval,info]=newton(@(x) x*x-4,@(x) 2*x,0.01)
[x,fval,info]=newton(@(x) x*x-4,@(x) 2*x,1e-9)
[x,fval,info]=newton(@(x) x*x-4,@(x) 2*x,1e6)

W ostatnim przypadku ilość iteracji w warunku stopu okazała się niewystarczająca.

Możemy też przetestować rząd zbieżności dla równania, którego rozwiązanie znamy, np. rozpatrzmy ostatnie równanie $ x^2=4 $:
Poniższy kod testuje rząd zbieżności metody Newtona dla dodatniego pierwiastka
funkcji $ f(x)=x^2-4 $:

M=6;
x=4;
e=0;
it=0;
for k=1:M,
 x=x-(x*x-4)/(2*x);
 it++;
 ep=e;
 e=x-2;
 if(it>1)
  printf("[%d] e/ep=%3.4e e/(ep*ep)=%3.4e \n",...
   it,e/ep,e/(ep*ep));
 endif
endfor

mo-nonlin-wyk.jpg

Otrzymaliśmy, że rząd zbieżności dla tej funkcji jest równy dwa, co jest zgodne z teorią, która mówi, że jeśli funkcja jest klasy $ C^2 $ na otoczeniu swojego pierwiastka $ x^* $, oraz wartość pochodnej funkcji w pierwiastku jest różna od zera to metoda jest lokalnie zbieżna kwadratowo (zbieżna z rzędem równym dwa): tzn. istnieje takie $ \epsilon>0 $ i stała $ C\geq 0 $, że jeśli $ x_0 $ startowe przybliżenie metody Newtona spełnia: $ |x^*-x_0|\leq \epsilon $, to

$$<br />
     |x_{n+1}-x^*|\leq C |x_n- x^*|^2,<br />
$$

por. np. [Kincaid:2006:AN].

Odwracanie funkcji

Podamy jeden przykład zastosowania metod rozwiązywania równań nieliniowych do znajdowania przybliżonych wartości funkcji odwrotnej do danej.
Załóżmy, że potrafimy policzyć wartości funkcji, i chcemy znaleźć wartości funkcji odwrotnej do danej funkcji, aby np. narysować wykres funkcji odwrotnej (zakładamy, że na danym odcinku ona istnieje). Zazwyczaj niestety nie znamy wzoru analitycznego na funkcję odwrotną.

mo-nonlin-wykod.jpg
Jak to zrobić z zastosowaniem funkcji fzero()?

mo-nonlin-wykod2.jpg
Gdyby chodziło o sam wykres możemy
zastosować proste odbicie:

x=linspace(a,b);
y=f(x);
plot(x,y,";wykres y=f(x);");
plot(y,x,";wykres x=f^{-1}(y);");

Na rysunku  [link] widać wykres $ f(x)=2*x+\sin(x) $ na $ [0,3] $, a na
rysunku  [link] wykres funkcji do niej odwrotnej.

Oczywiście wykres funkcji odwrotnej wygląda zgodnie z oczekiwaniami.
Rozpatrzmy inną funkcje:
$ f(x)=\tan(x) $ na $ [-1.57,1.57] $ i popatrzmy na wykresy funkcji odwrotnej otrzymanej w analogiczny sposób i wykres rzeczywistej funkcji odwrotnej, której wzór analityczny w tym przypadku znamy: $ (\tan)^{-1}(y)=\mathrm{atan}(y) $:
rysunek  [link]. Widać, że wykres funkcji odwrotnej jest trochę inny od wykresu uzyskanego naszym sposobem. Aby uzyskać wykres prawidłowy powinniśmy znać bardzo dokładne przybliżenia wartości funkcji odwrotnej w punktach siatki.

Aby obliczyć wartość $ f^{-1}(y) $ dla konkretnego $ y $ należy rozwiązać równanie:

$$<br />
  g(x)=y-f(x)=0,<br />
$$

co w praktyce można w sposób przybliżony dokonać przy użyciu np. funkcji fzero() w octave'ie.

mo-nonlin-wykod3.jpg

Poniższy kod oblicza wartości funkcji odwrotnej do $ f(x)=\tan(x) $ na siatce równomiernych punktów na $ [-f(1.57),f(1.57)] $ (korzystamy z tego, że funkcja jest nieparzysta i obliczamy jej wartości tylko dla argumentów dodatnich):

f=@tan;
c=f(-1.57);d=-c;
N=1000;
y=linspace(0,d,N);
x=zeros(N,1);
x(1)=0; #tan(0)=0
for k=2:N,
  x0=x(k-1);
  g=@(x) y(k) - f(x);
  x(k)=fsolve(g,x0);
endfor
y=[y-d y];x=[-flipud(x);x];
plot(y,x,";oblicz. f odwrotna;",y,atan(y),";atan(y);")
pause(2);
plot(y,x'-atan(y),";wykres bledu;")

%Użyliśmy funkcji fsolve(), która jet wywoływana analogicznie jak fzero() w drugim typie wywołania, tzn. drugi argument jest przybliżeniem zera funkcji $ f(x) $.

mo-nonlin-wykod4.jpg

Na rysunku  [link] widać, że wykresy funkcji odwrotnej oraz obliczonej przez nas pokrywają się. Potwierdza to także wykres błędu, por. rysunek  [link].

Rozwiązywanie równań liniowych

W tym rozdziale zajmiemy się algorytmami związanymi z algebrą liniową, a dokładnie rozwiązywaniem układów równań liniowych.

Interesuje nas zadanie znalezienia rozwiązania układu równań liniowych:

$$<br />
\label{eq:nal:url}<br />
 As=f<br />
$$

gdzie $ A $ to macierz nieosobliwa $ n \times n $ i $ f $ to wektor prawej strony.

Operator backslash

Podstawowym operatorem rozwiązującym ( [link]) jest operator:
octave'a \. Jego zastosowanie jest bardzo proste

 x=A\f

Przetestujmy ten operator dla kilku macierzy:

A=[1 1;-1 1];
s=[1; 3];
f=A*s;
x=A\f;
x-s
norm(x-s,2)
norm(A*x-f,2)

Weźmy macierz prawie osobliwą:

A=[1 1+2*eps;1 1];
s=[1; 3]
f=A*s;
x=A\f
x-s
norm(x-s,2)
norm(A*x-f,2)

Wywołanie operatora A\f zwróciło ostrzeżenie i dało niepoprawny wynik. Zauważmy, że norma residualna $ \|A*x-f\|_2 $ jest mała, podczas gdy rozwiązanie jest zupełnie złe: $ \|x-s\|_2>1 $.

Nasza macierz została specjalnie tak dobrana, aby odpowiedni algorytm używany przez octave'a nie zwrócił poprawnego wyniku - ale w praktycznych zastosowaniach istnieją takie macierze, dla których błędy zaokrągleń powodują, że rozwiązywanie układów równań liniowych z nimi jest obarczone ogromnym błędem.

Operator \ może też w jednym wywołaniu rozwiązać układ równań liniowych z prawą stroną będącą macierzą, tzn.:

$$<br />
  AX=F<br />
$$

gdzie $ A $ to macierz nieosobliwa $ n \times n $, $ X,F $ to macierze $ n \times m $.

Wywołujemy po prostu

X=A\F

Taki układ jest równoważny $ m $ układom z macierzą $ A $ i prawymi stronami - kolumnami macierzy $ F $, jako że:

$$<br />
 AX=(Ax_1,Ax_2,\ldots,Ax_m)=(f_1,\ldots,f_m)=F.<br />
$$

Pętle w octave'ie działają wolno, więc - jeśli to możliwe - należy unikać ich wywoływania. Na poniższym przykładzie porównamy czas działania operatora \ w pętli i wywołania bezpośredniego.
Rozpatrzmy układ z macierzą losową, ale silnie diagonalnie dominującą:

n=100;
m=100;
A=2*eye(n,n) + (1/n)*rand(n,n);
S=rand(n,m);X=XX=zeros(n,m);
F=A*S;
tic;X=A\F;t1=toc
#a teraz w petli
tic;
for k=1:m,
  XX(:,k)=A\F(:,k);
endfor
t2=toc
t2/t1

Na moim komputerze czas obliczeń w pętli był ponad $ 50 $ razy dłuższy.

Sprawdźmy, czy wyniki są poprawne:

norm(X-XX,1)
norm(X-S,1)
norm(A*X-F,1)

Wszystkie błędy są małe.

Macierz odwrotna

W octave'ie funkcja inv(A) zwraca macierz odwrotną (o ile ona istnieje.)

Sprawdźmy, jak to działa:

A=[1 1;-1 1];
iA=inv(A);
A*iA
iA*A
norm(A*iA- eye(2),1)
norm(iA*A- eye(2),1)

Sprawdźmy, jak działa ta funkcja dla macierzy prawie osobliwej:

A=[1 1+2*eps;1 1];
iA=inv(A);
A*iA
iA*A
norm(A*iA- eye(2),1)
norm(iA*A- eye(2),1)

Macierz odwrotna jest poprawna. Teraz sprawdźmy, czy możemy rozwiązać układ równań liniowych z wykorzystaniem macierzy odwrotnej:

s=[1;3];
f=A*s;
x=iA*f;
x-s
norm(x-s,2)
norm(A*x-f,2)

Wynik jest znów nieprawidłowy.

Macierz odwrotną można też obliczyć jednym wywołaniem operatora
\:

iA=A\eye(size(A))

Przetestujmy ten sposób:

A=[1 1;-1 1];
Id=eye(size(A));
iA=A\Id
iA2=inv(A)
norm(iA*A,1)
norm(iA-iA2,1)

Norma różnicy obu macierzy odwrotnych jest identyczna, czyli oba sposoby są równoważne.

Współczynnik uwarunkowania macierzy

mo-NALcondHilb.jpg

mo-NALtestHres.jpg
mo-NALtestHer.jpg

Współczynnik uwarunkowania macierzy w danej normie definiujemy jako:

$$<br />
 \mathrm{cond}(A,\|\cdot\|)=\|A\|\|A^{-1}\|.<br />
$$

Jak wiadomo, por. [Kincaid:2006:AN], jeśli macierz jest źle uwarunkowana, tzn. jeśli ten współczynnik jest duży względem dokładności arytmetyki, to rozwiązywanie układu równań liniowych z macierzą $ A $ może być obarczone dużym błędem.

Znaną źle uwarunkowaną macierzą jest tzw. macierz Hilberta:

$$<br />
 H_N=\left(\frac{1}{i+j-1}\right)_{i,j=1}^N<br />
$$

Jest to macierz symetryczna i dodatnio określona.

W octave'ie funkcja H=hilb(N) tworzy macierz Hilberta tylko z zamienioną kolejnością wierszy i kolumn (ale permutowanie wierszy i kolumn nie zmienia uwarunkowania macierzy).

Sprawdźmy jej współczynnik uwarunkowania w normie drugiej:

M=25;
c=zeros(M,1);
for j=1:M,
 c(j)=cond(hilb(j));
endfor
semilogy(c,";wsp uwar macierzy Hilberta;")

Uwarunkowanie wzrasta wykładniczo - jak widać na wykresie w skali półlogarytmicznej, por. rysunek  [link].

Przetestujmy, jak działa operator octave'a \ zastosowany do rozwiązywania układu równań liniowych z macierzą Hilberta.

M=11;
er=err=zeros(M,1);
for j=1:M,
 H=hilb(j);
 s=ones(j,1);
 f=H*s;
 x=H\f;
 er(j)=norm(x-s,2);
 err(j)=norm(H*x-f,2);
endfor
plot(1:M,er,";blad;")
plot(1:M,err,";blad residualny;")

Błąd residualny, tzn. $ \|H*x-f\|_2 $ jest na poziomie błędu zaokrągleń, por. rysunek  [link], ale
błąd w normie drugiej gwałtownie rośnie dla $ N>9 $, por. rysunek  [link].
Dla wymiaru $ N\geq12 $ operator octave'a zwraca ostrzeżenie, że macierz jest numerycznie osobliwa, i że nie da się rozwiązać tego układu równań.
Octave wtedy rozwiązuje to równanie w sensie LZNK, por. Rozdział  [link].

Proszę zauważyć, że w praktyce nie znamy rozwiązania, więc możemy obliczyć tylko błąd residualny, tzn. normę $ H*x-f $. Norma błędu residualnego może być na poziomie błędu arytmetyki zmiennopozycyjnej, podczas gdy realny błąd może być o kilka rzędów, czy nawet kilkanaście rzędów wielkości większy.

Rozkłady macierzy

Podstawowymi algorytmami bezpośrednimi rozwiązywania układu ( [link])
są tzw. rozkłady macierzy: rozkład $ LU $, rozkład $ L^TL $ - o ile $ A $ jest symetryczna, dodatnio określona i rozkład $ QR $. Tutaj $ L $ oznacza macierz dolnotrójkątną, $ U, R $ - macierze górnotrójkątne, $ Q $ - macierz ortogonalną ($ Q^TQ=QQ^T=I $).

Rozkład LU

Podstawowym algorytmem jest tzw. eliminacja Gaussa, która dla dowolnej
macierzy nieosobliwej zwraca - w wersji z częściowym wyborem elementu głównego (czyli inaczej z częściowym osiowaniem) - następujący rozkład typu $ LU $:

$$<br />
 PA=LU<br />
$$

gdzie $ P $ to macierz permutacji, $ U $ to macierz nieosobliwa górnotrójkątna, $ L $ to macierz dolnotrójkątna z jedynkami na diagonali.

Znając ten rozkład możemy rozwiązać układ $ Ax=b $ w trzech krokach:

  1. Policz $ f=Pb $
  2. Rozwiąż układ z macierzą dolnotrójkątną $ Ly=f $ ($ PAx=LUx=Pb $)
  3. Rozwiąż układ z macierzą górnotrójkątną $ Ux=y $

Możemy też obliczyć wyznacznik $ A $ korzystająć z tego, że
$ \mathrm{det}( A)= \mathrm{det}(P) \mathrm{det}(U)=\mathrm{det}(P)\: \Pi_k u_{k k} $.
Wyznacznik $ P $ wynosi jeden albo minus jeden, w zależności od tego, czy permutacja jest parzysta, czy nieparzysta.
Jeśli macierz $ A $ jest symetryczna i dodatnio określona, tzn.

$$<br />
 A=A^T \qquad x^TA x>0 \quad x\not=0<br />
$$

to algorytm Choleskiego zwraca macierz dolnotrójkątną $ L $ taką, że

$$<br />
 A=L L^T.<br />
$$

Z kolei jeśli macierz $ n\times n $ $ A $ jest nieosobliwa, to istnieje rozkład $ QR $ tej macierzy, tzn. istnieją macierze tego samego wymiaru co $ A $: ortogonalna $ Q $ (tzn. $ Q^TQ=I $) oraz macierz nieosobliwa górnotrójkątna $ R $ takie, że

$$<br />
  A=QR.<br />
$$

Aby rozwiązać $ Ax=b $ wykorzystując ten rozkład należy rozwiązać równoważny układ:

$$<br />
  Rx=Q^Tb<br />
$$

z macierzą górnotrójkątną.

Do znalezienia rozkładu $ LU $ macierzy służy funkcja octave'a:

[L,U,P]=lu(A);

Argumentem jest macierz kwadratowa $ A $, a funkcja zwraca $ L,U,P $
- odpowiednie macierze rozkładu $ LU $.

Przetestujmy te rozkłady na kilku prostych macierzach. Na początek rozpatrzmy macierz symetryczną i znajdźmy jej rozkład $ LU $:

 A=[3 1;1 3]
[L,U,P]=lu(A)
 L*U-P*A
norm(L*U-P*A,1)

Wynik jest prawidłowy. Zastosujmy do rozwiązania układu równań liniowych z macierzą $ A $ i znanym rozwiązaniem:

 s=[1;-2];
 b=A*s;
 f=P*b;
 y=L\f;
 x=U\y
 x-s
 norm(x-s,1)
 xx=A\f
 xx-s
 norm(xx-s,1)

Do znalezienia rozkładu Choleskiego $ A=LL^T $ macierzy symetrycznej dodatnio określonej służy funkcja octave'a
R=chol(A), której argumentem jest macierz $ A $, a funkcja zwraca macierz górnotrójkątną $ R $ taką, że $ R^T*R=A $, czyli $ R $ to transponowana macierz $ L $ z naszego rozkładu Choleskiego.

Porównajmy rozkład Choleskiego tej samej macierzy co w poprzednim przypadku.

 R = chol (A)
 Er=R'*R-A
norm(Er,1)

Przetestujmy ten rozkład do znalezienia rozwiązania równania.

 y=R'\b;
 xc=R\y
 xc-s
 norm(xc-s,1)

Otrzymaliśmy poprawne wyniki.

Rozkład QR

Ostatni typ rozkładu macierzy, służący rozwiązaniu układu równań liniowych, to rozkład $ QR $. Do otrzymania takiego rozkładu w octave'ie służy funkcja:

[Q,R]=qr(A)

Argumentem jest macierz $ A $, a zwracanymi wartościami macierze rozkładu.

Przetestujmy tę funkcję na przykładzie tej samej macierzy $ A $, co w powyższym przykładzie:

 [Q,R]= qr (A)
 Er=Q*R-A
 norm(Er,1)
 norm(Q*Q'-eye(size(A)),1)

Zastosujmy ten rozkład do rozwiązania układu równań liniowych $ Ax=b $.

 f=Q'*b;
 xq=R\f
 xq-s
 norm(xq-s,1)

Metoda, jak widać, działa.

Przetestujmy rozkłady $ QR $ i $ LU $ dla macierzy nieosobliwej niesymetrycznej, lub symetrycznej ale nieokreślonej.
Rozpatrzmy macierz niesymetryczną nieosobliwą i jej rozkład $ LU $:

 A=[1 1;2 1]
[L,U,P]=lu(A)
 L*U-P*A
norm(L*U-P*A,1)

Przetestujmy rozkład $ QR $:

 [Q,R]= qr (A)
 Er=Q*R-A
 norm(Er,1)
 norm(Q*Q'-eye(size(A)),1)
 norm(Q'*Q-eye(size(A)),1)

Co się stanie jeśli zastosujemy metodę Choleskiego?

 clear R
 R = chol (A)

Powtórzymy to samo dla macierzy symetrycznej ale nieokreślonej:

Najpierw metoda Choleskiego pozwoli nam sprawdzić, czy macierz rzeczywiście nie jest określona:

 A=[-4 1; 1 5]
 R = chol (A)

Czyli macierz nie jest dodatnio określona.
Zauważmy, że algorytm Choleskiego jest najprostszą obliczeniowo metodą sprawdzenia, czy macierz jest dodatni określona.

Teraz znajdźmy rozkład $ LU $ macierzy $ A $:

 [L,U,P]=lu(A)
 L*U-P*A
norm(L*U-P*A,1)

a następnie rozkład $ QR $:

 [Q,R]= qr (A)
 Er=Q*R-A
 norm(Er,1)
 norm(Q*Q'-eye(size(A)),1)
 norm(Q'*Q-eye(size(A)),1)

Rozpatrzmy teraz przykład źle uwarunkowanej macierzy Hilberta dla większego $ N=5,15 $.
Najpierw przypadek $ N=5 $:

 N=5;
 H=hilb(N);
 R = chol (H);
 norm(R'*R-H,1)
 [L,U,P]=lu(H)
 norm(L*U-P*H,1)
 [Q,R]= qr (H);
 norm(Q*R-H,1)
 norm(Q*Q'-eye(size(H)),1)
 norm(Q'*Q-eye(size(H)),1)

Spróbujmy wykonać te same obliczenia dla $ N=20 $:

 N=20;
 H=hilb(N);
 R = chol (H);

Czy inne rozkłady można przeprowadzić?

 [L,U,P]=lu(H);
 norm(L*U-P*H,1)
 [Q,R]= qr (H);
 norm(Q*R-H,1)
 norm(Q*Q'-eye(size(H)),1)
 norm(Q'*Q-eye(size(H)),1)

Spróbujmy rozwiązać $ Ax=b $ za pomocą tych rozkładów:

 s= rand(N,1);
 f=H*s;
 y=L\(P*f);
 x=U\y;
 norm(H*x-f,1)
 norm(x-s,1)

Niestety wynik jest bardzo zły.
A co z wynikiem uzyskanym przy pomocy rozkładu $ QR $?

 xx=R\(Q'*f);
 norm(H*xx-f,1)
 norm(xx-s,1)

Uwarunkowanie macierzy tego układu jest bardzo duże. Żaden algorytm bezpośredni rozwiązywania układów równań liniowych, przy takiej dokładności obliczeń w arytmetyce zmiennopozycyjnej podwójnej precyzji, nie jest w stanie znaleźć rozwiązania.

Liniowe zadania najmniejszych kwadratów

Liniowe zadania najmniejszych kwadratów polega na znalezieniu $ x^*\in \mathbb{R}^n $, który minimalizuje

$$<br />
 \|Ax-b\|_2<br />
$$

dla danej macierzy $ A\in \mathbb{R}^{m,n} $ i wektora $ b\in \mathbb{R}^m $.

Zauważmy, że jeśli $ m=n $ i $ A $ jest nieosobliwa, to LZNK jest równoważne rozwiązaniu układu $ Ax=b $.

Funkcją octave'a rozwiązującą dowolne LZNK jest operator \.

Wywołanie

 x=A\b

powinno zwrócić rozwiązanie LZNK.

W przypadku niejednoznaczności LZNK rozwiązanie jest wybrane według odpowiedniego kryterium, o którym napiszemy poniżej w rozdziale  [link].

Przetestujmy ten operator na kilku bardzo prostych przykładach macierzy
wymiaru $ 3\times 2 $.

Pierwszy przykład z $ b $ będącym pierwszą kolumną:

 A=[1 1 1; 1 1 2]'
 b=[1 1 1]'
 x=A\b

Tutaj rozwiązanie jest jednoznaczne.

Weźmy macierz o niezerowym jądrze (przestrzeni zerowej), wtedy oczywiście dodanie wektora z jądra macierzy do $ x $ nie zmieni normy $ Ax-b $.

 A=[1 1 1; -1 -1 -1]'
 b=[1 1 1]'
 x=A\b

Jądro macierzy:

 y=null(A)
 norm(A*(x+y)-b,2)

Regularne LZNK

Jeśli macierz $ A $ jest kolumnami regularna (tzn. kolumny stanowią układ niezależny liniowo), to LZNK nazywamy regularnym i tylko wtedy ma ono jednoznaczne rozwiązanie.

Rozwiążmy najbardziej klasyczne zadanie najmniejszych kwadratów. Dla danych punktów $ (x_k,y_k)_{k=1}^m $ dopasuj krzywą postaci $ \sum_{j=1}^na_jf_j(x) $
dla znanych funkcji $ \{f_j\} $ tak, aby

$$<br />
  \min_{a_j}\sum_k |y_k-\sum_ja_jf_j(x_k)|^2<br />
$$

była minimalna.

To jest oczywiście LZNK z macierzą:

$$<br />
  A=(f_j(x_k))_{k,j}<br />
$$

i wektorem prawej strony $ b=(y_k) $.

mo-NALrLZNKdop1.jpg
Rozwiążmy to zadanie dla dowolnych punktów zadanych poprzez dwa pionowe wektory $ x=(x_k) $ i $ y=(y_k) $ oraz krzywych postaci $ a_1\cos(x)+a_2\sin(x) $:

 function a=dopasuj(x,y,N=length(x))
   A=[cos(x),sin(x)];
   a=A\y;
 end

Przetestujmy na punktach leżących na wykresie funkcji $ \sin $.

 x=linspace(0,2*pi,10)';
 y=sin(x);
 a=dopasuj(x,y)

Czyli $ a_1=0 $ i $ a_2=1 $ (tak, jak się spodziewaliśmy).

Teraz rozpatrzmy losowe punkty blisko wykresu $ \cos $, por. rysunek  [link]:

 y=cos(x)+0.5*(rand(10,1)-0.5);
 a=dopasuj(x,y)
 z=linspace(0,2*pi);
 plot(x,y,"+",z,a(1)*cos(z)+a(2)*sin(z),";...
 dopasowana krzywa;",z,cos(z),";cos;")

Teraz możemy np. dodać duże dodatkowe zaburzenie np. $ y_1 $ i zobaczyć jak zmieni się dopasowywana krzywa. Pozostawimy to jako zadanie.

Nieregularne LZNK i rozkład SVD

(#)

W programie Matematyki Obliczeniowej nie ma treści związanych z nieregularnym LZNK, niemniej my w tym skrypcie pokażemy jak można rozwiązywać takie LZNK w octave'ie.

W ogólności nie ma miejsca wtedy jednoznaczności rozwiązania - jak to już zauważyliśmy.

Do rozwiązywania nieregularnego LZNK służy rozkład macierzy SVD, czyli rozkład względem szczególnych wartości macierzy.

Rozkład SVD macierzy $ A $ (dowolnej) to

$$<br />
 A=U\Lambda V^T,<br />
$$

gdzie $ U,V $ to macierze ortogonalne, a $ \Lambda $ to macierz diagonalna
z nierosnącymi nieujemnymi wartościami na diagonali. Wszystkie te wartości z diagonali, które są nieujemne, oznaczmy przez $ \sigma_1,\ldots,\sigma_r $. Nazywamy je wartościami szczególnymi macierzy $ A $. Wyznaczone one są jednoznacznie, ale już sam rozkład SVD nie jest jednoznaczny.

Widać, że
\begin{eqnarray*}
\|Ax-b\|_2^2=\|U\Lambda V^Tx-b\|_2^2=\|\Lambda y- U^T b\|_2^2
\end{eqnarray*}
dla $ y=V^Tx $.
Jeśli więc znajdziemy rozwiązanie $ y $ LZNK dla macierzy diagonalnej $ \Lambda $ z prawą stroną
$ g=U^Tb $, to od razu $ x=Vy $.

Dalej otrzymujemy, że
\begin{eqnarray*}
\|\Lambda y- U^T b\|_2^2=\sum_{k=1}^r |\sigma_k y_k -g_k|^2 +
\sum_{k=r+1}^n |g_k|^2.
\end{eqnarray*}
Biorąc minimum po $ y_k $ widzimy, że rozwiązaniem tego LZNK jest dowolne $ y $ spełniające:

$$<br />
  y_k=\frac{g_k}{\sigma_k } \qquad k=1,\ldots,r.<br />
$$

$ y_{r+1},\ldots,y_n $ mogą przyjąć dowolne wartości.

Otrzymujemy więc, że rozwiązaniem wyjściowego LZNK jest podprzestrzeń:

$$<br />
    \left\{V\left(<br />
 \begin{array}{c}<br />
  w \\ z<br />
 \end{array}<br />
\right): z\in \mathbb{R}^{n-r}  \right\}=<br />
V\left(<br />
 \begin{array}{c}<br />
  w \\ 0<br />
 \end{array}<br />
\right)+\left\{V\left(<br />
 \begin{array}{c}<br />
  0\\ z<br />
 \end{array}<br />
\right): z\in \mathbb{R}^{n-r}  \right\}<br />
$$

dla $ w\in  \mathbb{R}^r $ takiego, że $ w_k=y_k=\frac{g_k}{\sigma_k } $.

Zbiorem rozwiązań jest podprzestrzeń afiniczna wymiaru $ n-r $.

Operator \ zwraca element zbioru rozwiązań o najmniejszej normie drugiej, czyli

$$<br />
V\left(<br />
 \begin{array}{c}<br />
  w \\ 0<br />
 \end{array}<br />
\right).<br />
$$

Do znalezienia rozkładu SVD macierzy służy funkcja SVD:

  [U, S, V] = svd(A)

gdzie $ U,V $ to macierze ortogonalne, $ S $ to macierz diagonalna z wartościami szczególnych macierzy $ A $ takie, że $ U*S*V^T=A $.
Można też wywołać tę funkcję jako:

  S = svd (A)

Wtedy funkcja zwróci tylko wartości szczególne macierzy w wektorze.

Przetestujmy tę funkcję:

 A=hilb(3);
 S=svd(A)
 [U, S, V] = svd(A)
 norm(U*S*V'-A,1)

Przetestujmy tę funkcję dla macierzy niekwadratowej:

 A=[1 1 1; 1 2 3];
 B=A';
 [U, S, V] = svd (A)
 norm(U*S*V'-A,1)
 [U, S, V] = svd (B)
 norm(U*S*V'-B,1)

Rozpatrzmy macierz kwadratową osobliwą:

 A=[1 1 ;2 2];
 [U, S, V] = svd (A)

Formalnie wynik jest fałszywy. Druga własność szczególna powinna być równa zero. Z uwagi na błędy zaokrągleń nie możemy się spodziewać tego, że otrzymamy zero, a wynik na poziomie $ 1e-16 $ możemy uznać za numeryczne zero.
Rząd macierzy nie jest numerycznie dobrze uwarunkowany, najmniejsze zaburzenie macierzy może zmienić rząd, a więc
badanie rzędu maciery w arytmetyce zmiennopozycyjnej jest właściwie niemożliwe.

Przetestujmy, jak działa operator \ oraz funkcja svd() zastosowane do rozwiązania nieregularnego LZNK. Na początek rozważymy LZNK z macierzą o maksymalnym rzędzie, ale o większej ilości kolumn niż wierszy:

 A=[1 1 1; 1 2 3]
 [U, S, V] = svd(A)
 s=[1;0;0];
 f=A*s;
 xx=A\f
 y=V(:,1:2)*(S(1:2,1:2)\(U'*f))
 norm(A*xx-f,2)
 norm(A*y-f,2)

Rozpatrzmy $ f $ nie będące w przestrzeni kolumn macierzy $ A $.

 f=rand(2,1);
 xx=A\f
 y=V(:,1:2)*(S(1:2,1:2)\(U'*f))
 norm(xx-y,2)

Zajmijmy się macierzą, która nie ma maksymalnego rzędu:

 A=[1 1 1 1; 2 2 2 2; 1 -1 2 3]'
 [U, S, V] = svd(A)
 s=[1;0;0];
 f=A*s;
 xx=A\f
 g=U'*f;
 x=V(:,1:3)*(S(1:3,1:3)\g(1:3))
 norm(A*xx-f,2)
 norm(A*x-f,2)

Znów rozpatrzmy $ f $ nie będące w przestrzeni kolumn macierzy $ A $:

 f=rand(4,1);
 xx=A\f
 g=U'*f;
 x=V(:,1:3)*(S(1:3,1:3)\g(1:3))
 norm(A*xx-f,2)
 norm(A*x-f,2)

Wynik jest inny niż tego oczekiwaliśmy. Przyjmijmy, że trzecia wartość szczególna wynosi zero:

 S(3,3)=0;
 x=V(:,1:2)*(S(1:2,1:2)\g(1:2))
 norm(A*xx-f,2)
 norm(A*y-f,2)
 norm(xx-x,2)

Otrzymaliśmy te same rozwiązania.

W praktyce pojawia się pytanie, jak ustalić epsilon takie, że
jeśli $ |\sigma_k|\leq \epsilon $, to przyjmiemy, że $ \sigma_k $ jest równa zero.
Jakiś próg musimy przyjąć, ale zawsze wrowadzi to pewien błąd do naszego przybliżonego rozwiązania LZNK.

Przy pomocy rozkładu SVD wyznacza się też współczynnik uwarunkowania macierzy. Porównajmy wyniki obliczania przybliżonego uwarunkowania macierzy otrzymane przez svd(A) i funkcję cond(A):

 H=hilb(15);
 S=svd(H);
 max(S)/min(S)
 cond(A)

Otrzymaliśmy ten sam wynik. Współczynnik uwarunkowania macierzy Hilberta dla $ N=15 $ jest równy $ 2.7878e+17 $.

W przypadku, gdy $ A $ jest symetryczna i dodatnio określona, to $ U=V $ w rozkładzie SVD tej macierzy, tzn. rozkład SVD jest równoważny
znalezieniu bazy złożonej z układu ortonormalnych wektorów własnych $ A $.

 H=hilb(4);
 [U, S, V] = svd(H);
 norm(U-V,1)

Zadanie własne

Numeryczne zadanie własne polega na znalezieniu niektórych lub wszystkich wartości oraz wektorów własnych macierzy $ A $ wymiaru $ N\times N $ rzeczywistej lub zespolonej. Na wykładzie z Matematyki Obliczeniowej omawia się kilka najprostszych metod dla symetrycznego zadania własnego, tzn. dla zadania własnego z macierzą rzeczywistą symetryczną.

Funkcją octave'a służącą do znajdowania wartości i wektorów własnych dowolnej macierzy jest
eig().

Najprostsze jej wywołanie to:

w=eig(A)

po którym funkcja powinna zwrócić wszystkie wartości własne w wektorze $ w $.

Jeśli potrzebujemy poznać wektory własne, to musimy wywołać tę funkcję z dwoma zwracanymi wartościami:

[V,D]=eig(A)

W tym przypadku $ D $ będzie macierzą diagonalną, na której diagonali będą wartości własne macierzy $ A $, natomiast $ V $ będzie macierzą ortogonalną, w której kolumnach będą wektory własne dla odpowiednich wartości na diagonali $ D $, czyli

$$<br />
  AV=VD.<br />
$$

Sprawdźmy na losowej macierzy diagonalnej jak działa funkcja eig(),
najpierw przy prostszym wywołaniu:

A=rand(5,5);
A=A+A'; # symetryzujemy macierz
w=eig(A)
det(A-w(1)*eye(size(A)))

Wyznacznik $ A-w(1)*I $ jest na poziomie $ 10^{-14} $, czyli de facto jest równy zero.

Teraz sprawdźmy drugą formę wywołania funkcji eig():

[V,D]=eig(A)
norm(V*D-A*V,1)

Sprawdźmy, jak działa funkcja dla macierzy większego wymiaru:

N=10;
for k=1:3,
  A=rand(N,N);
  A=A+A'; # symetryzujemy macierz
  [V,D]=eig(A);
  norm(V*D-A*V,1)
  N*=10
endfor

Wniosek: funkcja działa dobrze nawet dla większych $ N $.

Przetestujmy ją dla macierzy ogromnego wymiaru o dużej ilości zer, np. macierzy trójdiagonalnej. Wykorzystamy też funkcje sparse() służące utworzeniu macierzy w formacie rzadkim:

N=1000;
x=ones(N,1);
A=sparse(diag(x))-sparse(diag(x(1:N-1),1));
A=A+A';
tic;[V,D]=eig(A);toc
norm(V*D-A*V,1)

Sprawdźmy, jak działa funkcja dla macierzy niesymetrycznej:

A=[2,-1;1,2]
w-eig(A)
[V,D]=eig(A)
norm(V*D-A*V,1)

Istnieje też funkcja eigs(), która oblicza tylko kilka wartości własnych macierzy, np. największych co do modułu wartości.

Przetestujmy funkcję dla losowej macierzy wymiaru $ 10\times 10 $:

A=rand(10,10);
A=A+A'; # symetryzujemy macierz
wa=eig(A)
w=eigs(A)

Funkcja domyślnie znajduje sześć największych co do modułu wartości własnych, ale może znaleźć mniej lub więcej:

w3=eigs(A,3)
w9=eigs(A,8)

lub bazując na innym kryterium, np. najmniejsze wartości własne co do modułu:

wm=eigs(A,3,'sm')
wd=eigs(A,3)

Możliwe opcje to:

  • 'lm'
    największe wartości co do modułu (default).

  • 'sm'
    najmniejsze wartości co do modułu

  • 'la'
    największe wartości (tylko o ile $ A $ jest symetryczna)

  • 'sa'
    najmniejsze wartości (tylko o ile $ A $ jest symetryczna)

  • 'be'
    wartości z końców widma macierzy, z jedną więcej z górnej części widma, o ile chcemy znaleźć nieparzystą liczbę wartości własnych (tylko o ile $ A $ jest symetryczna)

  • 'lr'
    największe części rzeczywistych wartości własnych (tylko o ile $ A $ jest niesymetryczna lub zespolona).

  • 'sr'
    najmniejsze części rzeczywistych wartości własnych (tylko o ile $ A $ jest niesymetryczna lub zespolona).

  • 'li'
    największe wartości części urojonych wartości własnych (tylko o ile $ A $ jest niesymetryczna lub zespolona).

  • 'si'
    najmniejsze wartości części urojonych wartości własnych (tylko o ile $ A $ jest niesymetryczna lub zespolona).

Spróbujmy policzyć dwie największe i dwie najmniejsze co do wartości własne macierzy trójdiagonalnej przy pomocy tej funkcji:

N=1000;
x=ones(N,1);
A=sparse(diag(x))-sparse(diag(x(1:N-1),1));
A=A+A';
tic;w=eigs(A,4,'be');toc

Niestety metoda zawiodła, a funkcja eig() - choć wolna, zwróciła wszystkie wartości i wektory tej macierzy.

Zaimplementujmy najprostszą wersję metody potęgowej:

Niech $ \hat{x}_0=x_0 $ będzie wektorem o normie drugiej równej jeden. Wtedy
ciąg iteracji metody potęgowej jest zdefiniowany następująco:
\begin{equation}(#)
\begin{array}{lclr}
\hat{x}_k&=&Ax_{k-1} &\qquad \mathrm{iteracja}\\
x_k&=&\frac{ \hat{x}_k}{ \|\hat{x}_k\|_2} & \qquad \mathrm{normowanie} \\
r_k&=& x_k^TAx_k \qquad &\mathrm{iloraz      Rayleigha}
\end{array} \qquad k>0
\end{equation}

W octavie kod może wyglądać następująco:

x=rand(N,1); #losujemy wektor
x=x/norm(x,1); #normowanie x0
y=A*x;
r=x'*y; rp=r+2*TOL;
it=0;
while(abs(r-rp)>TOL)
 x=y/norm(y,2);
 y=A*x;
 rp=r;
 r=x'*y;
 it++;
endwhile
x=y/norm(y,2)
r

Za warunek stopu przyjęliśmy $ |r_k-r_{k-1}|<TOL $ zadanej tolerancji.

Przetestujmy ten algorytm dla macierzy symetrycznej A=[4,1;1,4], która ma wartości własne $ 5,3 $.

Dla TOL=1e-12 otrzymaliśmy, że po $ 24 $ iteracjach iloraz Rayleigha wynosi
$ r =  4.99999999999963 $, czyli błąd $ 3.71e-13 $.
Z kolei przybliżenie wektora własnego to

x =

   0.707106628756140
   0.707106933616922

Możemy policzyć $ \|Ax-5*x\|_2 $ w octave'ie, czyli norm(A*x-5*x,2), która wynosi
$ 4.31e-07 $.

Po prostych rachunkach otrzymujemy, że wektorem własnym dla wartości własnej $ 5 $ jest wektor postaci $ x=(x_1,x_1)^T\not=0 $.

Sprawdźmy, dla naszego wektora:

octave:111> abs(x(1)-x(2))
ans =  3.04860781730198e-07
octave:112>

Oczywiście ilość iteracji i błędy zależą też od wektora $ x_0 $.

Powtórzmy obliczenia drukując na ekranie błąd $ |5-r_k| $ i $ |(1,-1)^Tx_k| $ dla $ x_k $ iteracji metody potęgowej i $ r_k $ ilorazu Rayleigha:

A =

   4   1
   1   4

[it]|  Rayleigh  |   |x-y|   |
-----------------------------
[ 1] | 2.5008722 | 0.0241133 |
[ 2] | 4.9995814 | 0.0204587 |
[ 3] | 4.9998493 | 0.0122760 |
[ 4] | 4.9999457 | 0.0073658 |
[ 5] | 4.9999805 | 0.0044195 |
[ 6] | 4.9999930 | 0.0026517 |
[ 7] | 4.9999975 | 0.0015910 |
[ 8] | 4.9999991 | 0.0009546 |
[ 9] | 4.9999997 | 0.0005728 |
[10] | 4.9999999 | 0.0003437 |
[11] | 5.0000000 | 0.0002062 |
[12] | 5.0000000 | 0.0001237 |
[13] | 5.0000000 | 0.0000742 |
[14] | 5.0000000 | 0.0000445 |
[15] | 5.0000000 | 0.0000267 |
[16] | 5.0000000 | 0.0000160 |
[17] | 5.0000000 | 0.0000096 |
[18] | 5.0000000 | 0.0000058 |
[19] | 5.0000000 | 0.0000035 |
[20] | 5.0000000 | 0.0000021 |
[21] | 5.0000000 | 0.0000012 |
-----------------------------
x =

   7.07106556762792e-01
   7.07107005610232e-01

Wartosc wlasna:   4.999999999999440

Powtarzając eksperyment dla macierzy $ A-2*I $, czyli macierzy o wartościach własnych $ 3 $ i $ 1 $, otrzymujemy po $ 14 $ iteracjach iloraz Rayleigha taki, że błąd $ |3-r| $ wynosi $ 3.99e-14 $ przy tym samym warunku stopu:

p =  2
A =

   2   1
   1   2

[it]| Rayleigh+p |   |x-y|   |
-----------------------------
[ 1] | 3.5005623 | 0.0335338 |
[ 2] | 4.9997501 | 0.0158070 |
[ 3] | 4.9999722 | 0.0052693 |
[ 4] | 4.9999969 | 0.0017564 |
[ 5] | 4.9999997 | 0.0005855 |
[ 6] | 5.0000000 | 0.0001952 |
[ 7] | 5.0000000 | 0.0000651 |
[ 8] | 5.0000000 | 0.0000217 |
[ 9] | 5.0000000 | 0.0000072 |
[10] | 5.0000000 | 0.0000024 |
[11] | 5.0000000 | 0.0000008 |
-----------------------------
x =

   7.07106736568327e-01
   7.07106825804765e-01

Wartosc wlasna:   4.999999999999928

Widać już z tych dwóch eksperymentów, że ilość iteracji zależy od stosunku $ \lambda_2/\lambda_1 $ ( gdzie $ \lambda_1 $ to wartość własna o maksymalnym module, $ \lambda_2 $ to druga z kolei wartość własna ze względu na moduł) i że szybkość zbieżności ilorazu Rayleigha jest dużo większa niż
wektorów iteracji do odpowiedniego wektora własnego.
Potwierdza to oszacowania teoretyczne, w których można pokazać, że

$$<br />
  v_1=x_k+O((\lambda_2/\lambda_1)^k), \qquad  \lambda_1=r_k+O((\lambda_2/\lambda_1)^{2k}).<br />
$$

Widać, że iloraz Rayleigha zbiega do wartości własnej dużo szybciej w tym przypadku.

Jeśli $ (\lambda_k,v_k) $ to para własna dla macierzy $ A $, to $ ((\lambda_k-p)^{-1},v_k) $ jest parą własną dla macierzy $ (A-pI)^{-1} $. Ta obserwacja pozwala nam skonstruować tzw. odwrotną metodę potęgową, która jest zdefiniowana jako metoda potęgowa zastosowana do macierzy $ (A-pI)^{-1} $ dla $ p $ parametru będącego dobrym przybliżeniem $ \lambda_k $ dowolnej wartości własnej $ A $. Warunkiem zbieżności iteracji odwrotnej metody potęgowej do $ v_k $ wektora własnego dla $ \lambda_k $, przy założeniu, że wektor startowy nie jest ortogonalny do $ v_k $ jest to, aby

$$<br />
|\lambda_k-p|>|\lambda_j -p| \qquad \forall j\not= k.<br />
$$

Szybkość zbieżności wektorów iteracji dla odwrotnej metody potęgowej zależy wtedy od
$ \max_{j\not =k}\frac{|\lambda_k-p|}{|\lambda_j-p|} $.

W octave'ie łatwo zaimplementować odwrotną metodę potęgową: wystarczy mnożenie przez macierz $ A $ w metodzie potęgowej zastąpić przez mnożenie wektora przez macierz odwrotną, które jest równoważne rozwiązywaniu układu równań. W octave kod odwrotnej metody potęgowej może wyglądać następująco:

A=A-p*eye(size(A)); # tworzymy macierz A-p*I

x=rand(N,1); #losujemy wektor
x=x/norm(x,1); #normowanie x0
y=A\x;
r=x'*y;
rp=r+2*TOL;
it=1;

while(abs(r-rp)>TOL)
 x=y/norm(y,2);
 y=A\x;
 rp=r;
 r=x'*y;
 it++;
endwhile

x=y/norm(y,2) #przyblizenie wektora wlasnego A
1/r+p #przyblizenie wartosci wlasnej A

Zobaczmy jak szybko działa metoda dla naszej testowej macierzy $ A $ i $ p=5.01 $:

p =  5.0100
A =

  -1.0100   1.0000
   1.0000  -1.0100

[it]|   1/r+p    |   |x-y|   |
-----------------------------
[ 1] | 4.9900001 | 0.0273613 |
[ 2] | 5.0000000 | 0.0001925 |
[ 3] | 5.0000000 | 0.0000010 |
[ 4] | 5.0000000 | 0.0000000 |
-----------------------------
x =

  -7.07106781186489e-01
  -7.07106781186606e-01

Wartosc wlasna:   5.000000000000000

A teraz dla $ p=2.8 $ - tu wektor własny ma postać $ a*(1,-1)^T $:

p =  2.8000
A =

   1.2000   1.0000
   1.0000   1.2000

[it]|   1/r+p    |   |x+y|   |
-----------------------------
[ 1] | 7.1576796 | 1.0000000 |
[ 2] | 3.9194920 | 1.3442338 |
[ 3] | 3.0139676 | 0.3789639 |
[ 4] | 3.0001162 | 0.0357476 |
[ 5] | 3.0000010 | 0.0032508 |
[ 6] | 3.0000000 | 0.0002955 |
[ 7] | 3.0000000 | 0.0000269 |
[ 8] | 3.0000000 | 0.0000024 |
[ 9] | 3.0000000 | 0.0000002 |
-----------------------------
x =

   7.07106782104050e-01
  -7.07106780269045e-01

Wartosc wlasna:   3.000000000000000

Przesuwanie widma macierzy poprzez dodanie do niej $ p*I $ jest bardzo ważne w różnych uogólnieniach metody potęgowej takich jak np. metoda $ QR $ z przesunięciami, por. np. [Kincaid:2006:AN], [Trefethen:1997:NLA], [Jankowska:1982:MN].