Definiowanie funkcji

W poprzedniej części kursu wyjaśniłem, 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 sin5, 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 sin5.m
  2. W pliku tym wpisujemy: 
    function y = sin5(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:

>> sin5 (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 sin5 stała się odwzorowaniem, wystarczy poprzedzić kropką wszystkie występujące w jej definicji operatory potęgowania:

function y = sin5(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 sin5). Stąd też ostateczna zawartość pliku sin5.m może wyglądać następująco:

# Funkcja sin5 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 = sin5(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

  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.+2?
  12. Skąd wiadomo, kiedy m-plik zawiera skrypt, a kiedy definicję funkcji?

Zadania

  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:
    wykres funkcji

  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):
    wykres funkcji
    Możesz posłużyć się następującym poleceniem:

    >> pole = quad("f", 0, z)
  5. Zapoznaj się z dokumentacją funkcji quad i fsolve. Jakie jest znaczenie ich parametrów? Jak można oszacować dokładność rozwiązań?