Modelování 1: Programování funkcí a skriptů

V tomto článku bych chtěl ukázat, jak v MATLABu vytvářet vlastní funkce a vysvětlit základní koncepty nezbytné k další práci s programem Matlab v předmětu “Modelování”.

Úvodem o funkcích

Koncept celého systému Matlab je postaven na funkcích. Funkce je několik řádek kódu a jejím úkolem je uskutečnit konkrétní operaci, například provést nějaký výpočet. Síla MATLABu je v tom, že již v základní verzi existuje velké množství funkcí, které umožňují řešit širokou oblast problémů. K základní verzi MATLABU je možné doinstalovat také řadu rozšíření – toolboxů, což není nic jiného než knihovny dalších funkcí. V našem předmětu budeme hlavně používat funkce, jejichž úkolem je například integrace soustavy obyčejných diferenciálních rovnic (ODE) nebo hledání řešení soustavy algebraických rovnic (AE). Tyto funkce jsou již obsaženy v základní sestavě Matlabu a nebudeme je muset programovat. Pro každou úlohu, kterou budeme řešit, si ale budeme muset vytvořit naše vlastní funkce popisující konkrétní problém, jako je například matematický model reaktoru.

Aby systém Matlab mohl funkci provést, tak ji musí najít. Systém hledá soubory s funkcemi na disku jen v předem v určených složkách a také v aktuální / pracovní složce. Pro přehlednost doporučuji pro každou úlohu vytvořit vlastní složku, nastavit pracovní složku jako tuto složku a všechny vlastní funkce potřebné pro konkrétní úlohu mít uložené v této složce.

Ukažme si to na příkladech. Napišme v příkazovém řádku postupně

>> which sin
>> which ode45
>> which mojefunkce
>> pwd

Příkaz which ukáže, kde se nachází daná funkce. Kód funkce se nachází v souboru, který musí mít stejné jméno jako je jméno funkce s koncovkou “.m”. Vidíme, že funkce sin je interní (build-in) funkce, což znamená, že její kód neuvidíme. Naopak funkce ode45 je napsána přímo v programovacím jazyce Matlabu a její kód bychom si, pokud by nás to zajímalo, mohli přečíst v uvedeném souboru. Třetí příkaz pravděpodobně vrátí informaci, že funkce mojefunkce nebyla nalezena. Poslední příkaz, pwd ukáže aktuální nastavení pracovní složky, tj. ve které složce na disku Matlab právě pracuje.

Změna aktuální / pracovní složky

Vytvořte si někde na disku složku, kam budete ukládat dnes vytvořené funkce: například Z:\matlab2013\week01\example01. Potom nastavte pracovní složku do této vytvořené složky. To se může udělat příkazem cd

>> cd 'Z:\matlab2013\week01\example01'

ale pohodlněji to lze udělat v panelu “Current Folder”, který by se měl nacházet v levé části grafického prostředí. Klikáním přepněte pracovní složku a příkazem pwd se přesvědčte, že se nacházíte ve správné složce. V panelu “Current Folder” se zároveň ukazují soubory (funkce) přítomné v pracovní složce. V tuto chvíli bude složka pravděpodobně prázdná.

Jak fungují funkce v Matlabu?

Jistě jste se již v Matlabu nebo jiném programovacím jazyce setkali s pojmem proměnná. Proměnná je pojmenovaný objekt, ke kterému je přiřazena nějaká hodnota. Pokud vytvoříme matici a uložíme ji do proměnné a, potom s maticí můžeme pracovat, protože ji Matlab bude znát, např.

>> a = [1, 2; 3, 4]
>> a = a + 4

Předpokládejme, že teď zavoláme nějakou funkci a program začne vykonávat příkazy v této funkci zapsané. Po celou dobu, kdy se bude nacházet uvnitř této funkce, ale proměnnou a ani jinou proměnnou přiřazenou z příkazovém řádku “nevidí”, takže ji nelze číst, přepsat jinou hodnotou nebo zrušit. Uvnitř funkce lze vytvořit proměnnou, která bude mít stejné jméno, tj. a, ale bude se jednat o naprosto jiný objekt. Toto a můžeme v rámci funkce používat, ale po ukončení funkce a návratu kontroly do příkazového řádku bude systém “vidět” opět původní proměnnou a. Tomuto chování budeme říkat, že každá funkce má vlastní jmenný prostor (name space). To znamená, že hodnoty proměnných do funkce nebo z funkce lze předávat pouze přes rozhranní (interface) funkce. Každá funkce má rozhranní nadefinované na svém prvním řádku. Například první řádek souboru “mojefunkce.m” bude

function [out1, out2] = mojefunkce(in1, in2, in3)

