Skrypty i funkcje

# W pliku programik.m znajduje się trywialny program testowy Octave
1/7
2/7
3/7

> programik

ans =  0.14286
ans =  0.28571
ans =  0.42857

> help programik

octave:6> help programik
programik is the file ./programik.m

W pliku programik.m znajduje się trywialny program testowy Octave
...

Definiowanie funkcji

W poprzedniej części kursu wyjaśniliśmy, jak łatwo w Octave można utworzyć własne programy (tzw. skrypty). Funkcje użytkownika najłatwiej można zdefiniować właśnie w formie skryptów.

Załóżmy, że naszym celem jest zdefiniowanie funkcji o nazwie sinusik, która dla argumentu x obliczać będzie wartość wyrażenia x - x3/6 + x5/120. Aby zdefiniować taką funkcję:

1. Tworzymy (w katalogu bieżącym) plik o nazwie sinusik.m
2. W pliku tym wpisujemy:

   function y = sinusik(x)
      y = x - x**3/6 + x**5/120;
   endfunction
   

Uwaga! Powyższa definicja nie jest optymalna i wkrótce ją poprawimy.

Funkcje tę wywołujemy w Octave w naturalny sposób:

>> sinusik (0.1)
ans =  0.099833

Ogólnie, podczas definiowania funkcji w pliku obowiązują następujące zasady:

• Plik musi mieć taką samą nazwę jak funkcja i rozszerzenie m.
• Plik musi być widoczny dla programu (np. znajdować się w katalogu bieżącym programu).
• Definicja funkcji musi zawierać się między poleceniami function i endfunction.
• Funkcja może przyjmować wiele argumentów, np.

  function y = f(a, b, c) 
  

• Funkcja może zwracać wiele wartości (będących składowymi wektora), np:

  function [y, eps, xi] = f(a, b, c) 
  
  

• Funkcja może składać się z wielu instrukcji, które warto kończyć średnikiem:

     y = x - x**3/6;
     y = y + x**5/120;
  

  Niezakończenie jakiejś instrukcji średnikiem spowoduje wyświetlenie wyniku  wykonywanych w niej obliczeń – ta cecha może być  przydatna np. podczas szukania błędów w skrypcie, ale rzadko kiedy ma zastosowanie w gotowej, przetestowanej funkcji.

Odwzorowania

W wielu językach programowania, w tym w Octave, istnieje możliwość wywoływania funkcji (np. matematycznych) z argumentami wektorowymi lub macierzowymi. Takie funkcje nazywa się odwzorowaniami (ang. maps lub mapping functions). Jedną z takich funkcji w Octave jest cos, czyli funkcja wyznaczająca cosinus kąta:

> cos([0, 0.1, 0.2])
ans =
   1.00000   0.99500   0.98007

Jak widać, jeżeli argumentem funkcji cos jest wektor, wynikiem jest wektor utworzony z wartości tej funkcji dla wszystkich elementów wektora [0, 0.1, 0.2], czyli

cos([0, 0.1, 0.2]) = [cos(0), cos(0.1), cos(0.2)].

Ogólnie mówimy, że f jest odwzorowaniem, jeśli dla dowolnej macierzy v o rozmiarze n na m wyrażenie f(v) też jest macierzą n na m, przy czym jeżeli podstawimy w = f(v), to spełniona jest zależność w(k,l) = f(v(k,l)). Innymi słowy, w k-tym wierszu i l-tej kolumnie macierzy f(v) znajduje się wartość funkcji f dla argumentu v(k,l) leżącego w k-tym wierszu i l-tej kolumnie macierzy v.

Informacja o tym, czy dana funkcja Octave jest odwzorowaniem, znajduje się w dokumentacji tej funkcji (help), w której odwzorowania określane są jako mapper function lub mapping function, zwykle już w pierwszym wierszu opisu.

Odwzorowania są w Octave bardzo ważną klasą funkcji, gdyż znacznie przyspieszają rozwiązywanie wielu podstawowych problemów numerycznych, takich jak znajdowanie pierwiastków równań z wieloma niewiadomymi, rozwiązywanie równań różniczkowych czy obliczanie pól powierzchni (całek), a także tworzenie wykresów funkcji.

Operatory „z kropką”

Jak można się łatwo przekonać, nasza funkcja sin5 nie jest odwzorowaniem i nie można jej wywoływać z argumentami innymi niż skalary, co znacznie ogranicza jej użyteczność. Można temu jednak łatwo zaradzić. Aby funkcja sinusik stała się odwzorowaniem, wystarczy poprzedzić kropką wszystkie występujące w jej definicji operatory potęgowania:

function y = sinusik(x)
   y = x - (x.**3)/6 + (x.**5)/120;
endfunction

