Úvod do Programování Lukáš Stříteský, Marek Pechal 20. prosince 2007 Abstrakt V článku se dozvíte, jak začít programovat v jazyce Pascal. Jediné předpokládané znalosti jsou, že dokážete z internetu stáhnout program a nainstalovat si jej na svůj počítač. Vše ostatní se zde dočtete. K úspěchu stačí chuť dozvědět se něco nového, kterou jistě máte. Naučíme se také používat „malováníÿ zvané Gnuplot. Závěrem si objasníme i některé základní numerické metody. A také jak naše výsledky elegantně publikovat.
1
Hned na to
1.1
První program
Určitě jste netrpěliví a rádi byste už začali psát svůj první program. Nic nám v tom nebrání. Jediné, co je k tomu potřeba stáhnout, například z FYKOSího webu1 je program FreePascal, takzvaný kompilátor (překladač) Pascalu. Jeho domovská stránka je http://www.freepascal.org/, kde se o něm případně dozvíte spoustu dalšího.2 Máte-li jej stažený, stačí si ho jednoduše nainstalovat (popřípadě předem ještě rozbalit) a spustit. Klepněte na první záložku nahoře „Fileÿ a potom klepněte na „Newÿ, jak je znázorněno na obrázku. Otevřelo se nám nové okno, napíšeme zde náš první program, opište do něj prosím tento text. Program Ahoj; Begin WriteLn(’Ahoj!’); ReadLn; End. Potom stiskněte Ctrl-F9 (popřípadě „Runÿ a opět „Runÿ). FreePascal vás vyzve k uložení Vašeho programu. Stačí napsat „ahojÿ a stisknout ENTER, popřípadě kliknout na „Okÿ. Náš první program je na světě, navíc se okamžitě spustil. Otevřela se konzole, ve které je napsáno něco jako Running c:\fcp\2.1.4\bin\i386-win32\ahoj.exe Ahoj! 1
http://fykos.mff.cuni.cz nebo odkudkoli jinde Stáhněte si prosím verzi FreePascalu 2.0.0 (ta je na webu FYKOSu). Současná verze tohoto programu 2.0.4 (červen 2007) v sobě bohužel obsahuje nepříjemnou chybu způsobující nefunkčnost unity Graph, která by vás třeba mohla lákat. Nebo by se dala použít ještě i verze 2.2.0-beta (tzv. 2.1.4 aka), tam už je opět vše v pořádku a stabilní verze by měla být hotova v srpnu 2007. 2
1
Obrázek 1: Nový soubor ve FreePascalu
2
Gratuluji! Naprogramovali jste svůj první program v Pascalu! Jednoduché, že? Teď si vysvětlíme, co se to přesně stalo, co náš program zatím umí a naučíme ho toho ještě trochu více. Pascalovské programy začínají klíčovým slovem Program, za kterým je uveden jeho název (hlavička programu)3 . Potom následuje deklarační část, kterou jsme prozatím vypustili, neboť ji zatím nepotřebujeme. Poslední část je příkazová označená slovy Begin a End. Důležitá vlastnost Pascalu je, že není takzvaně case-sensitive, čili nerozlišuje velká a malá písmena. Napsali jsme dva příkazy (každý příkaz v Pascalu je ukončen středníkem „;ÿ). První příkaz je WriteLn(’Ahoj!’); Tento příkaz vypíše na obrazovku Ahoj!. Druhý příkaz ReadLn; slouží k tomu, aby se do programu načetl vstup z klávesnice. Napište v příkazové řádce našeho programu cokoliv a stiskněte ENTER. Do programu se načetl váš vstup, čímž se příkaz ReadLn; vykonal. A protože je to příkaz pro tuto chvíli poslední, tak se program ukončil. Pojďme jej tedy naučit mnohem více, třeba počítat. Bude to opět velmi jednoduché.
1.2
Ciferný součet
Nyní si napíšeme program, který nám spočítá ciferný součet libovolného celého čísla. Například ciferný součet čísla 1987 je 1 + 9 + 8 + 7 = 25. Opět si prosím přepište4 zdrojový kód a zkuste si jej spustit. Vzápětí si jej vysvětlíme. Program CifernySoucet; Var X, Y: Integer; Begin WriteLn(’CIFERNY SOUCET’); WriteLn(’Napiste cislo a~stisknete Enter’); ReadLn(X); Y := 0; While X <> 0 do Begin Y := Y + X mod 10; X := X div 10 End; WriteLn(Y); ReadLn; end. Co je zde všechno nové? Řádek Var X, Y: Integer; nám takzvaně deklaruje dvě proměnné (angl. variables) X a Y. Co to znamená? Řekneme tím našemu programu, ať si zamluví určitou část paměti počítače pro nějaké dvě hodnoty, jež jsme si označili X a Y. Budou to dvě celá čísla, sdělili jsme to klíčovým slovem Integer. Příkaz ReadLn(X); načte vstup z klávesnice do naší již deklarované proměnné X. Čili vezme námi naťukané číslo a uloží ho do paměti počítače. Příkaz Y := 0; přiřadí proměnné Y hodnotu 0. Příkaz While X <> 0 do poručí našemu počítači: Prováděj níže uvedené výpočty (příkazy), dokud se X nebude rovnat nule! Je-li příkazů více, v našem případě dva, uzavřeme je opět mezi Begin a End;. 3
Je to nepovinný parametr. Asi jste si již všimli, že ve FreePascalu nám nefunguje oblíbená sekvence CTRL-C, CTRL-V, neboť se jedná o DOSovský program, nikoliv „Wokenníÿ. Namísto těchto klávesových zkratek zkuste CTRL+INSERT, resp. SHIFT+INSERT. Druhou možností je si náš zdrojový soubor, který máme uložený, třeba ahoj.pas, otevřít pomocí Poznámkového bloku. Provedeme v něm požadované změny a uložíme. Novou upravenou verzi načteme do FreePascalu stiskem CTRL-ALT-F a R nebo „Fileÿ a „Reloadÿ. 4
3
Tyto složené přiřazovací příkazy obsahují pro nás dva nové výrazy s operátory mod a div. Ty umí vydělit dvě čísla, avšak celočíselně se zbytkem, jak děti v první třídě5 . Prvním výrazem je X mod 10. Tento výraz nabude hodnoty (říká se spíše, že vrátí hodnotu) zbytku po celočíselném dělení čísla X desíti. Například pro X = 1987 je X mod 10 = 7 nebo 11 mod 3 = 2 apod. Naopak výraz X div 10 vrátí samotnou hodnotu celočíselného dělení desíti. Pro X = 1987 je X div 10 = 198 nebo 11 div 3 = 3 a pod. Ukažme si nyní na konkrétním příkladu, co náš program dělá. Dejme tomu, že uživatel naťuká 1987 a stiskne ENTER. Toto číslo se uloží do proměnné X. Proměnná Y má hodnotu 0. Poté vstoupíme do takzvaného podmíněného (While) cyklu. Cyklicky, tedy pořád dokola, budeme vykonávat příkazy Y := Y + X mod 10; a X := X div 10; dokud nenabude proměnná X hodnoty 0. A co že to ty příkazy provádí? Nejprve přiřadí do proměnné Y samotnou hodnotu Y (nyní nula) plus 1987 mod 10, čili sedm. Proměnná Y tedy nabyla hodnoty sedm. Do proměnné X se na dalším řádku přiřadí hodnota 1987 div 10, čili 198. Protože se X stále ještě nerovná 0, provedeme tyto příkazy znovu. Řádek Y := Y + X mod 10; nyní znamená Y =7 + 198 mod 10 = 15, čili proměnná Y má už hodnotu patnáct. Proměnná X vzápětí nabude hodnoty 198 div 10 = 19. Protože to však stále ještě není 0, tak musíme počítat zase znova. Y = 15 + 19 mod 10 = 24 a X = 19 div 10 = 1. Jednička už je hodně málo, avšak stále ještě ne 0, proto proběhne cyklus ještě naposledy a jeho výsledkem je, že Y = 24 + 1 mod 10 = 25 a X = 1 mod 10 = 0. Náš vydřený výsledek se vypíše příkazem WriteLn(Y); Příkaz ReadLn; nám nyní slouží vlastně pouze k tomu, aby se program neukončil okamžitě po vypsání výsledku, ale nechal nás si jej alespoň přečíst nebo vyfotografovat.
1.3
Rozklad na prvočísla
Nyní si napíšeme program, který nám rozloží číslo na jeho prvočísla. Například 180 = 2 · 2 · 3 · 3 · 3 · 5. Tento program vypadá následovně. Program Prvociselny_Rozklad; Var N: Integer; D: Integer;
// zkoumane cislo // delitele
Begin Write(’Zkoumane kladne cele cislo vetsi nez 1: ’); Readln(N); Write(’Prvociselny rozklad: ’); While N mod 2 = 0 do Begin N:=N div 2; Write(2, ’ ’) End; D:=3; // prvni lichy delitel While N > 1 do Begin While N mod D = 0 do Begin N:=N div D; 5
Přece jenom je náš program stále ještě takové dítě.
4
Write(D, ’ ’) End; D:=D+2 End; ReadLn; End.
1.4
// dalsi mozny delitel
Hřeb – naše nová kalkulačka
Následujícím programem je kalkulačka. Program Kalkulacka; Var vyraz: String; //najde polohu operatoru s~nejvetsi prioritou //a odstrani prebytecne vnejsi zavorky Function Poloperatoru(var vyraz: String): Integer; Var plus, krat: Integer; zavorky: Integer; i: Integer;
//poloha posledniho operatoru //bilance zavorek
begin plus:=0; krat:=0; zavorky:=0; for i:=1 to length(vyraz) do case vyraz[i] of ’(’ : zavorky := zavorky + 1; ’)’ : zavorky := zavorky - 1; ’+’, ’-’ : if zavorky = 0 then plus := i; ’*’, ’/’ : if zavorky = 0 then krat := i; end; if plus > 0 then Operator := plus else if krat > 0 then Operator := krat else if vyraz [1] <> ’(’ then Operator := 0 else //odstranime vnejsi zavorky begin vyraz := Copy(vyraz, 2, length(vyraz) - 2); Operator := Operator(vyraz) end end; //function Operator
5
//spocita rekurzivne hodnotu zadaneho vyrazu Function Hodnota(var vyraz: string): integer; Var poloha: vyraz1, vyraz2: vysledek:
Integer; String; Integer;
//poloha operatoru //levy a~pravy usek od operatoru //vysledna hodnota
begin poloha := Poloperatoru(vyraz); if poloha = 0 then Val(vyraz, vysledek); else begin vyraz1 := Copy(vyraz, 1, poloha-1); vyraz2 := Copy(vyraz, ploha + 1, Length(vyraz)-poloha); case vyraz[poloha] of ’+’ : vysledek := Hodnota(vyraz1) + Hodnota(vyraz2); ’-’ : vysledek := Hodnota(vyraz1) - Hodnota(vyraz2); ’*’ : vysledek := Hodnota(vyraz1) * Hodnota(vyraz2); ’/’ : vysledek := Hodnota(vyraz1) div Hodnota(vyraz2); end end; Hodnota := vysledek; end; Begin WriteLn(’Vitejte v~Kalkulatoru!’); WriteLn(’Prikazem {,,}+{‘‘} muzete scitat’); WriteLn(’Prikazem {,,}*{‘‘} muzete nasobit’); WriteLn(’Prikazem {,,}-{‘‘} muzete odcitat’); WriteLn(’Prikazem {,,}/{‘‘} muzete delit’); WriteLn(’Uzivejte zavorek {,,}({‘‘}, {,,}){‘‘} dle libosti’); WriteLn(’Jen at Vas vyraz neprekroci delku 256 znaku.’); ReadLn(vyraz); WriteLn(’Hodnota vyrazu je: ’); WriteLn(Hodnota(vyraz)); ReadLn; end.
2
Teorie
Nyní, poté co jsme si programování osahali na praktických příkladech, by vám mohla přijít vhod také nějaká ta teorie. Studujme na chvíli tedy jazyk Pascal.
6
2.1
Programovací jazyky
Nejstaršími programovacími jazyky jsou jazyky strojové a jejich symbolické verze (asemblery). Této skupině jazyků se říká souhrnně strojově orientované jazyky. Tyto programy se skládají z posloupnosti elementárních příkazů, které nazýváme instrukce a které může přímo provádět procesor. I dlouhé programy napsané takovýmto jazykem jsou tedy dobře zvládnutelné počítačem, ovšem o to obtížnější je práce programátora. Hlavním důvodem, proč se v nich tak špatně vytváří programy je, že nejsou strukturované. Tento nedostatek odstranily vyšší programovací jazyky, mezi které patří i Pascal. Povíme si nyní něco o jejich základních rysech. 2.1.1
Typy dat
Ve vyšších programovacích jazycích jsou vymezené typy dat, které určují jejich obor hodnot a možné operace nad nimi. Přirozeně to bude například typ celých čísel (v Pascalu Integer), reálných čísel (v Pascalu Real) nebo posloupnost alfanumerických znaků (v Pascalu String). 2.1.2
Konstanty
Konstanta je datový objekt, jehož hodnota se během celého výpočtu nemění. Podle toho, jak konstantu v programu identifikujeme, rozlišujeme dva druhy konstant. První jsou přímé konstanty neboli literály. Je jí přímý zápis hodnoty do programu, například 2,718282 nebo 3,14593 nebo ’posloupnost neboli řetězec alfanumerických znaků’. Pojmenované konstanty se identifikují svými svými jmény, takzvanými identifikátory a v programu se zavádějí deklarací. 2.1.3
Proměnné
Proměnná je datový objekt, jehož hodnota se může v průběhu výpočtu měnit. Proměnné se nejčastěji identifikují pomocí jmen, takzvaných identifikátorů proměnných. Ty se opět zavádí deklaracemi, čímž se určí též jejich datový typ. 2.1.4
Deklarace
Deklarace zavádí a pojmenovává určité součásti programu jako proměnné, konstanty, funkce, procedury a podobně. Deklarace se můžou lišit podle rozsahu jejich platnosti. Takto rozlišujeme deklarace globální, jenž platí v celém programu a deklarace lokální, které se uplatní pouze v jedné jeho části například pouze v jedné funkci. Deklarace lze rovněž dělit na statické a dynamické. Rozdíl je v tom, že dynamická deklarace může být ovlivněna za běhu programu například tím, jaký programu pošleme vstup. Naopak statická deklarace nezávisí na hodnotách proměnných a uskuteční se za všech okolností stejně. Spousta vyšší programovacích jazyků nepřipouští dynamické deklarace a Pascal mezi ně patří. 2.1.5
Výrazy
Výraz je operační struktura, pomocí které se popisuje výpočet hodnoty. V programovacích jazycích mají většinou podobný tvar jako matematické výrazy například X := ( 4 + 5 ) * 6. Skládají se z operátorů, operandů i závorek. Operandem může být konstanta, proměnná nebo hodnota funkce. Operátory se obvykle dělí na unární a binární podle počtu operandů a mohou působit pouze na přípustné typy operandů. Dále jsou pro ně určeny typy výsledných hodnot a také priority, podle kterých se určuje pořadí při provádění operací ve výrazu, v němž se vyskytuje více než jeden operátor. 7
Mezi operátory patří i relační operátory, kterými se tvoří relace. Relace je výraz, který nabývá jedné z logických hodnot6 podle toho, zda hodnoty operandů v něm uvedených jsou či nejsou ve vztahu, který udává relační operátor. Ve většině vyšších programovacích jazyků musí být výraz vždy součástí nějaké vyšší jazykové konstrukce, zpravidla příkazu, z níž je zřejmé, co má být s hodnotou výrazu provedeno, například přiřazena proměnné. 2.1.6
Příkazy
Pomocí příkazů se popisují jednotlivé výpočetní akce. Existují příkazy například přiřazovací (které přiřazují hodnotu proměnné), přikazy skoku, volání funkcí a procedur. Řídícími strukturami, kterými se vyjadřuje postupné, podmíněné či opakované provádění dílčích příkazů, jsou strukturované příkazy. V jazycích, které navazují na koncepci příkazů jazyka Algol 60, to je složený příkaz, podmíněné příkazy a příkazy cyklu. Většinou se skládají z dílčích příkazů a obsahují také výrazy udávající podmínky pro větvení či opakování. Ilustrativním příkladem je přiřazení proměnné Y absolutní hodnotu proměnné X. if X >= 0 2.1.7
then Y := X
else X := - Y
Syntaktický diagram
Syntaxe je souhrn pravidel udávajících přípustné tvary dílčích konstrukcí a celého programu. Každý syntaktický diagram definuje všechny přípustné tvary konstrukce, jejíž název je uveden v nadpisu diagramu. Každý přípustný tvar této konstrukce získáme tak, že projdeme nějakou cestou vedoucí od začátku syntaktického diagramu do jeho konce. identifikator
pismeno
pismeno
cislice
2.1.8
Podprogramy
Při hledání algoritmického řešení složitějšího problému si jej rozdělíme na několik podproblémů, kterým se pak věnujeme každému zvlášť a nakonec je všechny propojíme. Vyšší programovací jazyky umožňují aplikovat takovýto přístup při psaní programu jeho rozkladem na podprogramy. Podprogram je tedy popisem dílčího algoritmu, v němž konkrétní data mohou být zastoupena formálními parametry. Podprogramy se dělí na procedury a funkce. Procedura je popisem algoritmu, který nějakým způsobem transformuje data, funkce popisuje algoritmus výpočtu jedné hodnoty, která může být též použita jako operand ve výrazu. 6
pravda či nepravda
8
2.2
Jazyk Pascal
V této kapitolce si povíme něco konkrétně k jazyku Pascal. 2.2.1
Jazyk Pascal
2.2.2
Co s výstupem?
Jakožto fyzikům nám občas nebude stačit výstup z našeho programu na obrazovku, protože budeme potřebovat získaná data dále zpracovávat (např. z nich vytvářet grafy – v další kapitolce si ukážeme jak). Naštěstí nám operační systém (budeme teď chvíli mluvit konkrétně o Windows) umožňuje to, co program vypisuje, velmi jednoduše zapsat do textového souboru. Jak?7 Spustíme si příkazový řádek – v nabídce Start systému Windows zvolíme spustit, do textového pole zadáme cmd a stiskneme ENTER, případně klikneme na OK. Otevře se nám okno s černým pozadím – windowsovský příkazový řádek. Na posledním řádku je vypsán název aktuálního adresáře. Ten musíme nejdříve změnit na adresář s naším programem (předpokládejme, že je to třeba C:\mujadresar\programky). V adresářové struktuře se můžeme pohybovat pomocí příkazu cd. Zápisem cd.. se přesuneme o adresář výš (je-li to možné), příkazem cd adresar vstoupíme z aktuálního adresáře do jeho podadresáře s názvem adresar. Je-li počáteční aktuální adresář C:\Documents and Settings\Moje, budeme postupovat takto: C:\Documents and Settings\Moje>cd.. C:\Documents and Settings>cd.. C:\>cd mujadresar C:\mujadresar>cd programky C:\mujadresar\programky> Teď bychom chtěli přesměrovat výstup programu fykos.exe do textového souboru vystup.txt. Do příkazového řádku jednoduše zadáme fykos > vystup.txt. To, co by program normálně vypisoval na obrazovku, se nyní tiše zapíše do textového souboru. Tohle celé lze samozřejmě zařídit i programově – tedy přímo v pascalovském programu automaticky zapisovat do souboru (nemusíme se pak obtěžovat s příkazovým řádkem) – a není to až tak složité. Posuďte sami. Než začneme zapisovat, musíme si deklarovat proměnnou, která bude náš výstupní soubor reprezentovat. Tato proměnná bývá typu text či textfile (podle toho, jakou verzi Pascalu používáte). Této proměnné pak musíme přiřadit název souboru a ten pak otevřít. ... var f:text; begin assign(f,’nazevmehosouboru.txt’); rewrite(f); ... Příkaz rewrite přepíše veškerý dosavadní obsah našeho souboru (vyprázdní jej), a pokud soubor neexistuje, tak jej vytvoří. Pokud byste chtěli starý obsah souboru zachovat a pouze k němu připisovat další data, museli byste pro otevření použít příkaz append. V tomto okamžiku můžete do souboru zapisovat stejně jako na obrazovku – tedy příkazy write a writeln, pouze musíte jako jejich první parametr udat identifikátor souboru. Například: 7
Pokročilejším čtenářům se omlouváme, za trochu velmi základního opakování.
9
Obrázek 2: Okno právě spuštěného Gnuplotu verze 4.0 writeln(f,’Vypisu tam text a~cislo ’,x,’ a~zase text.’); Když už se souborem nehodláte dále pracovat, je dobré jej za sebou zavřít příkazem close(f).
3
Gnuplot
Program Gnuplot je velmi užitečný nástroj, určený – jak už jeho název napovídá – ke kreslení všemožných typů grafů. Předpona GNU nám také říká, že program je vyvíjen pod stejnojmennou opensourceovou licencí a je tedy zdarma ke stažení (i s obsáhlou dokumentací) na internetových stránkách http://www.gnuplot.info. Po jednoduché instalaci a spuštění zjistíte, že se Gnuplot chová podobně jako například DOSovský příkazový řádek – bliká na vás kurzor a čeká na váš příkaz, který je po stisknutí klávesy ENTER ihned vykonán. Nepíšete tedy celý program, který by se pak najednou provedl, jak je tomu v Pascalu (i když i to se dá v Gnuplotu dělat, jak uvidíme později). Dalšími důležitými rozdíly Gnuplotu oproti Pascalu, jejichž opomíjení se může škaredě vymstít8 , jsou jeho „case sensitivnostÿ (tzn. je důležité, jestli napíšete plot, nebo Plot, zatímco v Pascalu je to prašť jako uhoď) a jeho chování při dělení, které je pro uživatele neobeznámeného s jazykem C téměř zákeřné. Napíšete-li totiž příkaz plot [0:1] 1/3+x a budete očekávat, že vám Gnuplot poslušně naservíruje graf funkce 1/3 + x v intervalu [0, 1], překvapí vás vykreslením grafu funkce x a vy budete jistě dumat, jaký čert v tom vězí. Svízel je v tom, že Gnuplot používá znak / pro normální i pro celočíselné dělení a to, o které z nich se v tom kterém případě jedná, pozná z charakteru operandů. Jsou-li oba celočíselné, myslí si Gnuplot, že jde o dělení celočíselné, a tak konstantu 1/3 nevyhodnotí jako 0, 333 . . ., jak bychom zajisté chtěli, nýbrž jako prachsprostou nulu. To, že po něm žádáme jednu třetinu, oznámíme například zápisem9 8
Vzpomeňte si, až budete zoufale lovit chybu, která – jak se ukáže – spočívala v jediném velkém písmenu na místě malého. 9 Stejně dobře by fungovalo i 1/3.0 nebo dokonce 1./3 či 1/3..
10
Obrázek 3: Gnuplotí okna s výslednými grafy 1.0/3. Nyní už nejsou oba operandy celočíselné, a tak Gnuplot interpretuje zápis správně.
3.1 3.1.1
Základní příkazy Příkaz plot
Jistě nejdůležitějším příkazem, který budete v Gnuplotu používat, je příkaz plot, který vykreslí graf zadané funkce nebo data z externího souboru. Jeho povinným parametrem je interval nezávislé proměnné (označené implicitně jako x), v němž se graf vykreslí, a samozřejmě samotná funkce. Interval se udává ve tvaru [a:b], kde a a b jsou jeho krajní body, a funkce ve formě číselného výrazu, podobného tomu z Pascalu (používané operátory a funkce blíže popíšeme později). Použití funkce plot může vypadat například takto: plot plot plot plot
[0:1] x [-pi:pi] sin(x) [-1:1] sqrt(x**2+1) [2:3] sqrt(exp(x)/x**2+(1-x)*log(x))
Po zadání některého z těchto řádků do Gnuplotu a potvrzení ENTERem se otevře nové okno s výsledným grafem (viz obrázek 3) Konvence pro zápis výrazů v Gnuplotu je téměř stejná jako v Pascalu, pro elementární aritmetické operace se užívá operátorů +, -, *, /, jejich priorita je zcela obvyklá (tj. násobení a dělení má vyšší prioritu než sčítání a odčítání), stejně tak závorky fungují tak, jak jsme zvyklí. Oproti Pascalu existuje navíc binární operátor ** pro umocňování. Symboly matematických funkcí začínají malým písmenem, jejich argumenty se uzavírají do kulatých závorek. Několik nejzákladnějších uvádíme v tabulce 1. Pokud chceme, můžeme do jednoho grafu vykreslit i více funkcí. Stačí je od sebe oddělit čárkami, např. plot [0:1] x,x**2,x**3. Příkaz plot také umí vykreslovat do grafu data ze souboru, což jistě velmi oceníte, až se (možná) budete na matfyzu ve druhém semestru protloukat fyzikálním praktikem. Ukážeme si to na příkladu. Mějme textový soubor teplota.txt, obsahující ve dvou sloupcích10 data z měření průběhu venkovní 10
oddělených mezerou či tabulátorem
11
Tabulka 1: Základní matematické funkce v Gnuplotu druhá odmocnina exponenciála přirozený logaritmus sinus (resp. hyperbolický) kosinus (resp. hyperbolický) tangens (resp. hyperbolický) arkussinus (resp. hyperbolický) arkuskosinus (resp. hyperbolický) arkustangens (resp. hyperbolický)
sqrt exp log sin (resp. sinh) cos (resp. cosh) tan (resp. tanh) asin (resp. asinh) acos (resp. acosh) atan (resp. atanh)
Obrázek 4: Plot dat z externího souboru teploty během dne (v prvním sloupci čas uplynulý od půlnoci v hodinách, ve druhém teplota ve stupních Celsia). Chceme-li z našeho měření vykreslit graf, jednoduše Gnuplotu zadáme příkaz plot ’teplota.txt’. Než to ale uděláme, musíme správně nastavit pracovní adresář Gnuplotu na ten, ve kterém se soubor teplota.txt nachází. To uděláme pomocí příkazu cd (napíšeme např cd ’c:\mujadresar’, je-li náš soubor umístěn na disku c: v adresáři mujadresar). Pak už nám nic nebrání přistoupit k samotnému vykreslování. Výsledek ukazuje obrázek 4. Dejme tomu, že v souboru teplota.txt máme kromě sloupce s časem a teplotou ještě třetí sloupec s naměřenými hodnotami atmosférického tlaku, jehož závislost na denní době bychom také chtěli vykreslit. Musíme tedy Gnuplotu nějak sdělit, že nechceme do grafu vykreslovat data z prvního a druhého sloupce, nýbrž z prvního a třetího. To uděláme připojením textu using 1:3. Náš příkaz tedy bude plot ’teplota.txt’ using 1:3. Klíčové slovo using můžeme také využít v situaci, kdy nechceme na svislou osu grafu vynášet přímo data ze souboru, ale nějakou hodnotu na nich závislou, například jejich druhou mocninu. Pak napíšeme using 1:($2**2). Symbol $2 zde zastupuje data z druhého sloupce (a pozor – závorky zde nelze vynechat!). Stejný trik samozřejmě funguje i pro vodorovnou osu. Po zadání plot ’data.txt’ using (log($1)) :($1*$4) bychom získali graf, v němž na vodorovné ose je vynášen přirozený logaritmus čísel z prvního sloupce souboru data.txt, na svislé pak součin čísel 12
Obrázek 5: Plot dat z externího souboru s chybovými úsečkami z prvního a čtvrtého sloupce. Vraťme se nyní zpět k příkladu s výsledky měření teploty a tlaku. Jestliže bychom chtěli v grafu závislosti tlaku na čase vykreslené body pospojovat čárou (o to se ale v praktiku nikdy nepokoušejte, jinak skončíte ve fyzikálním pekle), jednoduše Gnuplotu řekneme plot ’teplota.txt’ using 1:3 with lines. Jsme-li jako experimentátoři alespoň trochu co k čemu, víme, že spolu s naměřenými hodnotami bychom měli uvést i jejich chyby. Mějme tedy ve čtvrtém sloupci našeho souboru zapsány chyby jednotlivých měření teploty. Ty dostaneme do grafu ve formě hezkých chybových úseček (errorbarů) příkazem plot ’data.txt’ using 1:2:4 with errorbars, kde tentokrát musíme specifikovat tři sloupce dat – první s nezávisle proměnnou, druhý se závisle proměnnou a třetí s její chybou. Výsledek bude vypadat jako na obrázku 5. Nakonec je ještě třeba zmínit se o souvisejícím příkazu replot, který – zadán bez parametrů – zopakuje poslední provedený plot. Dokáže však do něj také přidat další graf, zadáme-li mu jako parametr požadovanou funkci či název souboru s daty (případně doplněné o další specifikace jako např. using 1:3 či with lines). Všimněme si však, že tentokrát již nezadáváme interval nezávislé proměnné, protože ten je dán z posledního provedení příkazu plot. 3.1.2
„Pomocnéÿ příkazy
Samozřejmě se nemusíme spokojit se vzhledem grafu, který nám Gnuplot poskytne na první pokus. Nabízí se nám totiž řada možností, jak jej dále upravovat – můžeme přidat popisky os, legendu jednotlivých vykreslených závislostí, atd. To lze provádět nastavením nejrůznějších parametrů před samotným vykreslením grafu. Potřebný příkaz má většinou tvar set nazevparametru jehohodnota, chceme-li nastavit parametr nazevparametru na hodnotu jehohodnota, nebo unset nazevparametru, přejeme-li si jej vynulovat. Je-li tedy například naším cílem přidat popisek vodorovné osy t [ms], napíšeme příkaz set xlabel ’t [ms]’. Obdobně pro svislou osu užijeme parametr ylabel. Pokud se nám popisek nelíbí, můžeme jej opět odstranit příkazem unset xlabel. Gnuplot nám také umožňuje upravovat legendu vykreslovaných závislostí, je-li jich v jediném grafu více. Její text nastavíme při volání příkazu plot přidáním klíčového slova title následova13
Obrázek 6: Dvě závislosti v jediném grafu opatřené legendou ného požadovaným textovým řetězcem. Tak například příkaz plot [0:1] x title ’prvni’, x**2 title ’druha’ vykreslí v intervalu [0, 1] grafy funkcí x a x2 opatřené popisky prvni a druha, jak vidíme na obrázku 6. Zobrazený výsledek však má nepříjemnou vadu. Legenda se bohužel vykreslila na tak nevhodném místě, že zasahuje do grafů. Naštěstí nám ale Gnuplot umožňuje její polohu nastavovat. V levém horním rohu máme docela dost místa, napíšeme tedy set key left top a při příštím vykreslení (například příkazem replot) se nám legenda poslušně přesune tam, kam jsme chtěli. Způsobů, jak lze manipulovat s podobou výsledného grafu, existuje nepřeberné množství a není v našich silách je zde všechny rozebrat. Věříme, že váš hlad po dalších informacích ukojí oficiální dokumentace. A nyní se již pusťme do dalšího důležitého tématu, kterým je výstup do souboru. 3.1.3
Možnosti výstupu
To, co dělá z Gnuplotu pro nás tak zajímavý nástroj, je vedle jeho dostupnosti a bohatých vykreslovacích možností jistě i schopnost poskytovat výstup v široké škále různých formátů od rastrových obrázků přes vektorové (postscript, PDF) až po zdrojové texty pro MetaPost či LATEX. Pro nás bude důležitý zejména výstup do formátu postscript, protože jej můžeme velmi snadno začlenit do TEXovských dokumentů. Zájemci o informace o ostatních formátech je mohou snadno získat z volně dostupné dokumentace ke Gnuplotu. Přejeme-li si tedy výstup v již zmíněném postscriptu, zajistíme to příkazem set term postscript. Jestliže na konec přidáme ještě slovo eps, bude použit takzvaný encapsulated (zapouzdřený) postscript. Když už máme nastaven formát výstupu, zbývá ještě určit název výstupního souboru. Zadáme tedy příkaz set output ’mujsoubor.eps’. A nyní je už vše připraveno k tomu, abychom do souboru mujsoubor.eps kreslili. To můžeme provádět zcela přirozeně příkazem plot. Drobnou nevýhodou tohoto postupu však je skutečnost, že výsledný graf nevidíme (nekreslí se na obrazovku, ale do souboru). Zde nám může pomoci příkaz replot. Ještě před nastavením výstupu si graf vykreslíme na obrazovku a upravíme si jej postupně k obrazu svému. Když už jsme spokojeni, nastavíme výše popsaným způsobem výstup a necháme vše vykreslit znovu příkazem replot. Jakmile máme v souboru to, co chceme, měli bychom jej opět odpojit. Jinak totiž zůstane otevřený 14
a pokusíme-li se jej někam přesunout, přejmenovat či vymazat, nepodaří se nám to. Napišme tedy příkaz unset output a vše bude v pořádku. Náš první postscriptový obrázek z Gnuplotu je hotov! 3.1.4
Příkaz fit
Další velmi užitečnou schopností Gnuplotu je možnost fitování (prokládání) souboru dat funkční závislostí. Jednoduše řečeno jde o nalezení takové funkce z určité dané třídy funkcí (např. všech lineárních funkcí, polynomů čtvrtého řádu, exponenciálních funkcí, apod.), která v jistém smyslu co nejlépe vystihuje naměřená data. Zajímají-li vás bližší detaily, pak vězte, že Gnuplot používá nelineární metodu nejmenších čtverců a optimální parametry požadované funkční závislosti hledá metodou postupných iterací (tedy opakovaného provádění jistých operací za účelem postupného zpřesňování dosaženého výsledku). Výhodou je možnost fitování dat prakticky libovolnou třídou funkcí, nevýhodou pak to, že v určitých případech Gnuplot optimální parametry nemusí vůbec najít, případně s nimi dokonce může „utéctÿ do oblasti zcela nesmyslných hodnot11 . Příkaz fit se používá následujícím způsobem. Jako první parametr zadáme typ fitované funkční závislosti – například a*x pro přímou úměrnost, a*x+b pro lineární funkci, a*x**2+b*x+c pro kvadratickou, a*exp(b*x) pro exponenciálu, atd. Symboly a, b, c, apod. (lze však použít téměř libovolný identifikátor kromě x) zde představují parametry funkce, které bude Gnuplot určovat. Jako druhý parametr zadáme název souboru s daty. Dále pak můžeme pomocí klíčového slova using stejně jako v případě příkazu plot specifikovat, které dva sloupce nás zajímají. Potom Gnuplot při fitování bere všechna data se stejnou váhou. Můžeme však také udat místou dvou sloupců tři, kde třetí obsahuje chyby závislé proměnné. Pak se přesnější hodnoty berou s větší váhou než ty, jejichž chyba je větší12 . Na závěr musíme Gnuplotu sdělit, které parametry má hledat, což uděláme uvedením slova via následovaného jejich seznamem. Po potvrzení příkazu ENTERem Gnuplot provede hromadu výpočtů a o svém snažení (a – ve velké většině případů – o jeho úspěchu) nám podá obsáhlou zprávu, kterou zároveň uloží do logovacího souboru (implicitně se jmenuje fit.log). Zde najdeme finální hodnoty hledaných parametrů včetně jejich chyb a korelační matice (tj. do jaké míry jsou spolu tyto hodnoty svázány). Hodnoty parametrů také zůstanou v paměti Gnuplotu uloženy až do jeho ukončení pod námi zvolenými identifikátory. Můžeme je tedy dále využít a například do grafu vykreslit spolu s daty ze souboru i jimi proloženou závislost. Demonstrujme si vše na příkladu. Zkusíme data z měření teploty během dne proložit polynomem čtvrtého řádu. Samotný fit provedeme příkazem fit a*x**4+b*x**3+c*x**2+d*x+e ’teplota.txt’ using 1:2:4 via a,b,c,d,e Následně vykreslíme do grafu naměřené teploty a pak do něj přidáme proloženou závislost. plot ’teplota.txt’ using 1:2:4 with errorbars replot a*x**4+b*x**3+c*x**2+d*x+e Grafy nám poněkud nešťastně zasahují do legendy a tak si ji příkazem set key left top přesuneme do levého horního rohu a vše vykreslíme znovu zadáním replot. Teď když už jsme s výsledkem spokojeni, si jej necháme vykreslit do postscriptového souboru: 11
I tuto situaci je ale někdy možno řešit – můžeme Gnuplot trochu „postrčitÿ žádoucím směrem tak, že mu „napovímeÿ, jaké by zhruba hledané parametry měly být. K tomu je ovšem musíme umět nějakým způsobem odhadnout. 12 Přesněji – váha každé hodnoty je rovna převrácené hodnotě kvadrátu její chyby.
15
30 ’teplota.txt’ using 1:2:3 a*x**4+b*x**3+c*x**2+d*x+e 28 26 24 22 20 18 16 14 12 10 8 0
5
10
15
20
25
Obrázek 7: Naměřená data spolu s proloženou závislostí set term postscript eps set output ’teplota.eps’ replot Při bližším pohledu na získaný graf v obrázku 7 si všimneme, že námi navržená polynomiální závislost zvláště v podvečerních hodinách neaproximuje naměřená data příliš dobře. Stálo by tedy za to vyzkoušet i jiné typy závislostí.
3.2
„Programyÿ v Gnuplotu
Pokud vás trápí, že ať už chcete v Gnuplotu udělat cokoliv, musíte mu to po každém spuštění vždy ručně zadat znovu, máme pro vás dobrou zprávu. Také v tomto případě si můžete napsat jednoduchý program, který pak necháte Gnuplotem zpracovat, kdy budete chtít. Docela dobře vám k tomu poslouží i obyčejný poznámkový blok. Jednoduše napíšete požadovanou posloupnost příkazů, textový soubor uložíte (nechť se nachází například v adresáři C:\mujadresar pod názvem program.txt). Pokud nyní z příkazového řádku zavoláte Gnuplot s parametrem obsahujícím cestu k danému programu, bude v něm obsažená posloupnost příkazů postupně vykonána, načež se Gnuplot automaticky ukončí. Jak se to celé provede? Předpokládejme, že máme Gnuplot nainstalován v adresáři C:\gnuplot. Spustíme si příkazový řádek a dostaneme se do adresáře C:\gnuplot\bin. Skončili jsme v podadresáři bin gnuplotího kořenového adresáře, protože právě zde se nachází vlastní program Gnuplotu s názvem wgnuplot.exe. Kdybychom nyní do příkazového řádku zadali wgnuplot, spustil by se Gnuplot ve windowsovském okně tak, jak jej známe. Protože však chceme pouze zpracovat náš programu, napíšeme do příkazového řádku wgnuplot C:\mujadresar\program.txt
16
a potvrdíme ENTERem. Možná se vám bude zdát, že se vůbec nic nestalo. Pokud ale zkontrolujete adresář C:\mujadresar, najdete tam výsledný obrázek, který vám Gnuplot nakreslil13 . Tolik k úvodní teorii, jistě jste už ale nedočkaví si to vše vyzkoušet v praxi. Zde je náš první příklad. Mějme v souboru teplota.txt zaznamenány výsledky měření teploty ve dvou dnech. V prvním sloupci čas v hodinách, v druhém a třetím teploty z obou dnů, ve čtvrtém a pátém jejich chyby (to vše ve stupních Celsia). Následující sekvence příkazů slouží k automatickému zpracování tohoto měření do formy grafu. Výsledkem bude graf závislosti teploty v kelvinech na čase v minutách, včetně příslušných popisků a legendy, ve formátu encapsulated postscript (EPS). set key left top set xlabel ’t [min]’ set ylabel ’T [K]’ plot ’teplota.txt’ using ($1*60):($2+273.15):4 with errorbars title ’den 1.’ replot ’teplota.txt’ using ($1*60):($3+273.15):5 with errorbars title ’den 2.’ set term postscript eps set output ’teplota.eps’ replot Smysl jednotlivých řádků14 je poměrně jednoduchý. První nastaví pozici legendy grafů do levého horního rohu. Druhý a třetí určují popisky os. Ve čtvrtém se vykreslí graf závislosti teplot z druhého sloupce (první den), v pátém je přidán obdobný graf pro teploty z třetího sloupce (druhý den). Příkazem na šestém řádku je určen formát výstupu – enhanced postscript. Sedmý udává název výstupního souboru. Poslední řádek pak zajistí překreslení výsledných grafů, tentokrát již do postscriptového souboru (nezapomenout! – bez posledních tří řádků bychom neměli žádný výstup). Přímo výstavní výsledek, který získáme zpracováním výše uvedeného programu, můžete obdivovat na obrázku 8. A na závěr ještě jeden příklad, tentokrát ze života – konkrétně z měření viskozity škrobové suspenze rotačním viskozimetrem ve fyzikálním praktiku. Vzhledem k tomu, že jde o takzvanou nenewtonovskou kapalinu, závisí její viskozita na rychlosti otáčení hřídele měřícího přístroje. Tato závislost je přibližně mocninná a exponent této mocniny lze určit měřením. Soubor skrob_visc_L3.txt obsahuje následující naměřená data 12 2850 20 2089 30 1684 50 1299 60 1172 100 922
100 60 40 24 21 12
První sloupec udává rychlost otáčení měřící hřídele v otáčkách za minutu, druhý pak naměřenou viskozitu v mPa · s a třetí chybu měření opět ve stejných jednotkách. Skript pro Gnuplot, který tato data úhledně zpracuje a poskytne nám graf i parametry nafitované mocninné závislosti, může vypadat třeba takto: 13
Pokud jste samozřejmě nezapomněli v programu zajistit výstup do souboru – implicitní výstup na obrazovku by vám nebyl příliš platný vzhledem k tomu, že se Gnuplot po vykonání příkazů opět ihned vypne. 14 Pozor – čtvrtý a pátý řádek se nám nevešly na stránku, takže jsou oba rozlomeny do dvou řádků. V programu však musejí být celistvé.
17
305 den 1. den 2.
300
T [K]
295
290
285
280 0
200
400
600
800 t [min]
1000
1200
1400
1600
Obrázek 8: Výsledek demonstračního gnuplotího programu set key box set key right top set key title ’Legenda’ set xlabel ’otacky [rpm]’ set ylabel ’viskozita [mPa.s]’ fit a*x**b ’skrob_visc_L3.txt’ using 1:2:3 via a,b plot [10:110] ’skrob_visc_L3.txt’ using 1:2:3 title ’namerena data’ with errorbars replot a*x**b title ’mocninna zavislost’ set term postscript mono eps set output ’graf.eps’ replot První řádek Gnuplotu říká, že kolem legendy chceme vykreslit rámeček, druhý ji pak umisťuje do pravého horního rohu a třetí nastavuje její nadpis na Legenda. Další dva řádky se starají o popisky os, šestý provádí fitování měřených dat funkcí axb . Sedmý vykreslí naměřená data a osmý nafitovanou mocninnou závislost. Poslední tři řádky už jen zajistí výstup do souboru graf.eps v černobílém zapouzdřeném postscriptu. Výsledek vidíme na obrázku 9. Charakter mocninné závislosti proložené získanými daty můžeme následně vyčíst z logovacího souboru. Jeho nejdůležitější část, která nás především zajímá, vypadá asi takto: Final set of parameters ======================= a b
= =
9879.48 -0.517108
Asymptotic Standard Error ========================== +/- 449.7 +/- 0.01125
18
(4.552%) (2.175%)
3500 Legenda namerena data mocninna zavislost 3000
viskozita [mPa.s]
2500
2000
1500
1000
500 20
40
60 otacky [rpm]
80
100
Obrázek 9: Závislosti viskozity na rychlosti otáčení hřídele
correlation matrix of the fit parameters:
a b
a 1.000 -0.989
b 1.000
Jsou zde uvedeny hodnoty parametrů a a b nafitované funkce včetně jejich směrodatných i relativních odchylek a korelační matice. Co víc bychom si mohli přát?
3.3 3.3.1
Pár tipů závěrem Jednodušší zpracování programů v Gnuplotu
Jste-li již maličko zběhlejší v práci s počítačem, poradíme vám trik, jak si výrazně zjednodušit zpracovávání gnuplotích programů. Zcela tak odpadne práce s příkazovým řádkem. Zvolte si nějakou souborovou příponu, kterou zatím žádný z programů na vašem počítači nepoužívá – například .gpl, .gpt atd. Tuto příponu pak asociujte s Gnuplotem – ve Windows XP se to dělá v Průzkumníku přes Nástroje – Možnosti složky a záložku typy souborů. Pak již jen stačí, když gnuplotí programy místo standardní přípony textových souborů .txt ukládáte pod vámi zvolenou příponou a zpracováváte v Průzkumníku (nebo třeba Total Commanderu, chcete-li) jediným dvojklikem. Vaše touha po zjednodušování může jít ještě dál a vytváření souborů s programy pro Gnuplot můžete zcela přenechat vašemu programu v (například) Pascalu, který se tak postará nejen o numerickou simulaci daného fyzikálního problému a o výstup dat do souboru, ale i o zavolání Gnuplotu a předání vytvořeného instrukčního souboru. Jestliže si tedy trochu pohrajete, můžete jediným spuštěním vašeho pascalovského programu získat cenná numerická data i jejich pěkné grafické zpracování. 19
3.3.2
PSfrag
Jestliže vytváříte grafy v Gnuplotu za účelem jejich umístění do LATEXovského dokumentu a jste-li k tomu navíc esteticky citlivější povahy, může vám vadit fakt, že Gnuplot nedisponuje příliš pěknými fonty. A vy byste přitom tak moc chtěli, aby t v popisku vodorovné osy grafu bylo stejné jako t ve vzorci y = A sin(ωt). Gnuplot sice umožňuje výstup grafů v TEXovském formátu, ale žel bohu jsou vykreslovací možnosti samotného TEXu poměrně chabé (vyzkoušejte si to a zděsíte se zubatých grafů a podobných ošklivostí). Tudy tedy cesta nepovede. Naštěstí však existuje pomoc ve formě přídavného balíčku pro TEXjménem PSfrag. Použijeme-li jej v našem dokumentu, prohledají se veškeré texty ve vkládaných postscriptových obrázcích a je-li nalezen text ve formátu \tex{cokoliv}, vypíše se na jeho místě cokoliv, tentokrát ale TEXovským písmem (je možné i vkládání matematických formulí). Graf tedy můžeme s použitím PSfragu opatřit dobře vypadajícími popisky os například takto. set xlabel ’\tex{$\tau$ [s]}’ set ylabel ’\tex{$t$ [$^{\circ}$C]}’
4
Numerické metody
Z předchozích kapitol jste si snad přinejmenším teoreticky osvojili práci v programech Freepascal a Gnuplot. Aby se vám však otevřely dveře k jejich plnohodnotnému fyzikálnímu využití, budete potřebovat znát i pár chytrých numerických metod, které umožňují řešit matematické problémy (a tedy i problémy fyzikální, zapsané matematickým jazykem) numericky. Nejdůležitější z těchto metod si nyní stručně popíšeme. Úvodem však ještě několik slov o tom, kde může zvídavý čtenář utišit svůj hlad po dalších vědomostech, jež mu neposkytne tento studijní text. O numerických výpočetních metodách velmi obsáhle pojednává například anglickojazyčný text [1], poněkud stručněji pak již poněkud letitá česká kniha [2]. Pevně doufáme, že vám následující kapitoly poskytnou alespoň určitý základ a umožní vám začít psát první zajímavé programy.
4.1
Numerické řešení rovnic
Nejdříve se podívejme na to, jak se dají numericky řešit obyčejné rovnice. Někdy se zkrátka dostaneme do situace, kdy bychom chtěli řešit rovnici, která analyticky řešitelná není (možná s výjimkou nějakých pokročilejších technik, jako je použití nekonečných řad). Potom nám může pomoci některá z numerických metod. 4.1.1
Metoda půlení intervalů
Jednou z nejjednodušších, ale přesto poměrně účinných, je metoda půlení intervalů. Využívá toho, že „rozumnáÿ funkce (konkrétně spojitá), která v bodě a nabývá záporné a v bodě b kladné hodnoty, musí někde mezi body a a b projít nulou. Jestliže tedy máme řešit rovnici f (x) = 0 (uvědomme si, že každou rovnici o jedné neznámé lze na tento tvar převést) a jsou dány dva body a, b takové, že f (a)f (b) < 0, nachází se mezi nimi zcela jistě alespoň jeden kořen naší rovnice.
20
Obrázek 10: Numerické řešení rovnic – metoda půlení intervalu Podíváme se, jaké znaménko má funkční hodnota v bodě c = a+b . Pokud jsme se zrovna netrefili 2 bodem c do kořenu, pak alespoň jeden z intervalů (a, c), (c, b) bude jistě obsahovat kořen (který z nich to je zjistíme opět podle znamének funkčních hodnot v krajních bodech). Opakováním tohoto postupu získáváme stále menší a menší intervaly obsahující hledaný kořen. Jakmile délka intervalu klesne pod hodnotu požadované chyby, můžeme výpočet zastavit a za kořen považovat střed posledního získaného intervalu (samozřejmě s danou chybou). Čtyři první kroky právě popsaného algoritmu půlení intervalů schématicky znázorňuje obrázek 10. Body x1 , x2 , . . . jsou krajní body postupně získávaných intervalů, černé kroužky označují aktuální interval obsahující kořen v každém ze čtyř kroků. Pro ještě větší názornost uveďme krátkou část zdrojového kódu v pascalu, provádějící uvedené kroky. Hledáme přitom nulový bod funkce f(x:real):real (tu si samozřejmě musíme někde předtím definovat). ... // a,b...meze intervalu uzavirajiciho koren (plati f(a)*f(b)<0) // d...pozadovana presnost while (b-a>d) do if f((a+b)/2)*f(a)<=0 then b:= (a+b)/2 else a:= (a+b)/2; ... 4.1.2
Metoda sečen
Další často používanou metodou numerického řešení rovnic je takzvaná metoda sečen. Podobně jako v případě metody půlení intervalů se každá další aproximace xn+2 hledaného kořene určuje na základě 21
Obrázek 11: Numerické řešení rovnic – metoda sečen předchozích dvou aproximací xn a xn+1 . Metoda sečen však navíc zohledňuje i velikosti funkčních hodnot v xn a xn+1 15 . Principu metody sečen lze dát velmi jednoduchý geometrický význam. Máme-li graf funkce f (x), v každém kroku vedeme přímku body (xn , f (xn )) a (xn+1 , f (xn+1 )) (tj. body na grafu s x-ovými souřadnicemi xn a xn+1 ) a její průsečík s osou x nám dá hodnotu xn+2 . Graficky je tento postup znázorněn v obrázku 11. Matematicky můžeme jeden krok metody sečen zapsat jako xn+2 =
xn+1 f (xn+1 ) − xn f (xn ) . f (xn+1 ) − f (xn )
Podobný demonstrační program jako v případě metody půlení intervalu by nyní vypadal takto. ... // a,b...prvni dve aproximace korene // d...pozadovana presnost while (abs(b-a)>d) do begin x:= (a*f(a)-b*f(b))/(f(a)-f(b)); a:= b; b:= x; end; ... Možná vás při pohledu na vztah, kterým získáváme xn+2 z xn+1 a xn , napadla otázka, jestli se nemůže stát něco ošklivého, pokud náhodou budou funkční hodnoty f (xn+1 ) a f (xn ) hodně blízké. 15
Pokud je například f (xn ) = −0.001 a f (xn+1 ) = 1, je jistě výhodné volit xn+2 v blízkosti xn .
22
Obrázek 12: Numerické řešení rovnic – metoda regula falsi
Obrázek 13: Patologický případ pro metodu regula falsi Bohužel může a také se to často stane. Hodnota další aproximace nám pak vyjde v absolutní hodnotě nesmyslně velká, metoda tedy „ustřelíÿ kamsi zcela mimo původně zamýšlenou oblast. Tato vlastnost je velikou Achillovou patou metody sečen. Na rozdíl od metody půlení intervalů si totiž nehlídá, aby funkční hodnoty posledních dvou aproximací měly vždy opačné znaménko (pak by totiž případ s „ustřelenímÿ vůbec nemohl nastat). 4.1.3
Metoda regula falsi
Z tohoto důvodu se častěji používá její drobná modifikace, takzvaná metoda regula falsi. Ta kombinuje dobré vlastnosti obou dosud zmíněných metod. Funguje v podstatě stejně jako metoda půlení intervalů, jen místo středu intervalu bereme za dělící bod hodnotu získanou ze vztahu pro metodu sečen. Graficky to můžeme vidět znázorněno na obrázku 12. Výhodou metody regula falsi je, že odstraňuje nepěkné chování metody sečen, na druhou stranu však může konvergovat k výsledku mnohem rychleji než metoda půlení intervalů. Protože ale nic není dokonalé, musíme říct, že pro některé funkce f (x) regula falsi konverguje pomaleji než půlení intervalů (které bychom mohli trefně vystihnout rčením „pomalu, ale jistěÿ). Příkladem může být funkce, jejíž graf ukazuje obrázek 13. 23
Obrázek 14: Numerické řešení rovnic – Newtonova metoda Z tohoto důvodu je výhodné vybavit algoritmus řešící rovnice touto metodou nějakým druhem kontrolního mechanismu, který v případě příliš pomalého zkracování intervalu provede zkusmo jeden či více dalších kroků metodou půlení, než opět přejde k metodě regula falsi. 4.1.4
Newtonova metoda
Další metodou z poněkud jiného soudku je metoda Newtonova. Ta již na rozdíl od předchozích tří metod nepotřebuje uchovávat dvě poslední aproximace, aby z nich určila následující, ale stačí jí pouze jedna. Pomocí hodnoty derivace funkce f v bodě xn se graf funkce nahradí v jeho okolí přímkou danou vztahem y = f (xn ) + (x − xn )f 0 (xn ) a průsečík s osou x určí hodnotu xn+1 (viz obrázek 13). Jednoduchou úpravou dostaneme vztah xn+1 = xn −
f (xn ) . f 0 (xn )
Vidíme, že Newtonova metoda trpí podobným neduhem jako metoda sečen. Pokud bude derivace f blízká nule, můžeme se dočkat ošklivého překvapení. Navíc je nutné mít funkci, která vrací derivaci funkce f v libovolném bodě. Tu můžeme do programu vložit buď ručně, nebo ji počítat numericky. Každopádně je ovšem nutné hlídat v každém kroku hodnotu f 0 (xn ) a hrozí-li komplikace, přejít k bezpečnější metodě (půlení nebo regula falsi). Výhodou Newtonovy metody naopak je, že dokáže v jistých situacích hledat nulové body funkcí i v celé komplexní rovině. Na závěr našeho bleskového výkladu metod numerického řešení rovnic uveďme jeden praktický fyzikální příklad. Jak je známo, je spektrální rozdělení intenzity záření černého tělesa (tj. výkon připadající na infinitezimální interval frekvencí dω dělený tímto dω) dáno Planckovým vyzařovacím
24
zákonem
Aω 3 , e¯hω/kT − 1 kde A je určitá konstanta úměrnosti (přímo úměrná velikosti povrchu tělesa), která nás však nebude příliš zajímat, h ¯ je Planckova konstanta a k Boltzmannova konstanta. Zajímalo by nás, na jaké frekvenci dané černé těleso o teplotě T září nejintenzivněji, jinými slovy pro jakou frekvenci ω je I(ω) maximální. Abychom to zjistili, vypočteme derivaci I 0 (ω) a položíme ji rovnu nule, tedy I(ω) =
I 0 (ω) =
Aω 3 h ¯ ¯hω/kT h ¯ ω ¯hω/kT 3Aω 2 − e = 0 ⇒ 3(e¯hω/kT − 1) = e . h ¯ ω/kT h ¯ ω/kT 2 e − 1 (e − 1) kT kT
Zavedeme-li si substituci x =
hω ¯ , kT
přejde tato rovnice na tvar 1 − e−x =
x . 3
Pokud nyní nalezneme řešení této rovnice (samozřejmě různé od nuly), můžeme pro hledanou frekvenci napsat vztah xk T. ω= h ¯ Řešení budeme hledat metodou regula falsi. Zavedeme si funkci f (x) = e−x + x3 − 1, jejíž nulový bod chceme najít. Snadno se přesvědčíme, že f (1) < 0 a f (3) > 0, hledaný kořen tedy leží v intervalu (1, 3). Náš program bude vypadat takto
... a:= 1; // pocatecni nastaveni mezi intervalu b:= 3; N:= 0; while (abs(b-a)>1E-12) do // smycka vypoctu begin N:= N+1; fa:= f(a); fb:= f(b); writeln(’a = ’,a,’; b = ’,b); // vypis stavu vypoctu writeln(’f(a) = ’,fa,’; f(b) = ’,fb); if (fa=0) and (fb=0) then begin b:= a; break; end; c:= (a*fb-b*fa)/(fb-fa); // vypocet noveho deliciho bodu... fc:= f(c); // ...a funkcni hodnoty v~nem if fc=0 then begin a:= c; b:= c; break; end; if ((fa*fc<0) and ((b-a)>10*(b-c))) or // pokud by bylo zmenseni ((fa*fc>0) and ((b-a)>10*(c-a))) // intervalu moc male, then begin c:= (a+b)/2; fc:= f(c); end; // pouzijeme radeji puleni if fa*fc<0 // zmenseni intervalu then b:= c else a:= c; end; writeln(’x = ’,(a+b)/2); // vypis vysledku writeln(’chyba = ’,(b-a)/2); 25
writeln(N,’ kroku’); readln; ... Posledních několik řádků z mnoha, které nám program vypíše, bude vypadat asi takto ... f(a) = -4.28259858131774E-0017; f(b) = 9.28232100558257E-0013 a = 2.82143937212208E+0000; b = 2.82143937212377E+0000 f(a) = -4.28259858131774E-0017; f(b) = 4.64094583076113E-0013 x = 2.82143937212250E+0000 chyba = 4.23661106196960E-0013 42 kroku . Zjistili jsme tedy, že x = 2,82143937 . . .. Dosazením pak dostaneme vztah . ω = 3,6939 · 1011 s−1 · K−1 · T, neboli zapsáno pomocí obyčejné frekvence f GHz . · T. f = 58,790 K
4.2
Minimalizace a maximalizace funkcí
Další úlohou, před kterou jsme často postaveni, je najít minima, resp. maxima určité dané funkce. 4.2.1
Metoda zlatého řezu
Jednoduchou avšak poměrně spolehlivou metodou je takzvaná metoda zlatého řezu, která je v jistém smyslu analogií metody půlení intervalů z numerického řešení rovnic. Zatímco půlení intervalů využívá „uvězněníÿ kořene mezi dvěma body, metoda zlatého řezu uzavírá minimum resp. maximum trojicí bodů. Pokud totiž máme tři body a > b > c a přitom f (b) ≤ f (a) a f (b) ≤ f (c), funkce f (o které budeme opět předpokládat, že je spojitá) musí někde na intervalu (a, c) nabývat lokálního minima (podobně pro maximum). Pokud tedy na začátku výpočtu máme dány takové tři body, můžeme dále postupovat podobně jako při půlení intervalu. Zvolíme si nový bod x například mezi body b a c. Pokud je f (x) > f (b), bude naší novou trojicí bodů a, b, x. Pokud naopak f (x) ≤ f (b), bude jí b, x, c. Zbývá jen určit, jakým způsobem budeme volit v každém kroku dělící bod x. Ukazuje se, že z jistých důvodů je optimální, když x vždy volíme ve větším z intervalů (a, b), (b, c) a to tak, že jej tento nový bod rozdělí√v poměru zlatého řezu (přičemž bod x bude blíže bodu b), tj. v poměru γ : (1 − γ), kde γ = (3 − 5)/2. Rozdíl oproti numerickému hledání kořenů je v tom, jaké přesnosti můžeme numerickými metodami dosáhnout. Zatímco při hledání kořenů je tato přesnost řádově srovnatelná s√přesností použité počítačové aritmetiky a datových typů, při určování extrémů je to jen přibližně ε (tedy používá-li náš reálný typ 16 platných číslic, získáme polohu extrému maximálně s přesností asi na 8 číslic). Program pro hledání minima funkce f metodou zlatého řezu tedy může vypadat například takto ... a:= -1; b:= 0;
// nastaveni pocatecni trojice vymezujici // minimum na -1,0,1 26
c:= 1; g:= 0.382; // konstanta pro zlaty rez while c-a>1E-6 do // dokud nedosahneme dostatecne presnosti begin R:= (a+c>2*b); // do promenne R (boolean) ulozme, zda bude //delici bod x v~intervalu (b,c) if R // podle toho jej nastavime then x:= (1-g)*b + g*c else x:= (1-g)*b + g*a; fb:= f(b); // fx:= f(x); if (fx<=fb) // then begin if R then a:= b else c:= b; b:= x; end else // pripad if R then c:= x else a:= x; end;
vypocet funkcnich hodnot v~b a~x pripad f(x)<=f(b)
// nove nastaveni trojice bodu
f(x)>f(b) // nove nastaveni trojice bodu
... 4.2.2
Metoda konjugovaného gradientu
Další metodou, kterou lze využít i v případě funkcí více proměnných, je takzvaná metoda konjugovaného gradientu. Jde o vylepšení jednoduché metody, která by vás možná napadla, totiž vydat se z daného bodu směrem nejrychlejšího klesání, zastavit se, až narazíme na minimum (vzhledem k přímce, po které se pohybujeme) a pak krok libovolněkrát opakovat (této metodě se říká metoda nejrychlejšího sestupu). Uvedený algoritmus je poměrně neobratný, jelikož má tendenci v „úzkých údolíchÿ postupovat po malých krůčcích. Tuto nepříjemnost odstraňuje právě sofistikovanější metoda konjugovaného gradientu. Ozřejmění principu celého algoritmu by vyžadovalo určité znalosti lineární algebry, jistě nám tedy odpustíte, pokud v zájmu zachování stručnosti tohoto studijního textu přejdeme přímo ke kuchařce (zájemci mohou nahlédnout do Numerických receptů [1]). Jednotlivé iterace při hledání minima funkce f (~x) budou sestávat z těchto kroků: • Z bodu p~i vyjdeme ve směru ~hi a zastavíme se v minimu, které označíme p~i+1 . • Vektor ~gi+1 položíme roven mínus gradientu f v bodě p~i+1 . ·~gi+1 ~ hi (tzv. Fletcherova-Reevesova metoda) nebo • Vypočítáme ~hi+1 jako ~gi+1 + ~gi+1 ~gi ·~gi
~gi+1 +
(~gi+1 − ~gi ) · ~gi+1 ~ hi ~gi · ~gi 27
(Polak-Ribiere). Počáteční hodnota p~1 může být zvolena prakticky libovolně, ~g1 a ~h1 položíme rovny mínus gradientu v bodě p~1 Jak bude vypadat implementace tohoto postupu v programu? Pro jednoduchost mějme funkci pouze dvou proměnných, zobecnění na případ více proměnných je zřejmé. Máme tedy f(p:TVec):real, kde TVec je námi nově definovaný datový typ představující dvojrozměrný vektor. type TVec = record x,y: real; end; Dále budeme potřebovat umět vypočítat gradient této funkce (tj. vektor o složkách ∂f /∂x a ∂f /∂y). Dejme tomu, že jej počítá funkce16 gradf(x,y:real):TVec Dále budeme mít proměnné p, h, g, gs typu TVec. Proměnná p bude představovat aktuální polohu v prostoru proměnných, zbývající dvě jsou pomocné (proměnné g a gs uchovávají poslední dvě hodnoty ~g ). Jeden krok pak bude vypadat takto p:= LineMin(p,h); gs:= g; g:= MultVec(gradf(p),-1); gamma:= MultVecs(g,g)/MultVecs(gs,gs); h:= AddVecs(g,MultVec(h,gamma)); Tento kousek kódu sice vypadá stručně, ale obsahuje několik neznámých funkcí, a sice LineMin, MultVecs, MultVec a AddVecs. Funkce LineMin(p:TVec; h:TVec):TVec zajišťuje hledání minima ve směru daném vektorem ~h počínaje bodem p~. To můžeme udělat napřílad výše popsanou metodou zlatého řezu. Jediný problém je v tom, že na začátku nemáme minimum „chycenoÿ v nějakém intervalu. Funkce LineMin si tedy musí nejdříve zkusmo popojít ve směru ~h a najít trojici bodů ta < tb < tc takových, že f (~p + t1~h) > f (~p + t2~h) a f (~p + t3~h) > f (~p + t2~h). Funkce MultVecs(a,b:TVec):real, MultVec(a:TVec; k:real):TVec a AddVecs(a,b:TVec):TVec zajišťují po řadě operace skalárního násobení dvou vektorů, násobení vektoru skalárem a součtu dvou vektorů. Popsaný algoritmus samozřejmě není ideální, lze na něm mnohé vylepšit, obzvláště pak ve funkci hledající minimum na přímce. Pro ilustraci však jistě stačí. Hotový demonstrační program conjgrad si můžete stáhnout na webu FYKOSu. 4.2.3
Simulované žíhání
Další velmi zajímavou metodou hledání extrémů funkcí (přinejmenším z fyzikálního hlediska) je takzvané simulované žíhání. Tato metoda se hodí obzvláště v případech, kdy má minimalizovaná funkce diskrétní (a přitom velmi obsáhlý) definiční obor. Její princip je inspirován statistickou fyzikou. Pokud kov zahřátý na vysokou teplotu ochlazujeme dostatečně zvolna, dostaneme na konci velmi pravidelnou strukturu krystalické mřížky, která minimalizuje její energii. 16
Buď pomocí explicitního vztahu, nebo numericky.
28
Obrázek 15: Provedení náhodné fluktuace při řešení TSP simulovaným žíháním Můžeme si tedy představovat, že zadaná funkce, kterou máme minimalizovat, představuje energii nějakého systému. Zavedeme si „teplotuÿ T tohoto systému, která bude ovliňovat intenzitu simulovaných tepelných fluktuací v něm, a tu budeme v průběhu výpočtu snižovat. Na začátku výpočtu, kdy je teplota vysoká, se systém může dostat prakticky do všech přípustných stavů. Se snižující se teplotou však jeho ochota podnikat kroky do energeticky vyšších stavů klesá a začíná se zdržovat spíše v oblastech s nižší energií. Nakonec pro velmi nízké teploty systém „zamrzneÿ v minimu energie. Obecně je toto minimum pouze lokální, ale díky pomalému ochlazování máme docela dobrou šanci, že nalezneme globální minimum. To je výhoda simulovaného žíhání proti primitivnějším algoritmům, které pouze mechanicky sledují směr nejrychlejšího klesání uvažované funkce a jsou tak poměrně náchylné k uváznutí v mělkém lokálním minimu. Tolik tedy k teorii a teď už konečně k tomu, jak simulované žíhání prakticky realizovat. Vše si popíšeme na konkrétním případu, a sice na slavném problému obchodního cestujícího (travelling salesman problem – TSP). Ten má uskutečnit okružní pracovní cestu po N daných městech a otázkou je, jak zvolit jejich pořadí, aby nacestovaná vzdálenost byla co nejmenší. O tomto problému je známo, že při jeho řešení nemáme prakticky jinou možnost než vyzkoušet všechny možné trasy17 . Při větším počtu měst je tedy hledání přesného řešení nad síly dnešních počítačů. Naštěstí se však ukazuje, že některé algoritmy dokáží tento problém řešit s poměrně dobrými výsledky přibližně a hlavně v rozumném čase. Jedním z takových algoritmů je právě simulované žíhání. Než začneme psát vlastní program, musíme si rozmyslet, co bude naším konfiguračním prostorem (tj. množinou možných stavů). V případě TSP to bude zřejmě množina všech okružních cest po N městech. Dále musíme zvolit „energiiÿ, kterou budeme minimalizovat. Nejjednodušší volbou bude položit ji rovnu celkové délce cesty. Ještě si musíme ujasnit, jak budeme provádět v našem systému „fluktuaceÿ. Potřebujeme v každém kroku provádět malé změny konfigurace. Přirozenou volbou (použitou například v [1]) je vybrat náhodně jistý úsek cesty a změnit směr, ve kterém jím cestující projde (viz obrázek 15). Když takovou změnu provedeme, vypočítáme novou energii (délku), zjistíme její změnu ∆E oproti energii v předchozím stavu a vypočítáme pravděpodobnost P , s jakou tuto fluktuaci přijmeme (jinak se vrátíme k předchozímu stavu) jako ∆E . P = exp − T 17
Samozřejmě si můžeme pomoci nějakými chytrými úvahami, kterými předem vyloučíme zjevně nesmyslné možnosti, ale ani tak svou situaci moc nezlepšíme.
29
Obrázek 16: Srovnání řešení TSP nalezených simulovaným žíháním (vlevo) a „metodou slepé honby za minimemÿ (vpravo) Je-li P ≥ 1, změnu ponecháme vždy. V tomto vztahu T značí naši uměle zavedenou „teplotuÿ. Popsaný krok budeme provádět mnohokrát za sebou a přitom budeme postupně snižovat teplotu, dokud se nestanou další změny v hodnotě energie natolik malými či málo pravděpodobnými, že by další pokračování výpočtu již nemělo smysl. V předešlé větě, ačkoliv se může zdát zcela zřejmá a jasná, se ve skutečnosti skrývá největší problém při používání simulovaného žíhaní. Určení optimálního průběhu snižování teploty (jejíž počáteční hodnota by měla být řádově rovna maximálním změnám energie, které lze v systému při jednom kroku očekávat) a celkového počtu kroků většinou vyžaduje jisté experimentování. Pokud totiž není průběh žíhání volen vhodně, můžeme touto metodou dostat dokonce horší výsledky než by nám poskytla primitivní strategie vždy zvolit krok snižující energii. Program zihaniTSP řešící problém obchodního cestujícího je pro vaši inspiraci k dispozici ke stažení na FYKOSím webu. S jeho pomocí jsme provedli pro několik náhodně vygenerovaných rozmístění padesáti měst srovnání výsledků simulovaného žíhání s primitivní metodou podnikající vždy kroky k nižším hodnotám energie (délky). Výsledné cesty včetně jejich délek l můžete pro dvě konkrétní rozmístění měst vidět na obrázku 16. Metoda simulovaného žíhání tedy při vhodném nastavení nachází i řešení, která bychom při pouhém neustálém snižování délky nenalezli jednoduše proto, že by algoritmus uvízl v lokálním minimu. Rozdíl v délkách cest získaných oběma metodami by se sice mohl zdát zanedbatelný, ale relativně může představovat až několikanásobné přiblížení k optimálnímu řešení. Netvrdíme samozřejmě, že je simulované žíhání všemocné, ale pro určitou třídu problému je jistě 30
velmi zajímavou a přitom jednoduchou alternativou deterministických algoritmů.
4.3
Numerická derivace
Často se ocitneme v situaci, kdy máme možnost vypočítat hodnotu dané funkce a chtěli bychom také hodnotu její derivace v určitém bodě. Přitom ji ale nedokážeme derivovat analyticky. Potom přichází ke slovu derivování numerické. Jeho princip vychází z definice derivace jakožto limity f (x + h) − f (x) , h→0 h
f 0 (x) = lim
kterou však z jistých důvodů nahradíme následující (s předchozí ekvivalentní v případě, že derivace existuje) f (x + h/2) − f (x − h/2) f 0 (x) = lim . h→0 h Důvodem tohoto zdánlivě bezdůvodného obratu je fakt, že je tento tvar symetrický a tedy se s ním v určitých případech lépe pracuje. Místo matematicky abstraktního limitního přechodu h → 0 si zvolíme konkrétní dostatečně malé18 h. Máme-li tedy v Pascalu definovanou funkci mojefce(x:real):real, můžeme její derivaci velmi jednoduše definovat například takto. function dmojefce(x:real):real; const h = 0.001; begin dmojefce:= (mojefce(x+h/2) - mojefce(x-h/2))/h; end; Volba konstanty h hodně závisí na charakteru derivované funkce. Je-li jí například periodická funkce s periodou 10−6 , pak volba h = 0.001 jistě nebude vhodná a bude třeba sestoupit ještě minimálně o čtyři (lépe však o více) řády níže. Tato metoda samozřejmě funguje pouze u funkcí, které jsou derivovatelné, jinak nám budou vycházet nesmysly a při zmenšování hodnoty h bude výsledek místo konvergence nejspíš divoce poskakovat. Máme-li zálusk také na druhou či ještě vyšší derivaci, lze tuto metodu také uplatnit. Uvědomímeli si, že druhá derivace je vlastně první derivace z první derivace atd., snadno uvěříme následujícím vztahům (platným v případě, že příslušná derivace v daném bodě existuje) f (x + h) − 2f (x) + f (x − h) , h→0 h2 f (x + 3h/2) − 3f (x + h/2) + 3f (x − h/2) − f (x − 3h/2) f 000 (x) = lim , 3 h→0 h Pn i n f (x − (2i − n)h/2) i=0 (−1) (n) i f (x) = lim . h→0 hn f 00 (x) = lim
18
Ve skutečnosti to ovšem nemusí být až tak jednoduché. Hodnota h musí být dostatečně malá, aby v rámci změny závislé proměnné o h bylo možno f považovat za lineární. Bude-li ale h příliš malé, mohou nám výpočty zhatit zaokrouhlovací chyby, které se objeví vždy když od sebe odečítáme hodně blízké hodnoty.
31
Obrázek 17: Numerické integrování – obdélníková metoda
4.4
Numerické integrování
Stejně tak jako lze numericky derivovat, můžeme na počítači provádět i operaci opačnou, čili integrovat. Numerické integrování je ještě mnohem mocnějším nástrojem než numerická derivace, neboť zderivovat danou funkci lze ve velké většině případů analyticky, zatímco při integrování velmi často narazíme na funkce, jejichž integrál není vyjádřitelný pomocí elementárních matematických funkcí (a někdy ani pomocí těch „divočejšíchÿ, jako jsou Besselovy, Neumannovy, Airyho či jiné funkce). Jak tedy postupovat, máme-li dánu funkci f (x) a chtěli bychom znát hodnotu jejího určitého integrálu od a do b? Napoví nám klasická představa určitého integrálu jako plochy omezené grafem funkce. Nejjednodušeji můžeme tuto plochu aproximovat velkým počtem vertikálních proužků o stejné šířce a o výšce dané hodnotou funkce v příslušném intervalu (je poměrně jedno, zda jde o hodnotu průměrnou, či nějak libovolně zvolenou); přitom připouštíme i zápornou výšku. Čím více obdélníků použijeme, tím přesnější bude výsledek. Tuto takzvanou obdélníkovou metodu vystihuje následující vztah19 Z
b
f (x) dx = lim a
N →∞
N X b−a
N
i=1
f (a + i(b − a)/N ).
Pokrývání plochy pod grafem obdélníky je zjevně dosti hrubé a asi vás napadá, jestli by se to nedalo provést lépe. Můžeme samozřejmě použít místo obdélníků lichoběžníky a přesnost integrace tak poněkud vylepšíme. Není těžké si odvodit, jak vypadá vzorec pro lichoběžníkovou metodu. ! Z b N −1 b − a f (a) + f (b) X + f (a + i(b − a)/N ) . f (x) dx = lim N →∞ N 2 a i=1 Existují samozřejmě ještě přesnější metody numerického integrování. Zájemce odkazujeme na již zmíněné zdroje: volně dostupný elektronický text [1] či odbornou literaturu [2]. 19
Tento vztah není jediný možný. V tomto případě bereme funkční hodnoty vždy na konci intervalu – pokud bychom však sčítali od i = 0 do i = N − 1, budou to naopak hodnoty na začátku intervalu a metoda bude fungovat stejně dobře. Podobně by posloužil třeba i střed intervalu.
32
Obrázek 18: Numerické integrování – lichoběžníková metoda Následující příklad ilustruje, jak by mohla vypadat funkce, která zintegruje předem danou funkci mojefce(x:real):real v zadaných mezích. function intmojefce(a,b:real; N:integer):real; var i: integer; s: real; begin s:= (mojefce(a)+mojefce(b))/2; for i:= 1 to N-1 do s:= s+mojefce(a+i*(b-a)/N); intmojefce:= s*(b-a)/N; end;
4.5
Řešení ODR Eulerovou metodou
Než přejdeme přímo k jádru tématu, jímž bude numerické řešení obyčejných diferenciálních rovnic, uvědomme si jejich souvislost s integrováním. Řešením jednoduché rovnice y 0 (x) = f (x), y(a) = 0 je totiž právě určitý integrál Z
x
y(x) =
f (t) dt. a
Jak bychom ale tuto rovnici řešili bez přímého použití numerického integrování? Mohli bychom například derivaci nahradit přibližným vzorcem a dostat tak y(x + h) − y(x) = f (x) ⇒ y(x + h) = y(x) + hf (x) ⇒ h ⇒ y(a + kh) = h(f (a) + f (a + h) + f (a + 2h) + . . . + f (a + (k − 1)h)). (1) 33
Volbou h = (b − a)/k pak dostaneme vzorec pro obdélníkovou metodu. Výhodou tohoto postupu je jeho schopnost poskytnout řešení rovnice i v případě, že je její pravá strana závislá kromě x i na y, pak už ovšem řešení obecně nemá tvar určitého integrálu. Nahrazením derivace její diskrétní aproximací získáme vzorec, který svazuje hodnotu y(x) s y(x + h). Jeho opakovaným použitím pak po krocích o velikosti h dokráčíme až k té hodnotě nezávislé proměnné, pro niž chceme nalézt hodnotu y. Tak například máme-li rovnici y 0 (x) = sin y(x), y(0) = 1, dostaneme výše popsaným způsobem vztah y(x + h) = y(x) + h sin y(x). Zvolíme-li si tedy například h = 0,01, můžeme přibližně určit hodnotu y(1) tak, že odvozený rekurentní vztah aplikujeme stokrát za sebou, což pro počítač není žádný problém. ... y:= 1; // do y vlozime pocatecni hodnotu y(0) for i:= 1 to 100 do y:= y+0.01*sin(y); // a~uz tam mame y(1) ... Takto si můžeme po libosti řešit diferenciální rovnice tvaru y 0 (x) = f (x, y(x)), kde f je nějaká daná funkce dvou proměnných. Co ale diferenciální rovnice vyššího řádu? Ty se přece jen ve fyzice vyskytují mnohem častěji. Budeme pro ně muset vyvíjet nějakou novou metodu? Naštěstí ne! Můžeme totiž použít následující trik. Celá až dosud popsaná teorie a postupy pro rovnice prvního řádu jsou stejně dobře použitelné i pro soustavy rovnic. Jinými slovy – chceme-li hledat funkce y1 (x), . . . , yk (x), splňující rovnice y10 (x) = f1 (x, y1 (x), y2 (x), . . . , yk (x)), .. . yk0 (x) = fk (x, y1 (x), y2 (x), . . . , yk (x)), uděláme to úplně stejně jako v případě rovnice jediné – nahrazením derivace přibližným vzorcem. Zavedeme-li si substituci y(x) = y1 (x), y 0 (x) = y2 (x), . . . , y (n−1) (x) = yn (x), pak můžeme rovnici n-tého řádu y (n) (x) = f (x, y(x), y 0 (x), . . . , y (n−1) (x)) přepsat jako soustavu rovnic prvního řádu. y10 (x) = y2 (x), y20 (x) = y3 (x), .. . 0 yn−1 (x) = yn (x), yn0 (x) = f (x, y1 (x), y2 (x), . . . , yn (x)) 34
Libovolnou diferenciální rovnici, z níž umíme vyjádřit nejvyšší derivaci v závislosti na nižších, tedy dokážeme numericky vyřešit! Uvažme například rovnici y 00 (x) = (y 0 (x))2 + y x+1 (x), y(0) = 0, y 0 (0) = 1. Tu si výše popsaným postupem převedeme na soustavu y10 (x) = y2 (x), y20 (x) = y22 (x) + y1x+1 (x), kterou přibližně vyřešíme použitím rekurentních vztahů y1 (x + h) = y1 (x) + hy2 (x), y2 (x + h) = y2 (x) + h(y22 (x) + y1x+1 (x)). Následující program (pro stručnost uvádíme jen tělo programu) vypíše tabulku závislosti y (tj. y1 ) na x od 0 do 1. ... x:= 0; y1:= 0; y2:= 1; for i:= 1 to 100 do begin writeln(x,’ ’,y1); y1n:= y1+0.01*y2; y2n:= y2+0.01*(sqr(y2)+power(y1,x+1)); y1:= y1n; y2:= y2n; x:= x+0.01; end; ... Na závěr si ještě povězme, proč je právě popsané metoda řešení diferenciálních rovnic (běžně zvaná metodou Eulerovou) použitelná jen v nejvyšší nouzi. Jde totiž o tzv. metodu prvního řádu; nebudeme se zde pouštět do matematických podrobností – zjednodušeně řečeno se při snížení kroku na polovinu relativní chyba metody sníží také na polovinu (což není nic moc, uvědomíme-li si, že některé jiné metody při stejném zmenšení kroku sníží svou chybu 16krát či 32krát). Dále tedy popíšeme poněkud použitelnější metody.
4.6
Runge-Kuttovy metody řešení ODR
Jedním z faktorů snižujících použitelnost Eulerovy metody je to, že využívá pro aproximaci derivace nesymetrický vztah. Pokud bychom jej chtěli nahradit výše uvedeným symetrickým vztahem, pak bychom pro získání y(x+h) z y(x) potřebovali hodnotu derivace v bodě x+h/2, kterou ovšem – pokud pravá strana diferenciální rovnice závisí na y – neznáme. Můžeme si však pomoci úvodním krokem z y(x) k y(x + h/2) pomocí asymetrického vztahu pro derivaci. Takto přibližně určíme y(x + h/2) a dosazením do pravé strany rovnice zjistíme y 0 (x + h/2), kterou užijeme pro výpočet y(x + h). Princip této tzv. Runge-Kuttovy metody druhého řádu pro rovnici y 0 (x) = f (y(x), x) ukazuje následující schéma. k1 = hf (y(x), x), k2 = hf (y(x) + k1 /2, x + h/2), y(x + h) = y(x) + k2 . 35
Zobecnění této metody na soustavu rovnic prvního řádu je poměrně zřejmé. Matematickými metodami (konkrétně rozvojem do Taylorovy řady) lze ukázat, že tato metoda je již na rozdíl od Eulerovy metodou druhého řádu. Podobně je možné nalézt i Runge-Kuttovy metody vyšších řádů, z nichž pravděpodobně nejpoužívanější (a pro naše účely většinou bohatě postačující) je Runge-Kuttova metoda čtvrtého řádu, reprezentovaná schématem. k1 = hf (y(x), x), k2 = hf (y(x) + k1 /2, x + h/2), k3 = hf (y(x) + k2 /2, x + h/2), k4 = hf (y(x) + k3 , x + h), y(x + h) = y(x) + k1 /6 + k2 /3 + k3 /3 + k4 /6. Ukažme si na příkladu, jak Runge-Kuttovou metodou čtvrtého řádu numericky vyřešit rovnici y (x) = y 2 (x) + x, y(0) = 1 pro x od 0 do 10: 0
... y:= 1; x:= 0; while x<=10 do begin writeln(x,’ ’,y); k1:= h*(sqr(y)+x); k2:= h*(sqr(y+k1/2)+x+h/2); k3:= h*(sqr(y+k2/2)+x+h/2); k4:= h*(sqr(y+k3)+x+h); y:= y+k1/6+k2/3+k3/3+k4/6; x:= x+h; end; ...
4.7
Adaptivní délka kroku
Nevýhodou dosud popsaných metod numerického řešení ODR je to, že si musíme předem pevně zvolit délku kroku. Přitom některé úseky hledaného řešení mohou být „divočejšíÿ a vyžadovat tedy kratší krok, který ale naopak bude zbytečně krátký pro jiné úseky. Bylo by tedy dobré nějakým způsobem umět délku kroku v průběhu výpočtu kontrolovat. To lze nejjednodušeji realizovat tak, že průběžně monitorujeme velikost chyby, jíž se dopouštíme, a podle ní pak regulujeme velikost kroku (je-li chyba příliš velká, musíme krok zmenšit; v opačném případě si můžeme dovolit krok zvětšit a ušetřit tak cenný výpočetní čas). Jakým způsobem ale lze sledovat velikost chyby? Můžeme například příslušný krok udělat dvakrát pomocí dvou různých metod, z nichž jedna bude řádově přesnější než druhá. Rozdíl výsledků nám dá přibližnou chybu kroku. Pokud budeme s dosaženou přesností spokojeni, přijmeme výsledek přesnější z obou metod a pokračujeme dále – s případným zvětšením kroku, je-li přesnost zbytečně příliš vysoká. Jestliže se nám chyba nelíbí, výpočet zopakujeme s menším krokem. Nejjednodušší způsob, jak právě popsaný princip aplikovat, je udělat každý krok dvakrát – jako jediný krok velikosti h a jako dva kroky velikosti h/2 (to bude naše „řádově přesnějšíÿ metoda). Další možností je využít chytrou metodu zvanou anglicky „embedded Runge-Kutta methodÿ (vložená Runge-Kuttova metoda). Ta je založena na skutečnosti, že můžeme poměrně malým počtem 36
Tabulka 2: Konstanty pro Cash-Karpovu ERK metodu i 1
ai
2
1 5
1 5
3
3 10
3 40
9 40
4
3 5
3 10
9 − 10
6 5
5
1
− 11 54
5 2
− 70 27
35 27
6
7 8
1631 55296
175 512
575 13824
44275 110592
253 4096
1
2
3
4
5
j=
bij
ci 37 378
c∗i 2825 27648
0
0
250 621
18575 48384
125 594
13525 55296
0
277 14336
512 1771
1 4
vyčíslení funkce f (y(x), x) na pravé straně diferenciální rovnice získat několik hodnot, z nichž pomocí speciálních vzorců „vykouzlímeÿ velikosti kroku y(x + h) − y(x) pro metody různých řádů (např. pátého a čtvrtého). Ty pak můžeme užít pro odhad chyby. Uveďme zde jednu z těchto metod, objevenou Cashem a Karpem, převzatou sem (jako ostatně i mnoho dalších informací) z [1]. Pro naše účely ji můžeme považovat za „spadlou z nebeÿ. k1 = hf (y(x), x), k2 = hf (y(x) + b21 k1 , x + a2 h), .. . k6 = hf (y(x) + b61 k1 + . . . + b65 k5 , x + a6 h), y(x + h) = y(x) + c1 k1 + . . . + c6 k6 , y ∗ (x + h) = y(x) + c∗1 k1 + . . . + c∗6 k6 . Zde y(x + h) je hodnota podle metody pátého řádu a y ∗ (x + h) podle metody čtvrtého řádu. Hodnoty použitých konstant uvádí tabulka 2 (také převzata z [1]). Odhad chyby kroku tedy udává výraz |y(x + h) − y ∗ (x + h)| = |(c1 − c∗1 )k1 + . . . + (c6 − c∗6 )k6 |. Při upravování délky použitého kroku většinou postupujeme tak, že si na začátku výpočtu pevně stanovíme hodnotu chyby ∆0 a pak předpokládáme, že chyba je přímo úměrná h5 (jde o metodu pátého řádu), tedy zjistíme-li po provedení kroku odchylku y(x + h) od y ∗ (x + h) rovnou ∆, změníme krok h na 0,2 ∆0 h, ∆ abychom chybu co nejvíce přiblížili ∆0 . Přitom bylo-li před touto změnou ∆ > k∆0 , kde k > 1 je vhodné malé číslo20 , vrátíme se a poslední krok provedeme znovu. 20
Důvod, proč je v této podmínce přítomen „bezpečnostníÿ koeficient k, je následující nepříjemnost. Jestliže bude
37
Pokud neřešíme jedinou rovnici prvního řádu, nýbrž jejich soustavu (což se nám stane velmi často), pak můžeme za chybu kroku považovat například největší z chyb zjištěných výše uvedeným způsobem pro jednotlivé rovnice. Další velmi užitečné informace, ať už o tomto konkrétním tématu, odlišných metodách řešení ODR, či jiných problémech numerické matematiky se zaměřením na programovací stránku věci naleznete na stránkách [1].
5 5.1
Bonbónek: LATEX a MetaPost LATEX
LATEXje značkovací jazyk pro sazbu dokumentů programem TEX. Píše se v něm velká spousta vědeckých článků. Byl napsán v osmdesátých letech Lesliem Lamportem. Jeho oficiální stránka je http://www.latex-project.org/. TEXje volně šiřitelný počítačový program pro vysoce kvalitní elektronickou sazbu. Je poměrně velmi složitý, proto je dobré začínat s LATEXem. Chcete-li elegantně psát svá řešení a nezůstávat otroky Wordu, doporučujeme navštívit stránky československého sdružení uživatelů TEXu http://www.cstug.cz. Naleznete zde především velmi dobrou literaturu. Přijde mi vhodné začít s četbou Nepříliš stručného úvod do LATEX2e, Neboli LATEX2e v 73 minutách. Naleznete jej v sekci „Přehled dokumentace pro TEX, LATEX, Metafont, Metapostÿ. Adresa je http://www.cstug.cz/documentation-index.html.
5.2
MetaPost
MetaPost je programovací jazyk určený pro popis a kreslení obrázků. Napsán byl Johnem Hobbym. Je odvozený z jazyku Metafont, který vytvořil spolu s TEXem Donald Knuth. MetaPost je též překladač tohoto jazyka. Je to perfektní nástroj pro kresbu diagramů v Postscriptu pomocí vektorového geometrického popisu. Metapost bývá součástí TEXových distribucí. Najdete jej ve složce bin jako mpost.exe nebo jen mp.exe. MetaPostový výstup může generovat též Gnuplot. Vhodnou literaturu naleznete opět na CSTUGu. Alespoň si ukažme, že se dá velice jednoduše nakreslit pro fyzika vážně zajímavý graf. Jsou to náhodně umístěné body nejprve dle normálního Gaussova rozdělení, a potom podle rovnoměrného rozdělení ve 2D. Kód programu v MetaPostu je následující: prologues:=1; u=1mm; beginfig(1); pickup pencircle scaled .3u; draw (-50u, 0)--(50u, 0); draw (0, -50u)--(0, 50u); pickup pecircle scaled 3u; ∆ jen o trochu větší než ∆0 , pak bez zavedení koeficientu k se krok zopakuje, ale jelikož závislost chyby na velikosti kroku není zcela přesně taková, jak předpokládáme, může i při příštím pokusu vyjít ∆ > ∆0 a za určitých okolností by se program mohl dokonce zacyklit. Zavedením koeficientu k v řádu jednotek náchylnost metody k zacyklení snížíme. Přesto je vhodné přidat i nějaký další ochranný prvek v podobě počítadla za sebou jdoucích opakování téhož kroku, kterých by nemělo být příliš mnoho.
38
for i=1 upto 100: draw(15u*normaldeviate, 15u*normaldeviate); endfor; endfig; beginfig(2); pickup pencircle scaled .3u; draw (-50u, 0)--(50u, 0); draw (0, -50u)--(0, 50u); pickup pecircle scaled 3u; for i=1 upto 100: draw(15u*uniformdeviate(5), 15u*uniformdeviate(5)); endfor; endfig; end;
39
Reference [1] Numerical Recipes in C: The Art of Scientific Computing, http://www.nr.com, Cambridge University Press 1988-92 [2] Slavíček, L.: Výpočetní technika pro chemiky, SNTL, Praha 1983
40
[3] Rybička, J.: LATEX pro začátečníky, Konvoj, Brno 2003
41