Vidíme, že po klíčovém slově function je v hranatých [] závorkách uveden seznam výstupních proměnnych, dále je jméno funkce konzistentní se jménem souboru a za ním v kulatých závorkách () seznam vstupních proměnných. Názvy proměnných odpovídají volané části rozhranní, tj. názvům známým uvnitř funkce.

Na příkazovém řádku nebo v jiné nadřazené funkci, tj. na volací části rozhranní můžeme funkci mojefunkce zavolat takto

[a,b] = mojefunkce(d,sin(x),2)

Všimněte si, že proměnné ve volacím a volaném prostoru nemusí mít stejná jména, pro jejich propojení je důležité pouze jejich pořadí a počet.

[Seznam proměnných ve volací funkci může mít menší počet prvků než seznam ve volané funkci: Pokud chybí vstupní proměnné o kterých autor funkce rozhodl, že jsou nepovinné, může tyto nahradit přednastavenou hodnotou. Chybí-li výstupní proměnné, pak se předpokládá, že autor volací funkce nemá pro tyto žádné využití. Při volání námi vytvářených funkcí budeme zatím používat stejný počet proměnných na obou stranách rozhraní.]

Předávání proměnných mezi prostory probíhá v Matlabu voláním hodnotou, kdy aktuální hodnota proměnné ve volací funkci je předána do proměnné ve volané funkci. Vidíme, že na volací části múžeme pro vstupní proměnné použít libovolný výraz, včetně výsledku volání jiné funkce (sin(x)) nebo konstantu (2). Uvnitř prostoru funkce mojefunkce bude systém znát proměnné:

• in1 se stejnou hodnotou, jako měla proměnná d ve volacím prostoru,
• in2 s hodnotou odpovídající výsledku volání funkce sin ve volacím prostoru a
• in3 s hodnotou 2.

V průběhu operací ve volané funkci, múžeme vstupní proměnné in1/in2/in3 měnit, aniž by to způsobilo chybu nebo ovlivnilo hodnoty proměnných ve volacím prostoru, ale dle zásad slušného programování, doporučuji se vstupními proměnnými uvnitř funkce zacházet jako by byly jen pro čtení. Výsledky operací je nutné uložit do výstupních proměnných out1 a out2. Po opuštění volané funkce a návratu do volacího prostoru budou výsledné hodnoty uloženy v proměnných a a b a k dispozici pro další práci s nimi.

Po tomto únavném popisu konečně můžeme napsat první funkci.

Vytváříme první vlastní funkci

Udělejme si funkci, která po zadání teploty vypočítá tenzi par hexanu z Antoineho rovnice . Vstupním parametrem bude teplota, výsledkem funkce bude tenze par a funkce se bude jmenovat phexan. Napište na příkazovém řádku

>> edit phexan

nebo v menu vyberte “File > New > Function“. Měl by se otevřít editor, do kterého napíšeme vlastní kód funkce. Soubor potom nezapomeňme uložit tak, aby se nacházel v pracovním adresáři a jmenoval se “phexan.m”.) Celý program bude vypadat takto:

function [p0] = phexan(t)
%PHEXAN Calculate hexane vapour pressure 
%   Hexane vapour pressure calculation using Antoine equation.
%
%   P0 = PHEXAN(T)
%
%   T is temperature in Celsius and P0 is vapour pressure in kPa.
%
A = 13.8216;           % Antoine equation coefficients for hexane
B = 2697.55;
C = -48.78;
tK = t+273.15;         % convert temperature to Kelvins
p0 = exp(A-B/(tK+C));  % vapour pressure in kPa
end

Popišme si nyní jednotlivé řádky.

• Na prvním řádku musí být klíčové slovo function, které určuje, že se jedná o funkci (a ne o skript) a zároveň definuje rozhranní funkce, tak jak bylo vysvětleno dříve. Naše funkce má jednu vstupní proměnnou t a jednu výstupní proměnnou p0.
• Dále následuje blok řádků začínajících znakem procenta %. Obecně text na řádku za % je ignorován a slouží jako komentář pro lepší čitelnost kódu. Blok komentářových řádků na začátku funkce se dále využije u příkazů help, doc nebo lookfor, které umožní zobrazit návod ke každé funkci.
• U vlastního kódu si všimněte, že jednotlivé příkazy jsou zakončeny středníkem ;, aby funkce nevypisovala mezivýsledky na obrazovku.
• Při výpočtu Antoineho rovnice použijeme teplotu převedenou na Kelviny a konstanty, zapsané přímo v programu. Výsledek uložíme do výstupní proměnné p0, jejíž jméno jsme určili v hlaviččce funkce.

Funkci teď můžeme vyzkoušet na příkazovém řádku.