Ogólnie, kropką można poprzedzić dowolny operator “iloczynowy”, czyli operator mnożenia (.*, dzielenia (./ lub .\) i potęgowania (.** lub .^). To, kiedy stosować tę kropkę, zależy od następujących czynników:

• • Nie ma większego sensu stosować kropki przy mnożeniu lub dzieleniu przez liczby (ogólnie: skalary). Np. x/2 jest zawsze równoważne wyrażeniu x./2.
  • Jeżeli nie znasz rachunku macierzowego, to zawsze poprzedzaj operatory mnożenia, dzielenia i potęgowania kropką (z wyjątkiem mnożenie lub dzielenia przez liczbę).
  • Jeżeli znasz rachunek macierzowy, to pamiętaj, że mnożenie, dzielenie i potęgowanie “bez kropki” to operacje macierzowe. Jeżeli ani x, ani y nie jest skalarem, to:
    • w wyrażeniu x * y macierz x musi mieć tyle kolumn, ile wierszy ma macierz y.
    • wyrażenie x/y jest równoważne wyrażeniu ((y'-1)x')', ale nie wyznacza się wartości y'-1.
    • wyrażenie y\x ma wartość (y-1)*x, która jest wyznaczana szybkim algorytmem niewymagającym wyznaczania y-1.
    • w wyrażeniach x.*y oraz x./y macierze x i y muszą mieć ten sam rozmiar. Operacje z kropką wykonywane są “element po elemencie”. Np.

      >> m = [1 -1; -1 1]
      m =
         1  -1
        -1   1
      >> v = [1 2; 3 4]
      v =
         1   2
         3   4
      >> m .* v # mnożenie element po elemencie
      ans =
         1  -2
        -3   4
      >> m * v # mnożenie macierzowe
      ans =
        -2  -2
         2   2
      

  • Wyrażenie x**y (lub x^y) ma bardzo nieoczekiwaną interpretację, jeżeli x jest skalarem, a y czymś innym niż skalar (np. macierzą kwadratową). Dopóki nie zrozumiesz, jak to działa, używaj wyłącznie “potęgowania z kropką”: .** lub .^. Na przykład, aby wygenerować ciąg dziesięciu kolejnych potęg liczby 2, można użyć następującego wyrażenia:

    >> 2.^(1:10)
    ans =
         2     4     8    16    32    64   128   256   512  1024

  • Nie używaj operatorów  .+ i .–.

• Pamiętaj, że kropka może też oznaczać część ułamkową liczby. Jeżeli zachodzi konflikt tych znaczeń, Octave przyjmie, że kropka łączy się z operatorem. Np. wyrażenie

  >> 2./m	
  

  zostanie zinterpretowane następująco:

  2 ./ m
  

Na zakończenie warto dodać, że jeżeli w pliku zawierającym definicję funkcji umieścimy komentarz (najlepiej przed jej definicją), to będzie on wyświetlany przez Octave po wydaniu polecenia help nazwa_funkcji (np. help sinusik). Stąd też ostateczna zawartość pliku sinusik.m może wyglądać następująco:

# Funkcja sinusik przybliża wartość funkcji sin(x) za pomocą wielomianu stopnia 5.
#
# Przybliżenie to jest całkiem dobre dla małych wartości x, 
# np. dla |x| < 1.  
#

function y = sinusik(x)
   y = x - x.**3/6 + x.**5/120; 
endfunction

Pliki z funkcjami a pliki ze skryptami

Zarówno funkcje, jak i skrypty zapisuje się w m-plikach, czyli plikach o rozszerzeniu .m.  Octave odróżnia jedne od drugich na podstawie zawartości pierwszej instrukcji w pliku. Jeśli jest nią definicja funkcji (ze słowem kluczowym function), cały plik uznawany jest za definicję funkcji; w przeciwnym wypadku – za skrypt. Skrypty wykonywane są w całości, natomiast po kodzie definiującym funkcje można w pliku umieścić dowolny inny kod, np. definicje funkcji pomocniczych. Ten dodatkowy kod nie będzie widoczny z zewnątrz, ale może być wykorzystywany w kodzie funkcji.

Quiz z funkcji

1. Dlaczego definicje funkcji zwykle umieszcza się w osobnych plikach?
2. Jaka jest rola słów kluczowych function i endfunction?
3. Czy funkcja może zwracać wiele wartości?
4. Dlaczego zapis instrukcji funkcji z reguły kończy się średnikiem?
5. Co to są odwzorowania?
6. Jak sprawdzić, czy dana funkcja Octave jest odwzorowaniem?
7. Jak definiuje się własne odwzorowania?
8. Jaka jest różnica między x*y i x.*y?
9. Jaka jest różnica między x/y i x./y?
10. Jaka jest różnica między x**2 i x.**2?
11. Jaka jest różnica między x+y i x.+y?
12. Skąd wiadomo, kiedy m-plik zawiera skrypt, a kiedy definicję funkcji?

Zadania z definiowania funkcji

1. Zdefiniuj funkcję f(x) zadaną wzorem: f(x) = x \sqrt {|1 - x|} \sin^2(x)  + x\cos(x) .
   Pierwiastek kwadratowy można w Octave wyznaczyć funkcją sqrt, a wartość bezwzględną funkcją abs.
2. Wykonaj wykres f(x) dla 0 ≤ x ≤ 2.5. W tym celu możesz użyć następującego ciągu poleceń:

   >> x = 0 : 0.001 : 2.5; # zakres z krokiem 0.001
   >> y = f(x);
   >> plot (x,y)	
   

   Pierwsze z nich definiuje wektor “iksów”, drugie – wektor odpowiadających mu wartości (“igreków”), a trzecie powoduje wygenerowanie takiego wykresu jak ten:
   

3. Z wykresu wynika, że równanie f(x) = 0 ma pierwiastek z równy w przybliżeniu 2.2. Wyznacz jego wartość za pomocą polecenia fsolve:

   >> z = fsolve("f", 2.2)
   

4. Oblicz pole powierzchni pomiędzy wykresem funkcji f(x) dla 0 ≤ x ≤ z a osią rzędnych (tj. y = 0) (na poniższym rysunku zostało ono zaznaczone kolorem zielonym):
   
   Możesz posłużyć się następującym poleceniem:

   >> pole = quad("f", 0, z)
   

5. Zapoznaj się z dokumentacją funkcji quad i fsolve (przy okazji zapamiętaj te funkcje na przyszłe zajęcia). Jakie jest znaczenie używanych przy nich parametrów?