>> help phexan
>> lookfor 'pressure'
>> phexan(60)

Vidíme, že naše první funkce funguje a dokonce si k ní můžeme zobrazit dokumentaci.

Funkce s maticemi jako proměnnými

Matlab je navržen tak, aby pracoval s vektory a maticemi a proto i většina vestavěných funkcí je uzpůsobena k tomu, že její proměnné nemusí být pouze skaláry, ale také vektory nebo matice. Hodí se to například, když chceme nakreslit graf funkce sinus

>> x = linspace(0,2*pi,100);
>> y = sin(x);
>> plot(x,y,'r-o')

Prvním příkazem pomocí funkce linspace vytvoříme řádkový vektor se 100 čísly, ekvidistantně rozloženými mezi 0 a 2π. Vidíme, že tento vektor x může být vstupní proměnná funkce sin, která vrátí stejně veliký vektor y. Každému prvku z vektoru x odpovídá prvek z vektoru y, který má hodnotu sin(x). Třetí příkaz nakreslí graf. Řetězec 'r-o' znamená červené (‘r’) kroužky (‘o’) spojené plnou čárou (‘-‘).

Podobně bychom mohli chtít do grafu zakreslit závislost tenze par hexanu na teplotě v rozmezí 20 až 105 °C

>> t = linspace(20,105,100);
>> p0 = phexan(t);
>> plot(t,p0,'r-o')

Vidíme, že naše funkce, která fungovala pro jednu teplotu, teď pro vektor více teplot v proměnné t vrací chybu. Pokud chceme funkci opravit, aby fungovala i pro proměnné, které jsou vektory nebo matice, musíme si znovu projít kód funkce a uvažovat o proměnné t jako o vektoru / matici.

• Nejdříve k matici t přičteme skalár.

  tK = t+273.15;

  tK = t+273.15;

  I když jedné straně operátoru + je matice a na druhé skalár, Matlab tento výraz interpretuje tak, že ke každému prvku matice přičte tento skalár. Výsledkem je matice tK o stejném tvaru jako má proměnná t.

• Na dalším řádku k tK přičteme skalár C. Potom se ale pokoušíme dělit skalár B maticí (tK+C)

  p0 = exp(A-B/(tK+C));

  p0 = exp(A-B/(tK+C));

  Operátor /, který dělení maticí interpretuje jako násobení inverzní maticí, ale v našem případě vrátí chybu, že operaci nelze provést.

• Oprava kódu je jednoduchá. Námi požadovaného chování dosáhneme použitím “tečkového” operátoru ./, takže celý řádek bude vypadat

  p0 = exp(A-B./(tK+C));

  p0 = exp(A-B./(tK+C));

• Zbytek kódu pak už bude fungovat i pro matice: zbývá odečíst matici od skaláru A, čímž vznikne opět matice. Nakonec funkce exp je uzpůsobena na maticové proměnné podobně jako funkce sin v předcházejícím příkladu.

Přidáním jedné tečky do kódu naší funkce phexan jsme ji upravili tak, že bude fungovat i pro maticové proměnné a jedním voláním dostaneme tenzi par hexanu při různých teplotách. Poznamenejme ještě, že podobně jako na operátor dělení je potřeba dávat pozor i na operátor umocňování, kdy ve většině případů budete chtít použít tečkový .^ místo maticového ^.

Univerzální funkce pro výpočet tenze par

Námi vytovřená funkce phexan funguje skvěle, ale není příliš praktické mít pro každou látku zvláštní funkci. Proto si udělejme univerzální funkci ptension, která počítá tenzi par pomocí Antoineho rovnice, která bude vypadat takto

function [p0] = ptension(A,B,C,t)
%PTENSION Calculate vapour pressure
%   Vapour pressure calculation using Antoine equation.
%
%   P0 = PTENSION(A,B,C,T)
%
%   A,B,C are Antoine equation constants,
%   T is temperature in Celsius and
%   P0 is vapour pressure in kPa.
%
tK = t+273.15;          % convert temperature to Kelvins
p0 = exp(A-B./(tK+C));  % vapour pressure in kPa
end

a kterou musíme uložit do souboru ptension.m do pracovní složky. Vidíme, že kromě teploty jsme jako vstupní proměnné museli přidat koeficienty Antoineho rovnice. Protože hledání a zadávání těchto konstant před každým použitím funkce je pracné, vytvoříme si funkci ptension_coef, ve které budou tyto koeficienty uloženy

function [A, B, C] = ptension_coef()
%PTENSION_COEF Antoine equation coefficient database
%   Supply Antoine equation coefficients for PTENSION function
%
%   [A,B,C] = ptension_coef()
%
%   In upstream routine call:
%   >>[A, B, C] = ptension_coef();
%
%   source data: CHI-tabulky (Holecek), str. 53-55
A(1)=13.8216; % n-hexane (valid between 245-370 K)
B(1)=2697.55;
C(1)=-48.78;
 
A(2)=13.8744; % n-heptane (valid between 265-400 K)
B(2)=2895.51;
C(2)=-53.97;
 
A(3)=13.9276; % n-octane (valid between 292-425 K)
B(3)=3120.29;
C(3)=-63.63;
 
end

Tyto řádky kódu musí být uloženy do souboru ptension_coef.m. Všimněte si, že tato funkce nemá žádné vstupní proměnné, tj. seznam vstupních proměnných je prázdný “()”. Vidíme, že nám vrátí tři výstupní proměnné: A, B a C což jsou vektory, které obsahují příslušné konstanty pro tři složky, které máme zatím v “databázi”: hexan, heptan a oktan. Budeme-li chtít například tenzi par heptanu při 60 °C zadáme

>> [A, B, C] = ptension_coef()
>> ptension(A(2),B(2),C(2),60)


Připomeňme si, že výraz A(2) znamená druhý prvek vektoru A, ve kterém máme uložen koeficient pro heptan.

Když si prostudujete kód, zjistíte, že funkce ptension funguje i pro případ, kdy proměnné A, B a C jsou vektory a proměnná t skalár. Volání

>> ptension(A,B,C,60)

nám tedy vráti tenzi par všech tří složek v databázi při teplotě 60 °C.

[Zatím nemúžeme mít vektory jako vstupní proměnné pro konstanty i pro teplotu. V tuto chvíli nám to nevadí, příslušná úprava kódu pro ptension by vyžadovala použit například cykly, což už je mimo rámec tohoto článku.]

Dvojici funkcí ptension a ptension_coef budeme používat i v dalších příkladech. Všimněte si, že funkce ptension_coef slouží jako základ “databáze”, kterou můžeme v případě potřeby doplňovat a vylepšovat. Například bychom mohli přidat vstupní parametr, kterým bychom zadali jméno látky a dostali jí odpovídající koeficienty. Fungovalo by to takto

>> [A, B, C] = ptension_coef('hexan')

Taková úprava už ale jde mimo rozsah předmětu “Modelování”, v našich příkladech a úlohách si vystačíme s jednoduchou databází, kdy si pořadí jednotlivých složek budeme pamatovat.

A teď úplně něco jiného: skripty

Stejně jako funkce, tak i skripty se skládají z řádků kódu uložených v souboru s koncovkou “.m”. Jestliže funkce má vlastní jmenný prostor pro své lokální proměnné, tak skript používá stejný jmenný prostor jako prostředí ze kterého je spuštěn. Skripty, které vytvoříme, budeme vždy spouštět z příkazového řádku. Pokud ve skriptu definujeme nebo změníme nějakou proměnnou, tak hodnotu této proměnné “uvidíme” i po ukončení skriptu. Zároveň skript “uvidí” a může pracovat se všemi proměnnými, které jsme předtím vytvořili na příkazovém řádku nebo příkazem v předtím spuštěném skriptu. To znamená, například, že skript spuštěný podruhé “vidí” proměnné, které si sám při prvním spuštění vytvořil.

Pokud postupně zadáváme příkazy a spouštíme vestavěné i námi vytvořené funkce z příkazového řádku, tak můžeme vyřešit všechny úlohy i bez skriptů. Se skriptem si potřebné příkazy můžeme předem zapsat do souboru (tj. vytvořit skript) a tento skript potom spustit. Výhoda skriptu je v tom, že když uděláme chybu nebo chceme něco pozměnit, nemusíme znovu pracně zadávat příkazy do příkazového řádku, ale stačí pozměnit skript a znovu jej spustit. Další výhoda je, že si nemusíme přesně pamatovat příkazy potřebné k vyřešení konkrétní úlohy, protože je máme zapsány ve skriptu a můžeme se k ni vrátit i po delší době.

Jestli je kód v souboru s koncovkou “.m” skript nebo funkce poznáme z prvního řádku souboru. Funkce vždy na prvním řádku obsahuje klíčové slovo function, které definuje její rozhranní. Jinak je to skript. V předmětu “Modelování” bude každý problém tvořen jedním skriptem, kterému budeme říkat hlavní program a jednou nebo více námi vytvořenými funkcemi. Při zkoušení a ladění funkcí budeme zadávat příkazy na příkazovém řádku, ale jakmile vše bude fungovat, sepíšeme tyto příkazy do hlavního programu / skriptu. Spuštění hlavního programu / skriptu povede k vyřešení příslušné úlohy.