Programování pro fyziky 2015/16: cvičení s L. H.

Tato stránka shrnuje průběh cvičení k předmětu Programování pro fyziky cvičícím L. H. Primárním zdrojem informací je samotná přednáška, dvojí poznámky k ní (stránka L. H. a stránka paralelky T. L.) a učebnice Numerical Recipes. Zde se zabýváme více technickými otázkami a vůbec zábavou kolem praktického programování numerických aplikací. Používáme Pascal (Free Pascal s vývojovým prostředím Lazarus) jako tradiční výukový jazyk (a volbu, zda C nebo C++ nebo Fortran necháváme na pozdější čas). Omezujeme se na programování typické pro fyzika, tedy konzolové aplikace doprovázené následnou vizualizací jinými programy – zde gnuplotem a ParaView. Příležitostně se obracíme k systému GNU Octave (nevlastní bratr MATLABu), velmi užitečným nástrojem pro numerické aplikace. Tato stránka je vytvořena v systému Sphinx na bázi ještě jiného jazyka oblíbeného i mezi fyziky – Pythonu, na který nám čas během semestru nevyjde už vůbec. Odkazy na výše vyznačené názvy hledejte vespod tohoto textu.

e-mail pro zasílání domácích úkolů: ladislav.hanyk@mff.cuni.cz

Cvičení č. 1

Na prvním cvičení jsme se usadili ve vývojovém prostředí Free Pascalu Lazarus, http://www.lazarus-ide.org, a sestavili tři zdrojové texty pro součet aritmetické posloupnosti podle fortranských vzorů z úvodní kapitoly poznámek k přednášce, http://geo.mff.cuni.cz/~lh/NOFY056. Ony tři programy v Pascalu následují, u poznámek k přednášce je ostatně i link na zip s celou sérií, v Pascalu i jiných jazycích. “Domácí úkol”: najděte si cestu k Lazaru, přeložte a spusťte si ty tři programy. Mně nic posílat nemusíte, ale můžete.

Free Pascal a Lazarus

Free Pascal je patrně jediná aktuálně vyvíjená implementace jazyka Pascal, která je volně dostupná, navíc pro Windows, Linux i Mac. (Object Pascal v Delphi je komerční.) Její užitečnost násobí dostupnost integrovaného vývojového prostředí (“IDE”) – Lazaru. Postačí získat a instalovat Lazarus, Free Pascal je součástí balíčku.

Po prvním spuštění Lazaru je dobře možné ignorovat vše, co Lazarus podsouvá, a otevřít nový projekt (volba Projekt–Nový projekt–Jednoduchý program). Nabídnutou kostru budoucího programu lze doplnit na tvar:

// První program
program Project;
begin
  writeln('Ahoj');
  readln;
end.

uložit (volba Projekt–Uložit projekt nebo alespoň Soubor–Uložit neboli Ctrl+S), přeložit (volba Spustit–Přeložit neboli Ctrl+F9) a spustit (volba Spustit–Spustit neboli F9). Program si otevře okno, které díky příkazu readln setrvá na obrazovce i po doběhnutí programu, až do stisku klávesy Enter. Pokud je okénko programu (ve Windows) příliš malé, po kliknutí v jeho levém horním rohu lze volit jeho velikost v Properties–Font a pro příště i v Defaults–Font.

_images/console.png

Překladač Free Pascalu lze volat také z příkazového řádku, bez Lazaru. (Překladače programovacích jazyků to běžně umožňují, některé ani žádné IDE nemají.) Pokud je operačnímu systému známa (tj. je uvedena v proměnné PATH; Windows: Properties of Computer–Advanced system settings–Environment variables) cesta k překladači – v mých Windows je to C:\lazarus\fpc\2.6.4\bin\i386-win32, pak stačí otevřít příkazové okno (Windows: tlačítko Start–All Programs–Accessories–Command Prompt nebo v okénku Run zadat cmd nebo použít vhodnou nabídku souborového manažeru, jako je Total Commander nebo Altap Salamander), dostat se do adresáře se zdrojovým textem – má příponu .lpr (Windows: pro přemístění se k souboru např. D:\Programovani\cviceni01\project.lpr zadat v příkazovém okně D: a change directory, cd \Programovani\cviceni01) a spustit z příkazového řádku Free Pascal Compiler:

fpc project.lpr

(Časem, pro složitější programy, bude třeba doplňovat překladový řádek o různé volby, které Lazarus v tichosti přidává sám, např. pro režim syntaxe Object Pascal: fpc -Mobjfpc project.lpr.) Překladač vytvoří proveditelný soubor project.exe a ten lze spustit opět z příkazového řádku (už bez přípony):

project

Pokud je ve zdrojovém textu závěrečný řádek readln, je k ukončení programu třeba stisknout Enter i zde.

Součet aritmetické posloupnosti

V první ukázce pro řešení \(S_n=\sum_{i=1}^n i\) jsme naznačili, že numerickému programování by vždy měla předcházet úvaha, zda nelze zadání formulovat matematicky ekvivalentně, ale výpočetně úsporněji. Zde \(S_n=\frac{n(n+1)}2\) řeší úlohu 3 operacemi, zatímco předchozí vzorec vyžaduje počet operací úměrný \(n\):

// sum_{i=1}^n i
// č. 1: matematická optimalizace
program p01a;
var
  n : integer;
  s : real;
begin
  write('Zadejte prosim cele cislo: ');
  readln(n);         // prikaz vstupu
  s:=n*(1.+n)/2;     // prirazovaci prikaz s ochranou proti integer přetečení
  writeln(n:7,s:15); // prikaz vystupu
  readln;
end.

V ukázce jsme kladli důraz na zdvořilost a informativnost konverzace rozvinuté programem (tedy nikoliv pouhé readln(n) bez předchozí výzvy write) a na prevenci celočíselného přetečení (integer overflow) pro větší n. Při n : integer máme pro n vyhrazeny 4 byty neboli 32 bitů a rozsah hodnot \(-2^{31}\)\(2^{31}-1\). Výraz n*(n+1)/2 pak přeteče už při \(n\ge46341\), měli bychom tedy preventivně přejít k vyčíslování v reálné přesnosti s rozsahem (při 8 bytech) až k \(10^{307}\). Toho lze dosáhnout např. zápisem vhodného (tedy brzy použitého) literálu o reálném datovém typu (jako výše: místo 1 jsme napsali 1.), nebo explicitní konverzí vhodného operandu na reálný typ, např. real(n)*(n+1)/2, nebo deklarací reálné proměnné rn : real, přiřazením rn:=n a následným používáním této proměnné: rn*(rn+1)/2. Doplňme ještě, že se Pascal celkem exoticky brání zápisu (n+1.), snáší přitom (n+1.0), (n+1. ) i zde použité (1.+n). V jiných jazycích se s přecitlivělostí na zápis (1.) hned tak nesetkáme.

V druhé ukázce jsme předvedli a odmítli špagetové programování, tedy styl, kdy se namísto běžných příkazových konstrukcí užívají příkazy skoku na návěští. Můžeme si ovšem celkem oprávněně představovat, že právě takto nějak pracuje procesor, že přeložený kód se příkazy skoku hemží. Skoky nevadí procesoru, vadí (dříve či později) programátorovi.

// sum_{i=1}^n i
// č. 2: špagetové programování s příkazy skoku
program p01b;
label
  L111,L999;
var
  n,i : integer;
  s : real;
begin
  write('Zadejte prosim cele cislo: ');
  readln(n);
  i:=0;
  s:=0;
  L111: if i>=n then goto L999;
  i:=i+1; s:=s+i;
  goto L111;
  L999: writeln(n:7,s:15);
  readln;
end.

V třetí ukázce jsme použili jednu z hlavních příkazových konstrukcí strukturovaného programování, úlohu řešili pomocí indexovaného cyklu. Ten zabezpečuje automaticky jak inkrementaci indexu, tak testování, zda index nedosáhl hodnoty horní meze. Mimochodem, v této ukázce přetečení nehrozí, neboť ani pro \(n=2^{31}-1\) kumulativní vyčíslování do reálné proměnné, s:=s+i, nepřeteče. Samozřejmě, při zadání např. n=3000000000 k přetečení dojde, a to hned na řádku readln(n). \((3\times10^9 > 2^{31})\)

// sum_{i=1}^n i
// č. 3: strukturované programování
program p01c;
var
  n,i : integer;
  s : real;
begin
  write('Zadejte prosim cele cislo: ');
  readln(n);
// indexovany cyklus
  s:=0;
  for i:=1 to n do begin
    s:=s+i;
  end;
  writeln(n:7,s:15);
  readln;
end.

Poznamenejme závěrem, že při n=2000000000 program vykoná několik miliard operací, což při frekvenci procesoru v GHz zabere nejvýše několik sekund.

Cvičení č. 2

Pokračujeme s programováním součtů řad, tentokrát s líbivým výsledkem. Náš program je zatím nestrukturovanou sekvencí myšlenek, různá drobná úskalí, která musíme překonávat, nás dostatečně saturují.

\(\pi\): Machin, Leibniz, Euler

Na druhém cvičení jsme navázali výpočtem \(\pi\) podle vzorců vybíraných ze souhrnu v PDF, konkrétně podle Machina,

\[\frac\pi4=\arctan\frac12+\arctan\frac13=2\arctan\frac12-\arctan\frac17=2\arctan\frac13+\arctan\frac17=4\arctan\frac15-\arctan\frac1{239}\]

Leibnize

\[\frac\pi4=\sum_{n=1}^\infty\frac{(-1)^{n-1}}{2n-1}=1-\frac13+\frac15-\frac17+\cdots\]

a Eulera:

\[\frac{\pi^2}6=\sum_{n=1}^\infty\frac1{n^2}=1+\frac1{2^2}+\frac1{3^2}+\frac1{4^2}+\cdots\]
// Výpočty pí podle http://geo.troja.mff.cuni.cz/~lh/NOFY056/cviceni-2015/_static/piFormulas.pdf

program p02;
uses
  Math;  // pripojeni modulu s funkci power

var
  nmax,n,z : integer;  // celociselne promenne
  s,rn : real;         // realne promenne
  piArctan,piMachin,piLeibniz,piEuler : real;

begin
// standardni funkce prekladace pro pi
  writeln('pi podle Lazaru:':30,pi:20:15);

// pi vypoctem z arctan
  piArctan:=arctan(1)*4;
  writeln('pi z atan:':30,piArctan:20:15);

// jeden z Machinovych vzorcu
  piMachin:=(arctan(0.5)+arctan(1/3))*4;
  writeln('pi podle Machina:':30,piMachin:20:15);

// aproximace souctu Leibnizovy rady
  nmax:=1000000;
  s:=0;
  z:=1;
  for n:=1 to nmax do begin
//  umocnovani je prilis pracna cesta pro pocitacove vycisleni alternujiciho znamenka
//  s:=s+power(-1,n-1)/(n*2-1);
    s:=s+z/(n*2-1);
    z:=-z; // alternace znamenka
    end;
  piLeibniz:=s*4;
  writeln('pi podle Leibnize:':30,piLeibniz:20:15);

// aproximace souctu Eulerovy rady
  nmax:=1000000;
  s:=0;
  for n:=1 to nmax do begin
    rn:=n;
    s:=s+1/sqr(rn);
    end;
  piEuler:=sqrt(s*6);
  writeln('pi podle Eulera:':30,piEuler:20:15);

  readln;
end.

Plyne odtud např. následující:

  • Konstanta (konstantní funkce) pi je v jazycích často předdefinovaná.
  • Lazarus některé matematické funkce umí hned (jsou součástí vždy připojené jednotky System), jiné až po připojení jednotky Math (uses Math). Do první série patří např. (absolutní hodnota) abs, (zaokrouhlování) round, (řezání desetinné části) trunc, (druhá mocnina) sqr, (druhá odmocnina) sqrt, (exponenciála) exp, (přirozený logaritmus) ln, sin, cos a zde použitý arctan, do druhé série např. (desítkový logaritmus) log10, tan, arcsin, arccos a (umocňování) power. Název arctan není příliš univerzální, v některých jazycích se říká atan.
  • Funkce arctan (a nejen ona) očekává reálný argument, vložení celočíselného literálu, arctan(1), ovšem v Pascalu nevadí – provede se automatická konverze datového typu. V některých jazycích se konverze argumentů funkcí a procedur automaticky neprovádí a je pak třeba si zvyknout na zápis arctan(1.).
  • Celočíselné dělení: v některých jazycích se operátor dělení / při integer operandech chápe jako celočíselné dělení (s celočíselným výsledkem). Je pak 1/3 rovna 0, arctan(1/3) samozřejmě také, a nezbývá než dbát o reálný typ alespoň jednoho operandu, tedy zde otečkovat literály: arctan(1./3). V Pascalu stojí operátor / vždy pro reálné dělení (operátor celočíselného dělení je div) a stačí proto psát arctan(1/3).
  • Řady s alternujícím znaménkem: v matematice píšeme \((-1)^n\) a je jasné, že jde jen o střídání znaménka. Programovací jazyk si to typicky vyloží jako požadavek na logaritmování, násobení a vyčíslení exponenciály, případně jako n-1 násobení. Každopádně počítání součtu Leibnizovy řady při použití funkce power je řádově pomalejší než při zajištění střídání znaménka způsobem užitým v ukázce nebo např. také takto:
    • if n mod 2=0 then z:=1 else z:=-1;
    • z:=(n+1) mod 2*2-1;
  • A stejně jako na prvním cvičení: pozor na celočíselné přetékání. Zatímco typ integer bezpečně stačí k ukládání čísel řádu milionů až (dvou) miliard, při výskytu mocnic samozřejmě narazí mnohem dříve: pro var n : integer = 46340; se výrazy sqr(n) a sqr(n+1) vyhodnotí jako 2147395600 a -2147479015 (druhý výraz přetekl, záporný výsledek je jistě špatně).

Za domácí úkol zkuste doplnit uvedený zdrojový text tak, aby počítal \(\pi\) i aproximací Wallisova nekonečného součinu:

\[\frac\pi2=\prod_{n=1}^\infty\frac{(2n)^2}{(2n-1)(2n+1)} =\frac{2\cdot2}{1\cdot3}\cdot \frac{4\cdot4}{3\cdot5}\cdot\frac{6\cdot6}{5\cdot7}\cdot \frac{8\cdot8}{7\cdot9}\cdots\]

Více návodů, jak počítat \(\pi\), naleznete leckde: MathWorld, Wikipedia, Piphilology.

Cvičení č. 3

Počítáme další řady pro \(\pi\). Výsledný program už získá jistou strukturu, každá myšlenka bude zabalena do vlastní programové jednotky. Nazrál čas i na počty, které lze vizualizovat: Lissajousovy obrazce. Pro kreslení budeme teď i potom používat gnuplot a občas, jako bonbonek, ParaView.

\(\pi\): Viète, Ramanujan

Na dalším cvičení jsme pokračovali ve výpočtech \(\pi\) podle vzorců ze souhrnu v PDF, a to Viètovými odmocninami

\[\frac2\pi=\sqrt{\frac12}\sqrt{\frac12+\frac12\sqrt{\frac12}} \sqrt{\frac12+\frac12\sqrt{\frac12+\frac12\sqrt{\frac12}}}\cdots\]

a fascinující Ramanujanovou řadou.

\[\frac1\pi=\frac{\sqrt8}{9801}\sum_{n=0}^\infty\frac{(4n)!}{(n!)^4}\frac{1103+26390n}{396^{4n}} =\frac{\sqrt8}{9801}\left(1103+\frac{54986}{2049271488}+\cdots\right)\]

V obou případech už se nehodí programovat n-tý člen řady explicitně, lépe je vyjádřit jej pomocí (n-1)-ho členu:

\[\begin{split}\mbox{Viete}\quad & a_n=\sqrt{\frac12+\frac12 a_{n-1}},\quad a_1=\frac12 \\ \mbox{Ramanujan}\quad & a_n=a_{n-1}\frac{(4n-3)(4n-2)(4n-1)(4n)}{(396n)^4},\quad a_0=1\end{split}\]

V prvním případě ani jiné řešení nevidíme, v druhém vzorci nás k tomuto postupu vedla hrozba separátního přetékání v čitateli i jmenovateli při počítání podle zadání; nekonfliktní faktor (1103+26390n) jsme vyčlenili stranou. I tak se ovšem jistíme proti celočíselnému přetékání používáním reálných literálů, např. (4.*n-3).

Jednotlivé výpočty jsme nyní strukturovali jako (programové jednotky) funkce, parametrizované svými argumenty. Popisy použitých proměnných jsme vsunuli do těla funkcí, omezili jsme tak jejich platnost na blok funkce v zájmu prevence konfliktů jmen. Výsledný program:

// Výpočty pí podle http://geo.troja.mff.cuni.cz/~lh/NOFY056/cviceni-2015/_static/piFormulas.pdf
// Dílčí výpočty strukturovány do funkcí
program p03;
uses
  Math;  // pripojeni modulu s funkci power

// pi vypoctem z arctan
function piArctan : real;
begin
  piArctan:=arctan(1)*4;
end;

// dva z Machinovych vzorcu
function piMachin(const branch : integer) : real;
begin
  case branch of
  1: piMachin:=(arctan(0.5)+arctan(1/3))*4;
  2: piMachin:=(2*arctan(0.5)-arctan(1/7))*4;
  else
     piMachin:=0;
  end;
end;

// aproximace souctu Leibnizovy rady
function piLeibniz(const nmax : integer) : real;
// uses Math; // nutno pripojit zkraje celeho programu
var
  s,z : real;
  n : integer;
begin
  s:=0;
  z:=1;
  for n:=1 to nmax do begin
    s:=s+z/(n*2-1);
    z:=-z; // alternace znamenka
    end;
  piLeibniz:=s*4;
end;

// aproximace souctu Eulerovy rady
function piEuler(const nmax : integer) : real;
var
  s,rn : real;
  n : integer;
begin
  s:=0;
  for n:=1 to nmax do begin
    rn:=n;
    s:=s+1/sqr(rn);
    end;
  piEuler:=sqrt(s*6);
end;

// aproximace Vietova soucinu
function piViete(const nmax : integer) : real;
var
  s,a : real;
  n : integer;
begin
  a:=sqrt(0.5);
  s:=a;
  for n:=2 to nmax do begin
    a:=sqrt(0.5+0.5*a);
    s:=s*a;
    end;
  piViete:=2/s;
end;

// aproximace Ramanujanovy rady
function piRamanujan(const nmax : integer) : real;
var
  s,a : real;
  n : integer;
begin
  a:=1;
  s:=a*1103;
  for n:=1 to nmax do begin
    a:=a*(4.*n-3)*(4.*n-2)*(4.*n-1)*4*n/sqr(sqr(396.*n));
    s:=s+a*(1103+26390*n);
    end;
  result:=1/(s*sqrt(8)/9801); // result synonymem pro jmeno navratove hodnoty funkce
end;

begin
  writeln('pi podle Lazaru:':30,pi:20:15);
  writeln('pi z atan:':30,piArctan:20:15);
  writeln('pi podle Machina:':30,piMachin(2):20:15);
  writeln('pi podle Leibnize:':30,piLeibniz(1000000):20:15);
  writeln('pi podle Eulera:':30,piEuler(1000000):20:15);
  writeln('pi podle Vieta:':30,piViete(100):20:15);
  writeln('pi podle Ramanujana:':30,piRamanujan(10):20:15);
  readln;
end.

Ve funkci piRamanujan jsme si zkusili značit návratovou proměnnou nikoliv jménem funkce, ale univerzálním jménem result. Jazyky to podporují, Lazarus v default nastavení (tj. v režimu syntaxe Object Pascal) také.

Lissajousovy obrazce

Počítání jako výše, tedy výpočet jediného čísla – navíc předem přibližně známého – je ve fyzice celkem vzácnost. Typičtější je chrlení množství dat, které je pak třeba dále analyzovat, např. prostou vizualizací. To nyní uděláme.

Budeme generovat časový vývoj dvou souřadnic popsaný goniometrickými funkcemi a kreslit dráhu, po které se bude pohybovat bod \((x,y)\):

\[\begin{split}x(t) & = a_x \sin(f_xt+\varphi_x) \\ y(t) & = a_y \sin(f_yt+\varphi_y)\end{split}\]

Taková dráha, zvláště je-li podíl frekvencí v x a y směru racionálním číslem, se nazývá Lissajousův obrazec (Wikipedia). Program

// Příprava dat pro Lissajousovy obrazce
// https://en.wikipedia.org/wiki/Lissajous_curve
program lissajous;
const
  nmax = 1000; // pocet kroku
  ax = 1;      // x amplituda
  ay = 1;      // y amplituda
  fx = 3;      // x frekvence
  fy = 2;      // y frekvence
  phx = 0;     // x faze
  phy = 0;     // y faze
  vystup = 1;  // 1: na obrazovku, 2: do souboru
  soubor = 'lissajous.dat';
var
  n : integer;
  t,dt,x,y : real;
  s : text;
begin
  case vystup of
  2: begin assign(s,soubor); rewrite(s) end;
  end;
  dt:=pi*2/nmax;
  t:=0;
  for n:=0 to nmax do begin
    t:=t+dt;
    x:=ax*sin(fx*t+phx);
    y:=ay*sin(fy*t+phy);
    case vystup of
    1: writeln(x:8:4,y:8:4);
    2: writeln(s,x:8:4,y:8:4);
    end;
  end;
  case vystup of
  2: begin close(s); writeln('Vytvoren soubor ',soubor,'.') end;
  end;
  readln;
end.

si nese vstupní data ve svém těle (ideálně někde dobře na očích, zde hned nahoře) a vypisuje výstupní data volitelně buď na obrazovku nebo do souboru (lissajous.dat). Jakmile máte data v souboru, můžete je vykreslit např. v gnuplotu.

Získat vizualizační program gnuplot (http://www.gnuplot.info) není trnité. Po instalaci (na mysli máme dále Windows, v Linuxu to není moc jiné) se mezi programy nabídky Start objeví gnuplot s položkou gnuplot 4.6 (nebo jiná verze), která spustí wgnuplot.exe. Kde se gnuplot spustil, na to odpoví příkaz pwd (print working directory). Pokud jsou Vaše data (soubor lissajous.dat) v adresáři např. D:\User\Data, je třeba zadat příkaz cd 'D:\User\Data' (change directory). Je samozřejmě také možné nastavit si pracovní adresář v Properties položky gnuplot 4.6 v nabídce Start, nastavením adresáře Start in, nebo přidat adresář s wgnuplot.exe do proměnné PATH způsobem, který jsme popsali na prvním cvičení. Pak už by nic nemělo stát v cestě příkazu plot 'lissajous.dat' with lines neboli p 'lissajous.dat' w l, kterým se otevře vizualizační okno s výsledkem. Z toho se do příkazového okna lze vrátit mj. stisknutím mezerníku, a můžete pokračovat dalším příkazem nebo skončit zadáním quit neboli q.

_images/gnuplot.png

Sekvenci příkazů vedoucích až k obrázku v PNG formátu shrnuje následující skript:

# gnuplot: kreslení dat ze souboru na obrazovku a v PNG formátu
plot 'lissajous.dat' with lines
set terminal png
set output 'lissajous.png'
replot
set output

# zkracovaný zápis
# p 'lissajous.dat' w l
# set term png
# set o 'lissajous.png'
# rep
# set o

který buď spustíte v gnuplotu příkazem load 'lissajous.gn' nebo se obsahem skriptu inspirujte a v něm uvedené příkazy zadejte v gnuplotu ručně.

_images/lissajous-gnuplot.png

Lissajousův obrazec pro poměr frekvencí v x a y směru 3:2. Kresleno v gnuplotu.

A jistě, kdybychom si zrovna neprocvičovali programování v klasickém programovacím jazyce, mohli bychom Lissajousovy obrazce kreslit v gnuplotu rovnou. Gnuplot rozumí funkcím zadaným analyticky a umí i s proměnnými, takže např. omega=1; plot sin(omega*x) vykreslí goniometrickou funkci v default intervalu <-10,10>. (x je default symbol gnuplotu pro nezávisle proměnnou.) Gnuplot rozumí i parametrickému zápisu funkcí, takže výše předložený Lissajousův obrazec lze získat i na zelené louce tímto skriptem (při parametrickém zápisu je default symbolem pro parametr t):

# gnuplot: kreslení funkce zadané analyticky, zde navíc parametricky
set parametric
set samples 1000
fx=3
fy=2
plot sin(fx*t),sin(fy*t)

Za domácí úkol upravte zdrojový text ze cvičení tak, aby počítal trojrozměrné Lissajousovy obrazce (nebo rovnou Lissajousovy uzly). Pro jejich kreslení můžete opět využít gnuplot: splot 'lissajous.dat'. V interaktivním režimu gnuplotu si 3D obrázkem můžete pomocí myši otáčet. Alternativně se můžete podívat na výsledek ve směru jednotlivých os pomocí průmětů, tedy výběrem různých párů z vypočtených tří sloupců: plot 'lissajous.dat' using 1:2 with lines; p 'lissajous.dat' u 1:3 w l; p 'lissajous.dat' u 2:3 w l.

Můžete si také prohlédnout dvoustránkový výtah z dokumentace gnuplotu.

ParaView

Gnuplot je silný v dvourozměrných (2D) obrázcích, slabší (co do kvality výstupu) ve třech dimenzích (3D). Pro 3D vizualizaci je samozřejmě na výběr více cest. Jednou z nejhodnotnějších (a přitom volně dostupných, a pro Windows/Linux/Mac) je ParaView. Tak proč si 3D Lissajousův obrazec neprohlédnout v ParaView?

Vstupní data má ParaView nejraději ve VTK formátu. Pro nakreslení série 3D bodů nám půjde o variantu zvanou polygonální síť alias polydata. Vlastně bude stačit vložit na počátek datového souboru pro gnuplot, který máte získat jako domácí úkol, hlavičku:

# vtk DataFile Version 3.0
comment (max. 1023 chars)
ASCII
DATASET POLYDATA
POINTS 1001 float

a souboru dát příponu vtk. (Budete-li generovat jiný počet bodů než 1001, upravte údaj v hlavičce.) V ParaView pak můžete postupovat třeba takto:

  • File–Open–’lissajous.vtk’, vlevo uprostřed tlačítko Apply
  • Filters–Alphabetical–Glyph (nebo nahoře ve třetím řádku ikonek sedmá ikonka zleva), Apply
  • a dál v okně Properties... Glyph Type: Sphere, Apply
  • (pro zvětšení koulí) Scaling–Scale Factor, Apply
  • (pro pravidelný rozestup koulí) Masking–Glyph Mode–All Points, Apply
  • (pro obarvení koulí) Coloring–Edit a Background–Color

a získáte třeba toto:

_images/lissajous-paraview.png

3D Lissajousův obrazec pro poměr frekvencí v x, y a z směru 3:5:7. Kresleno v ParaView.

S obrázkem si v okně ParaView můžete točit, posouvat atd., můžete ho uložit jako bitmapu (File–Save Screenshot) a můžete si také třeba nastavit View–Animation View, pak u modrého křížku dole Camera–Orbit–klik na modrý křížek–OK, ještě Mode–Real Time–Duration: 10, nahoře zelenou šipkou Play a dál už si založit ruce a sledovat. Nebo ještě File–Save Animation a získat toto.

K ParaView a kreslení ve třech rozměrech se ještě vrátíme.

Cvičení č. 4

Cílem je vytvořit komplexní aritmetiku, zapouzdřit ji do modulu a zhotovit obrázek Mandelbrotovy množiny, podobný obrázku z Wikipedie.

Komplexní aritmetika

Čistý Pascal se (na rozdíl od C, Fortranu, Octave ad.) o komplexní datový typ a tím ani o počítání s komplexními čísly nesnaží. Snažit se budeme my. Odhlédneme od toho, že k Free Pascalu je přibalen (dodnes nedokumentovaný) modul UComplex (soubor ucomplex.pp), který se s tématem utkává, a připravíme si potřebný typ, operátory a funkce sami. Komplexní čísla jsou dvojicemi reálných čísel, můžeme je tak reprezentovat pomocí polí i záznamů:

type
  complexArr = array [1..2] of real;
  complexRec = record re,im : real end;

Zvolíme druhou cestu. Pro převody mezi reálnými a komplexními čísly se inspirujeme ve Fortranu a konverzní funkce nazveme cmplx, areal a imag. Přidáme funkce pro sčítání, odčítání a násobení, některé přetížíme (tj. násobně definujeme) tak, aby pracovaly s komplexními i reálnými argumenty, a provážeme je s běžnými aritmetickými operátory. Doplníme funkci pro výpočet absolutní hodnoty. To vše sbalíme do vlastního samostatného modulu (programová jednotka unit) tak, abychom se k tomu mohli snadno vracet.

// Komplexní datový typ, operátory a funkce
unit uMyComplex;

interface

type complex = record re,im : real end;
function cmplx(const a,b : real) : complex;
function areal(const x : complex) : real;
function imag(const x : complex) : real;
function cAdd(const x,y : complex) : complex; overload;
operator + (const x,y : complex) : complex; overload;
function cAdd(const x : complex; const y : real) : complex; overload;
operator + (const x : complex; const y : real) : complex; overload;
function cSub(const x,y : complex) : complex;
operator - (const x,y : complex) : complex;
function cMul(const x,y : complex) : complex; overload;
operator * (const x,y : complex) : complex;
function cMul(const x : complex; const y : real) : complex; overload;
operator * (const x : complex; const y : real) : complex; overload;
function cAbs(const x : complex) : real;

implementation

function cmplx(const a,b : real) : complex;
begin
  result.re:=a; result.im:=b;
end;

function areal(const x : complex) : real;
begin
  result:=x.re;
end;

function imag(const x : complex) : real;
begin
  result:=x.im;
end;

function cAdd(const x,y : complex) : complex; overload;
begin
  with result do begin re:=x.re+y.re; im:=x.im+y.im end;
end;

operator + (const x,y : complex) : complex; overload;
begin
  result:=cAdd(x,y);
end;

function cAdd(const x : complex; const y : real) : complex; overload;
begin
  with result do begin re:=x.re+y; im:=x.im end;
end;

operator + (const x : complex; const y : real) : complex; overload;
begin
  result:=cAdd(x,y);
end;

function cSub(const x,y : complex) : complex;
begin
  with result do begin re:=x.re-y.re; im:=x.im-y.im end;
end;

operator - (const x,y : complex) : complex;
begin
  result:=cSub(x,y);
end;

function cMul(const x,y : complex) : complex; overload;
begin
  with result do begin re:=x.re*y.re-x.im*y.im; im:=x.re*y.im+x.im*y.re end;
end;

operator * (const x,y : complex) : complex; overload;
begin
  result:=cMul(x,y);
end;

function cMul(const x : complex; const y : real) : complex; overload;
begin
  with result do begin re:=x.re*y; im:=x.im*y end;
end;

operator * (const x : complex; const y : real) : complex; overload;
begin
  result:=cMul(x,y);
end;

function cAbs(const x : complex) : real;
begin
  with x do begin
    if abs(re)>abs(im) then
      result:=abs(re)*sqrt(1+sqr(im/re))
    else if abs(im)>0 then
      result:=abs(im)*sqrt(1+sqr(re/im))
    else
      result:=0;
  end;
end;

end.

Vnímáme, že moduly mají část s veřejným rozhraním (interface) a část se soukromou implementací (implementation). Vně modulu lze použít jen to, co je uvedeno v jeho rozhraní. Vnímáme také, že u přetížených funkcí a operátorů uvádíme klauzuli overload, ve Free Pascalu nepovinnou (v Delphi povinnou). A ještě vnímáme, že se vyhýbáme výpočtu absolutní hodnoty přímočarým vzorcem result:=sqrt(sqr(x.re)+sqr(x.im)), kde by umocnění jednotlivých složek mohlo někdy způsobit přetečení. (Při var x : real se po x:=sqr(1e300); x:=sqrt(x);x vždy objeví Infinity, neboť jsme vynutili uložení 1e600 do 8bytové reálné proměnné, což přeteče. Po x:=sqrt(sqr(1e300)) bude v x zpátky 1e300, pokud jsou mezivýsledky ve výrazech vyčíslovány v 10bytové přesnosti extended – ve Free Pascalu jsou. S tím se potkáváme i v jiných jazycích: některé překladače nechávají mezivýsledky výrazů vyčíslovat v přesnosti operandů, jiné v maximální dostupné přesnosti.)

Modul můžeme částečně prověřit následujícím testovacím programem:

// Minitest modulu uMyComplex
program testComplex;
uses uMyComplex;
var c : complex;
begin
  c:=cmplx(0,1);
  writeln(areal(c),imag(c));                  // 0 1
  writeln(areal(c*c),imag(c*c));              // -1 0
  writeln(cabs(cmplx(1,1)),cabs(cmplx(0,0))); // 1.4142 0
  readln;
end.

Mandelbrotova množina

Populární fraktál zvaný Mandelbrotova množina je množinou komplexních čísel \(c\), pro které posloupnost komplexních čísel \(z_n\),

\[z_{n+1}=z_n^2+c,\quad z_0=0,\]

pro \(n\to\infty\) zůstane omezená. Neomezenost posloupnosti pro pevné \(c\) lze po úvaze ztotožnit s existencí konečného \(N\), pro které (poprvé) nastane \(|z_N|>2\). Bod komplexní roviny do Mandelbrotovy množiny buď patří nebo nepatří, obrázky jako tento navíc probarvují vyloučené body barvou přiřazenou indexu \(N\) (a ještě přidávají spojité barevné přechody mezi diskrétními hodnotami \(N\)).

My nyní vytěžíme náš modul uMyComplex a gnuplot k získání podobného obrázku. V následujícím programu se zvolený obdélník v komplexní rovině pokrývá náhodně vybíranými body, pro které se provede až nmax kroků uvedené posloupnosti a testuje se splnění uvedené podmínky. Na výstup se posílají souřadnice testovaného bodu a min (\(N\), nmax+1). Funkce random vrací (při každém volání jiné) pseudonáhodné reálné číslo z intervalu \(\langle 0,1)\), procedura randomize inicializuje generátor náhodných čísel (tj. nastaví náhodný výchozí člen jejich posloupnosti).

// Příprava bodů pro kreslení Mandelbrotovy množiny
program Mandelbrot;
uses uMyComplex;
const
  xmin=-2;
  xmax=0.6;
  ymin=-1.2;
  ymax=1.2;
  nmax=100;
var
  x,y : real;
  c,z : complex;
  n,nstop : integer;
begin
  randomize;
  while true do begin
    x:=xmin+(xmax-xmin)*random; // přepočet na <xmin,xmax)
    y:=ymin+(ymax-ymin)*random; // přepočet na <ymin,ymax)
    c:=cmplx(x,y);
    z:=cmplx(0,0);
    nstop:=0;
    for n:=1 to nmax do begin
      z:=z*z+c;
      if cAbs(z)>2 then begin nstop:=n; break end;
    end;
    if nstop=0 then nstop:=nmax+1;
    writeln(x:8:5,y:9:5,nstop:5);
  end;
end.

Ve vnitřním cyklu se mnohokrát provádějí 2 řádky zahrnující zhruba 10 elementárních operací (real sčítání a násobení) a 1 real odmocninu. Obojí může být svou časovou náročností srovnatelné a bylo by proto efektivní nahradit zbytnou odmocninu v cAbs(z)>2 podmínkou cAbsSqr(z)>4, kde cAbsSqr by na rozdíl od cAbs neodmocňovalo.

Domácí úkol č. 4: uložte si výstupní data do souboru a pomocí gnuplotu a přiloženého skriptu vytvořte PNG obrázek. Soubor s daty můžete zapisovat po doplnění programu tak, jak jsme provedli už v minulém úkolu, nebo můžete program (mandelbrot.exe) spustit v příkazovém okně a výpis směřující na obrazovku přesměrovat do souboru:

mandelbrot > mandelbrot.dat

Nekonečný běh ukázky můžete přerušit klávesovou kombinaci Ctrl+C nebo počet testovaných bodů omezte v programu. Obrázek bude znehodnocen bílými tečkami, neboť náhodně prováděným výběrem bodů se v rozumném čase celá oblast nepokryje. Upravte tedy (součást domácího úkolu) program tak, aby se zvolený obdélník pokrýval systematicky a s předepsaným rozlišením. A také zoomujte obdélníkem blíže k hranici Mandelbrotovy množiny, ať vidíme více. Při hlubokém zoomování pozor na nmax a také formátovací specifikace ve writeln; obojí může být potřeba zvětšit.

Skript pro gnuplot definuje spektrum a rozsah barevné palety, nastaví výstup do PNG souboru o zvoleném rozlišení a volá příkaz plot pro datový soubor, v jehož prvních dvou sloupcích se očekávají souřadnice kresleného bodu a ve třetím sloupci index barvy:

# gnuplot: vykreslení Mandelbrotovy množiny
set palette defined \
  (0 'dark-blue',6 'dark-blue',30 'white',70 'gold',100 'brown',100.001 'black',101 'black')
set cbrange [0:101]
set term png size 1024,768
set output 'mandelbrot.png'
unset key
plot [-2.:.6][-1.2:1.2] 'mandelbrot.dat' using 1:2:3 palette z with dots
set output

Skript se v gnuplotu spustí již známým způsobem, load 'mandelbrot.gn', a při systematickém průchodu obdélníkem by výsledný obrázek mohl vypadat třeba takto:

_images/mandelbrot.png

Posloupnost \(z_n\) vedoucí k Mandelbrotově množině lze zobecnit a kreslit např. Multibrot sets. Pro dobrovolníky, které při práci na minulém domácím úkolu uhranulo ParaView, je k mání 3D zobecnění zvané Mandelbulb.

Cvičení č. 5

Obrátíme se k celým číslům a připravíme program pro nalezení všech prvočísel menších než nebo rovných N. Kreslit budeme (jako domácí úkol) jen dvě křivky: závislost počtu prvočísel a počtu prvočíselných dvojčat na N. Dáme si však v gnuplotu záležet a řádně popíšeme obrázek i osy. Zocelíme se také v práci s velkými datovými soubory.

Eratosthenovo síto

Připravíme si pole o N prvcích. Prvky budou (na závěr) obsahovat informaci Je můj index prvočíslem? ano/ne. Na začátku do všech prvků vložíme s důvěrou ano a vyškrtneme neprvočíslo 1 (ne do prvku s indexem 1). Pak v cyklu vyhledáme nejbližší větší nevyškrtnuté číslo n (prvek s indexem n), kterému ponecháme ano, ale vyškrtáme všechny jeho násobky (2n, 3n, ...), a opakujeme, než se dostaneme k N. Násobky jsou vyškrtány, zůstala – prvočísla. Vizuální ukázky a varianty algoritmu jsou na Wikipedii. Hodnotná je tato optimalizace: z pozice prvočísla n není třeba vyškrtávat jeho násobky menší než n-násobky, neboť již byly vyškrtnuty dříve; začínáme tedy vyškrtávat od kvadrátu n. To zároveň znamená, že pole procházíme ne až k N, ale jen k jeho odmocnině, neboť nad ní už je vyškrtáno.

Free Pascal vzdoruje snaze vytvořit pole o více než 2 miliardách prvků. Datový typ prvků použijeme co nejmenší, 1bytový. Přímočaré je volit pro uložení ano/ne typ boolean. Pro menší N můžeme prvočísla vypisovat, pro větší zjistíme jen jejich počet. (Nakonec si ale prvočísla pro vykreslení závislosti jejich počtu na N a pro trénink zacházení s velkými daty stejně uložíme do souboru všechny.)

// Eratosthenovo síto pro nalezení všech prvočísel menších než nmax
program Eratosthenes;
const
  nmax = 1000000;                 // beztypova konstanta nmax
  vypis_prvocisel : boolean = false;  // psat ci nepsat?
var
  a : array [1..nmax] of boolean;
  n,nn,count : integer;
begin
  writeln('Eratosthenovo sito');
  writeln('... prosivaji se cela cisla v intervalu <1,',nmax,'>');
  for n:=1 to nmax do a[n]:=true; // inicializujeme
  a[1]:=false;                    // 1 neni prvocislo
  n:=2;
  while n<=sqrt(nmax) do begin
    nn:=n*n;                      // od n-nasobku n
    while nn<=nmax do begin
        a[nn]:=false;             // vyskrtavame nasobky
        nn:=nn+n;
      end;
    n:=n+1;
    while not a[n] do begin       // nejblizsi vyssi prvocislo
      n:=n+1;
    end;
  end;
  count:=0;                       // pocet prvocisel
  for n:=1 to nmax do if a[n] then count:=count+1;
  writeln('... pocet prvocisel: ',count);
  writeln('... pocet prvociselnych dvojcat: ','?');
  if vypis_prvocisel then begin
    writeln('... seznam prvocisel:');
    for n:=1 to nmax do if a[n] then write(n:8);
    writeln;
  end;
  readln;
end.

Program na mém desktopu pro nmax rovno miliardě zjišťoval počet prvočísel při default nastavení překladače pod 30 s, při zvýšení na druhou úroveň optimalizace zrychlil pod 15 s. To je výrazné poselství: při delších počtech mysleme na možnost nechat překladač optimalizovat. (Lazarus: Projekt–Volby projektu–Kompilace a linkování–Úrovně optimalizace, příkazový řádek: fpc -Mobjfpc -O2 project.lpr.) Poznamenejme, že při podmínění vyškrtávacího řádku, if a[nn] then a[nn]:=false, tedy při přeskočení zbytečných zápisů do paměti, se program bez optimalizace zrychlil z 28 s na 20 s, zatímco po optimalizaci tato úprava řádku už žádné zrychlení nepřinesla. Vidíme tak zřejmě jeden z optimalizačních postupů Free Pascalu. O více technikách se můžeme dočíst na Wikipedii.

A výpis dat? Ukázka je nachystána pro úsporný vícesloupcový výpis na obrazovku, upravte si ji pro výpis do souboru. Vhodné pro gnuplot bude umístit každé prvočíslo na vlastní řádek (číslo řádku pak odpovídá počtu prvočísel menších než nebo rovných prvočíslu na řádku), čímž zabereme až 11 bytů na prvočíslo (9 bytů pro 9ciferná prvočísla a ve Windows 2 byty pro symbol konce řádku). Při miliardě v nmax a zhruba 5 procentech prvočísel v ní bude soubor velký přes půl gigabytu. Pro úspornější ukládání mohou být užitečné binární soubory, v kterých jsou data uložena ve stejné podobě jako v paměti počítače. To pak stojí jen 4 byty na integer (8 bytů na real), tedy v našem případě asi 200 megabytů. Binární soubor si sice sami snadno nepřečteme, ale stačí, umí-li jej přečíst gnuplot, a to on umí. Navíc, programy zapisují a čtou binární soubory zpravidla rychleji než textové, což kupodivu neplatí pro Free Pascal, ale pro gnuplot ano. Zápis do textových (neboli ASCII) a binárních souborů v Pascalu:

// Zápis do textového souboru
var s_ascii : text;
    soubor : string = 'soubor.dat';
begin
assign(s_ascii,soubor);
rewrite(s_ascii);
writeln(s_ascii,123456789:9);
close(s_ascii);
end.

// Zápis do binárního souboru
var s_binary : file of integer;
    soubor : string = 'soubor.bin';
begin
assign(s_binary,soubor);
rewrite(s_binary);
write(s_binary,123456789);
close(s_binary);
end.

Domácí úkol: doplňte program s algoritmem Eratosthenova síta tak, aby kromě počtu prvočísel uměl zjistit i počet párů prvočíselných dvojčat (prvočísel, jejichž rozdíl je roven 2) menších než N. Obě funkce, p(N) a d(N), vykreslete pro kladná N do miliardy nebo (v případě problémů) alespoň do 100 milionů.

Zbývá nám pravidelná týdenní dávka gnuplotu. Chceme vlastně málo, jen kreslit dvě obyčejné čáry, pro které očekáváme data ve dvou jednosloupcových souborech, buď textových, primes.dat a twins.dat, nebo binárních, primes.bin a twins.bin. Tak alespoň ještě doplníme dvouřádkový nadpis s českou diakritikou, popisy os včetně horních indexů, přesuneme legendu do jiného rohu a nastavíme čtvercový formát obrázku:

# gnuplot: počet prvočísel (primes) a prvočíselných dvojčat (twin prime pairs) menších než N
xmax=1e9      # horni limit x-osy
ymax=xmax/20  # horni limit y-osy
xsize=800
ysize=800
binary=0      # typ vstupnich souboru: 0 ascii, 1 binary
set title " p(N): počet prvočísel menších než N \n d(N): počet prvočíselných dvojčat menších než N "
set xlabel 'N'
set xtics xmax/5 format "%.0t.10^%1T"
set mxtics 4  # minor xtics
set ytics ymax/5 format "%.0t.10^%1T"
set mytics 4  # minor ytics
set key left top
set size square
set terminal png enhanced size xsize,ysize
set output 'eratosthenes.png'
if (!binary) \
plot [0:xmax][0:ymax] 'primes.dat' using 1:0 with lines title 'p(N)', \
                      'twins.dat' using 1:0 with lines title 'd(N)'; \
else \
plot [0:xmax][0:ymax] 'primes.bin' binary format='%i' using 1:0 with lines title 'p(N)', \
                      'twins.bin' binary format='%i' using 1:0 with lines title 'd(N)'
set output
  • Pomáháme si proměnnými inicializovanými zkraje skriptu, abychom čelili výskytu magických konstant hlouběji ve skriptu.
  • Při popisu os příkazy set xtics, set ytics používáme formátové specifikace vyložené v nápovědě gnuplotu: help format specifiers. Pro horní indexy uvozené znakem ^ a podobná vylepšení je třeba podpora výstupního terminálu, zde: set terminal png enhanced.
  • Legendu lze přesouvat: set key left top, různě vylepšovat: help set key, i vypnout: unset key.
  • Větvíme podmíněným příkazem if podle hodnoty v proměnné binary (operátor ! neguje 0 na 1 a zpět) na kreslení z textových nebo binárních souborů. Použili jsme jednořádkový if (pokračovací řádky po ukončení předchozího řádku znakem \ se nepočítají), moderní gnuplot umí i blokovou variantu se složenými závorkami: help if.
  • Příkaz plot obsahuje explicitní nastavení dolní a horní meze v obou osách: [xmin:xmax][ymin:ymax], odkazuje na virtuální nultý sloupec čili číslo řádku: using 1:0, a mění default legendu: title string.
  • Příkaz plot s klauzulí binary čte data z binárních souborů. Jejich strukturu je vhodné popsat pomocí formátových specifikací podle návodu po help plot binary general a kliknutí na format. Stručně: píšeme format='%i' pro čtení pascalského integer, %f pro single a %lf pro double neboli real.
  • Skript s českými texty musíme pro gnuplot uložit v kódování UTF-8 without BOM (byte order mark neboli signature). To např. v Notepadu ve Windows 7 nelze. Alespoň o důvod víc zkusit některý programátorský editor, např. Notepad2 (Wikipedia) nebo Notepad++ (Wikipedia).
  • 32bitový gnuplot pro N rovno miliardě (tedy už pro 200MB-ový sloupec) selže. Použijte 64bitový gnuplot nebo N zmenšete.
_images/eratosthenes.png

Funkce p(N) má své jméno: prvočíselná funkce, své písmeno: \(\pi(n)\), i svou webovou stránku: Wikipedia. Tam je k porovnání i tabulka jejích hodnot.

Cvičení č. 6

Jsme zpátky u počítání s reálnými a komplexními čísly. Připravíme programy pro numerický výpočet určitého integrálu reálné funkce jedné proměnné a numerické hledání nulového bodu reálné funkce reálné proměnné a komplexní funkce komplexní proměnné. Hledač kořenů (kostrbatě z anglického root finder) vytěžíme v komplexní rovině pro obrázek dalšího fraktálu.

Lichoběžníkové pravidlo

Chceme numericky vyčíslit \(\int_a^b f(x)\,dx\). Interval \(\langle a,b\rangle\) rozdělíme na \(N\) podintervalů téže délky, \(h=(b-a)/N\), jejich hraniční body označíme \(x_0=a,\ x_1,\ \dots,\ x_N=b\) a hodnoty integrandu v nich \(f_n=f(x_n),\ n=0,\dots,N\). Body \((x_n,f_n)\) vedeme po částech lineární polynom (lomenou čáru). Hledaný integrál pak lze aproximovat plochou pod touto lomenou čarou, v případě podintervalu \(\langle x_n,x_{n+1}\rangle\) plochou lichoběžníka \(\frac h2(f_n+f_{n+1})\), na celém intervalu pak souhrnem těchto dílčích vzorců (tzv. složené lichoběžníkové pravidlo neboli trapezoidal rule):

\[\int_a^b f(x)\,dx \approx h\left(\frac{f_0}2+f_1+f_2+\dots+f_{N-1}+\frac{f_N}2\right).\]

Počet lichoběžníků \(N\) z nebe nespadne, vhodnější je proto poručit si pro výsledek požadovanou přesnost \(\varepsilon\). Pak můžeme začít s \(N=1\) a následně \(N\) zdvojnásobovat tak dlouho, než nastane \(|L_N-L_{N/2}|<\varepsilon\). Jistě je to jen aproximace požadované přesnosti, ale běžná. Současně limitujeme \(N\) vhodnou horní mezí. Běžné je také, že knihovní procedury pro numerickou integraci žádají \(f(x)\) jako uživatelem poskytnutou funkci (procedurální argument), meze intervalu a požadovanou přesnost. Jméno funkce quad (jako kvadratura) převezmeme od hodnotných vzorů (MATLAB, Octave) a máme tu řešení:

// Numerický výpočet určitého integrálu funkce f(x) na intervalu <a,b>
// Lichoběžníkové pravidlo a) na pevném počtu uzlů, b) s předepsanou přesností

program Lichobeznik;

type
  tF = function (const x : real) : real; // typ pro proceduralni argument

// polynom
function fpol(const x : real) : real;
begin
  result:=x*x;
end;

// goniometricka funkce
function fgon(const x : real) : real;
begin
  result:=sin(x);
end;

// integrator zalozeny na kvadraturnim vzorci na pevnem poctu uzlu
function quad(f : tF; const a,b : real; const nmax : integer) : real; overload;
var
  x,h : real;
  n : integer;
begin
  h:=(b-a)/nmax;
  result:=(f(a)+f(b))*0.5;
  x:=a;
  for n:=1 to nmax-1 do begin
    x:=x+h;
    result:=result+f(x);
  end;
  result:=result*h;
end;

// kvadraturni integrator s predepsanou presnosti
function quad(f : tF; const a,b : real; const eps : real) : real; overload;
const
  nmax = 1000000000;      // max. pocet uzlu
var
  L,Lold : real;
  n : integer;
begin
  n:=1;
  L:=-1e300;              // vychozi hodnota vysledku
  repeat
    Lold:=L;
    L:=quad(f,a,b,n);     // kvadraturni vzorec na pevnem poctu uzlu
    n:=n*2;               // zdvojnasobeni poctu uzlu
    if n>nmax then break; // zastavovaci podminka - pocet uzlu
  until abs(L-Lold)<eps;  // zastavovaci podminka - presnost
  result:=L;
end;

begin
  writeln('Integral: ',quad(@fpol,0,1,10000):18:15);  // presne: 1/3
  writeln('Integral: ',quad(@fgon,0,pi,10000):18:15); // presne: 2
  writeln('Integral: ',quad(@fpol,0,1,1e-10):18:15);
  writeln('Integral: ',quad(@fgon,0,pi,1e-10):18:15);
  readln;
end.

Bylo by odpovědné si uvědomit, že \(L_N\) lze při zahrnutí \(L_{N/2}\) vyčíslit s poloviční námahou než při vyčíslení na zelené louce (tak, jak to marnotratně dělá ukázka). Zvýrazněme ještě, že místo /2 používáme raději *0.5, neboť dnešní procesory násobí rychleji než dělí, ovšem překladače mohou tuto záměnu udělat za nás. Nepřehlédněme adresní operátor @ u procedurálního argumentu ve volání funkce quad, v Delphi byl nepovinný, Free Pascal nám ho neodpustí.

Domácí úkol bude mít dvě různě pracné části. Tou snazší je doplnit tuto ukázku o výpočet prvního integrálu z druhé stránky PDF se vzorci pro \(\pi\),

\[\pi=\int_0^1\frac{16x-16}{x^4-2x^3+4x-4}\,dx\,.\]

Newtonova metoda tečen

Chceme numericky nalézt jeden z nulových bodů funkce \(f(x)\). Metod je řada, jednou z výchozích je Newtonova metoda tečen (anglicky Newton-Raphson method). V bodě \(x_0\) poblíž nulového bodu (bod \(x_0\) je třeba zjistit předběžnou, třeba grafickou analýzou funkce) sestrojíme tečnu k \(f(x)\), \(y(x)=f(x_0)+f'(x_0)(x-x_0)\), nalezneme analyticky její nulový bod a ten ztotožníme s novým odhadem \(x_1\) nulového bodu \(f(x)\). Jde tedy o iterační proces:

\[x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)},\quad n=0,1,\dots\]

Metoda si žádá nejen funkční hodnotu, ale i hodnotu první derivace. Je pěkné, je-li pro derivaci dostupný analytický vzorec, získávat ji lze ovšem i numericky (jde pak o metodu sečen). Vnímáme, že derivace je ve jmenovateli a že je-li blízká nule, vzorec nás odvede k průsečíku téměř vodorovné tečny s vodorovnou osou, tedy patrně někam daleko od našeho cíle. (Metoda nemusí konvergovat k žádnému nulovému bodu.)

// Numerické hledání nulového bodu funkce f(x) v okolí bodu x0
// Newtonova metoda tečen
// f(x) reálná funkce reálné proměnné

program NewtonR;

type
  tA = array [0..1] of real;           // pro 0. a 1. derivaci
  tF = function (const x : real) : tA; // pro proceduralni argument

// polynom
function fpol(const x : real) : tA;
begin
  result[0]:=x*x*x-1;
  result[1]:=3*x*x; // analyticka 1. derivace
end;

// goniometricka funkce
function fgon(const x : real) : tA;
begin
  result[0]:=sin(x);
  result[1]:=cos(x); // analyticka 1. derivace
end;

// Newtonuv hledac korenu (root finder)
function newton(f : tF; const x0,eps : real) : real;
const
  nmax=20;                        // max. pocet iteraci
var
  ff : tA;
  x,xold : real;
  n : integer;
begin
  n:=0;
  x:=x0;                          // vychozi bod
  repeat
    xold:=x;
    ff:=f(x);                     // vycisleni funkce a jeji derivace
    x:=xold-ff[0]/ff[1];          // Newtonuv krok
    inc(n); if n>nmax then break; // zastavovaci podminka - pocet iteraci
  until abs(x-xold)<eps;          // zastavovaci podminka - presnost
  result:=x;
end;

begin
  writeln('Nulovy bod: ',newton(@fpol,2,1e-10):18:15); // presne: 1
  writeln('Nulovy bod: ',newton(@fgon,3,1e-10):18:15); // presne: pi
  readln;
end.

Je také pěkné, že i komplexní funkce komplexní proměnné mají tečny a že Newtonovu metodu lze tedy zobecnit pro tyto funkce:

// Numerické hledání nulového bodu funkce f(x) v okolí bodu x0
// Newtonova metoda tečen
// f(x) komplexní funkce komplexní proměnné

program NewtonC;
uses uMyComplex; // nase komplexni aritmetika

type
  tA = array [0..1] of complex;           // pro 0. a 1. derivaci
  tF = function (const x : complex) : tA; // pro proceduralni argument

// polynom
function fpol(const x : complex) : tA;
begin
  result[0]:=x*x*x-1;
  result[1]:=x*x*3; // analyticka prvni derivace
end;

// Newtonuv hledac korenu (root finder)
function newton(f : tF; const x0 : complex; const eps : real) : complex;
const
  nmax=20;                        // max. pocet iteraci
var
  ff : tA;
  x,xold : complex;
  n : integer;
begin
  n:=0;
  x:=x0;                          // vychozi bod
  repeat
    xold:=x;
    ff:=f(x);                     // vycisleni funkce a jeji derivace
    x:=xold-ff[0]/ff[1];          // Newtonuv krok
    inc(n); if n>nmax then break; // zastavovaci podminka - pocet iteraci
  until abs(x-xold)<eps;          // zastavovaci podminka - presnost
  result:=x;
end;

var
  x : complex;

begin
  x:=newton(@fpol,Cmplx(2,0),1e-10);
  writeln('Nulovy bod: ',aReal(x):18:15,Imag(x):18:15);
  readln;
end.

Ukázka odkazuje na modul uMyComplex s naší implementací komplexní aritmetiky, který jsme nachystali na jednom z předchozích cvičení. Tam jsme ovšem nepotřebovali komplexní dělení, které potřebujeme teď. Zvolíme tuto cestu: \(x/y=x.(1/y)\), kde \(1/(a+ib)=(a-ib)/(a^2+b^2)\). Modul uMyComplex tedy doplníme o toto:

interface
function cDiv(const x,y : complex) : complex; overload;
operator / (const x,y : complex) : complex; overload;
function cDiv(const x : complex; const y : real) : complex; overload;
operator / (const x : complex; const y : real) : complex; overload;

implementation
function cInv(const x : complex) : complex;
var
  factor : extended;
begin
  factor:=1/(sqr(x.re)+sqr(x.im));
  with result do begin re:=x.re*factor; im:=-x.im*factor end;
end;

function cDiv(const x,y : complex) : complex; overload;
begin
  result:=cMul(x,cInv(y));
end;

operator / (const x,y : complex) : complex; overload;
begin
  result:=cDiv(x,y);
end;

function cDiv(const x : complex; const y : real) : complex; overload;
begin
  with result do begin re:=x.re/y; im:=x.im/y end;
end;

operator / (const x : complex; const y : real) : complex; overload;
begin
  result:=cDiv(x,y);
end;

Modul bychom měli ještě doplnit o implementaci odčítání pro dvojici complex a real a také o přetížení jména abs pro complex argument sdružením s funkcí cAbs. Když to všechno dáme dohromady, vznikne zip s modulem uMyComplex se 3 konverzními funkcemi, 4 aritmetickými operátory a funkcí abs.

Domácí úkol (jeho druhá část) zahrnuje sloučení systematického průchodu obdélníkem v komplexní rovině, obdobně jako v předchozím úkolu s Mandelbrotovou množinou, a Newtonovy metody pro komplexní funkci. Zjistěte, jakou část komplexní roviny zabírají body, z nichž startovaná Newtonova metoda pro komplexní polynom \(f(x)=x^3-1\) konverguje k jednomu ze tří kořenů, \(x_j=e^{(2/3)ij\pi}\), \(j=0,1,2\). (Např. \(x_0=1\).) Můžete (nemusíte) sledovat číslo iterace, po které se Newtonova metoda přiblíží k vybranému kořenu na vzdálenost menší než \(\varepsilon\), a volit podle čísla této iterace barvu výchozího bodu v obrázku.

Výstupní data (dva sloupce souřadnic a třetí sloupec s indexem iterace) můžete vypisovat ASCII nebo binárně a vykreslit např. tímto skriptem v gnuplotu:

# gnuplot: vykreslení fraktálu
r=1.5            # polomer dat
xsize=800        # rozliseni ve smeru x
ysize=xsize      # rozliseni ve smeru y
binary=0         # typ vstupnich souboru: 0 ascii, 1 binary
set size square
set palette defined (  0     'dark-blue',   2 'dark-blue',  10 'blue',  15 'gold', 50 'brown',\
                      50.001 'black', 51 'black')
set cbrange [0:51]
unset border
unset xtics
unset ytics
unset colorbox
unset key
set term png enhanced truecolor size xsize,ysize
set output 'newton.png'
if (!binary) \
plot [-r:r][-r:r] 'newton.dat' using ($1**2+$2**2<r**2?$1:1/0):2:3 palette z with dots; \
else \
plot [-r:r][-r:r] 'newton.bin' binary format='%2f%i' \
                               using ($1**2+$2**2<r**2?$1:1/0):2:3 palette z with dots
set output

Ve skriptu můžeme zaměřit tyto dosud nekomentované postupy:

  • Vypínání vlastností, zde os včetně jejich popisů, sloupce s barevnou škálou a legendy, se provádí pomocí příkazu unset.
  • Binární čtení reálných čísel: příslušný formát pro 4bytový typ (Pascal: single) zní %f, pro 2 sloupce %2f, pro 8bytový typ (Pascal: real) je to %lf – zde ovšem můžeme narazit na zádrhel, kterým se nyní nezabývejme a raději reálná čísla v Pascalu binárně zapisujme jako 4bytová; deklarace výstupního souboru by pak mohla znít podobně jako tato: file of record x,y : single; n : integer end.
  • Analytické přepočty dat ze souborů: to je schopnost mimořádné síly. Hodnoty ve sloupci datového souboru lze v klauzuli using příkazů plot a splot vyjádřit symbolicky jako (zde pro 1. sloupec) $1 a použít v jakémkoliv výrazu, jakého je gnuplot schopen; to celé se ještě obklopí závorkami. Místo obyčejného using 1:2 můžeme tedy napsat ekvivalentní using ($1):($2) a v závorkách se pak dále rozmáchnout. To jsme udělali i v příkazu plot výše, kde jsme na místě prvního sloupce pomocí ternárního podmíněného operátoru ?: (Wikipedia) podmínili výstup jen pro body, které splňují podmínku v prvním operandu (leží v kruhu o poloměru r). Pro ostatní body se totiž nemá použít údaj $1, ale výraz 1/0, což je gnuplotovský způsob, jak říci, že bod nemá být vykreslen.

Obrázek malovaný níže, vytvořený skriptem uvedeným výše, zachycuje rychlost konvergence Newtonovy metody ke každému ze 3 kořenů polynomu \(x^3-1\). Barevné odlišení oblastí příslušných jednotlivým kořenům se mi nelíbilo, proto tam není, a kruhový formát obrázku připraveného z dat na obdélníku jsem zvolil jednak pro ukázku schopností gnuplotu, jednak pro zvýraznění invariance fraktálu vůči otočení kolem počátku o násobky \(\frac23\pi\).

_images/newton.png

Přikládám týž obrázek v rozlišení (10k)2 bodů pro podrobnější vhled. O fraktálech tohoto typu se na Wikipedii píše v heslech Newtonova metoda, Atraktor a Juliova množina.

Cvičení č. 7

Pokračujeme s numerickými metodami: budeme programovat výpočet určitého integrálu metodou Monte Carlo a obrátíme se k algoritmům lineární algebry – nachystáme modul s operátory pro maticové násobení a řešení soustav lineárních algebraických rovnic.

Integrování Monte Carlo

Pro numerický výpočet vícerozměrného určitého integrálu lze použít vzorec

\[\int_V f\,dV \approx V\langle f\rangle \pm V\sqrt{\frac{\langle f^2\rangle-\langle f\rangle^2}N}\,,\]

kde \(V\) je integrační oblast a úhlové závorky značí aritmetický průměr, \(\langle f\rangle=\frac1N\sum_{i=1}^N f_i,\ f_i=f(x_i)\). Jméno metodě Monte Carlo dává náhodný výběr bodů \(x_i\), v kterých se integrand vyčísluje. Chybový člen je třeba chápat ve statistickém smyslu; jasným poselstvím je však pomalý pokles chyby s \(1/\sqrt{N}\). Pro srovnání, chyba lichoběžníkového pravidla klesá s \(1/N^2\). Metoda Monte Carlo tak obvykle nesoutěží s kvadraturními vzorci a používá se až za hranicí jejich aplikovatelnosti. My si však zasoutěžíme a Monte Carlo pro 1D případ použijeme:

\[\int_a^b f(x)\,dx \approx (b-a)\langle f\rangle = \frac{b-a}N \sum_{i=1}^N f_i\,.\]

Použijeme přímočaře tento vzorec, a použijeme i 2D variantu s oblíbenou geometrickou interpretací: (1D) integrand pokryjeme (2D) obdélníkem \(\langle a,b\rangle\times\langle y_{\min},y_{\max}\rangle\) tak, aby \(y_{\min}\le f(x)\le y_{\max}\), v obdélníku vybereme \(N\) náhodných bodů \((x_i,y_i)\) a zjistíme počet \(K\) těch z nich, pro které platí \(y_i\le f(x_i)\). Je pak totiž

\[\int_a^b f(x)\,dx \approx \frac{SK}N+(b-a)y_{\min}, \quad S=(b-a)(y_{\max}-y_{\min})\,,\]

kde první sčítanec je plochou obdélníka sníženou úměrně relativnímu počtu “zásahů” pod integrand. Odvození tohoto vzorce ze vzorce pro vícerozměrný integrál je v poznámkách k přednášce.

// Numerický výpočet určitého integrálu funkce f(x) na intervalu <a,b>
// Metoda Monte Carlo v 1D a 2D variantě

program MonteCarlo;

type
  tF = function (const x : real) : real; // typ pro procedurální argument

// integrand
function f(const x : real) : real;
begin
  result:=(16*x-16)/(((x-2)*x*x+4)*x-4);
end;

// 1D Monte Carlo integrátor
function MonteCarlo(f : tF; const a,b : real; const nmax : integer) : real; overload;
var
  x,s : real;
  n : integer;
begin
  s:=0;
  randomize; // inicializace posloupnosti náhodných čísel
  for n:=1 to nmax do begin
    x:=a+(b-a)*random;
    s:=s+f(x); // střední hodnota <f>
  end;
  result:=(b-a)*s/nmax; // (b-a)*<f>
end;

// 2D Monte Carlo integrátor
function MonteCarlo(f : tF; const a,b,ymin,ymax : real; const nmax : integer) : real; overload;
var
  x,y : real;
  n,k : integer;
begin
  k:=0;
  randomize;
  for n:=1 to nmax do begin
    x:=a+(b-a)*random;
    y:=ymin+(ymax-ymin)*random;
    if y<=f(x) then k:=k+1; // leží bod pod integrandem?
  end;
  result:=(b-a)*(ymax-ymin)*k/nmax+(b-a)*ymin; // plocha*K/N + plocha mezi x-osou a obdélníkem
end;

begin
  writeln(MonteCarlo(@f,0,1,1000000));
  writeln(MonteCarlo(@f,0,1,0,4,1000000));
  readln;
end.

Jako integrand jsme volili racionální funkci z domácího úkolu k minulému cvičení. Dali jsme si záležet a k vyčíslení polynomu ve jmenovateli použili vytýkání před závorku alias metodu zvanou Hornerovo schéma, jejíž časová složitost je lineární vzhledem k stupni polynomu (tedy úspornější, než kdybychom vyčíslovali x*x, x*x*x a x*x*x*x navzájem nezávisle). Co má při integraci v intervalu \(\langle0,1\rangle\) vyjít, víme od minula, ostatně můžeme si to od oka ověřit v gnuplotu:

# gnuplot: plot with boxes
f(x)=(16*x-16)/(((x-2)*x*x+4)*x-4)
set samples 1000
set style fill solid
plot [0:1] f(x) with boxes linecolor 'green'
_images/pi_int.png

Nepohoršujme se, že pro \(N=10^6\) bodů získá metoda Monte Carlo pouze 3 platná místa výsledku, o její \(1/\sqrt{N}\) chybě víme od začátku.

Lineární algebra

Jsme u kapitoly č. 1 učebnic numerické matematiky. Maticové násobení,

\[c_i=\sum_k A_{ik}b_k,\quad C_{ij}=\sum_k A_{ik}B_{kj}\,,\]

se v nich sice nutně neobjevuje, to je (samozřejmě jen na první pohled) příliš jednoduché, ale řešení soustav lineárních algebraických rovnic ano,

\[\sum_j A_{ij}x_j=b_i,\quad \sum_j A_{ij}X_{jm}=B_{im}\,.\]

My si pro obě tyto operace přetížíme aritmetické operátory * a /, podobně jako to dělá Octave (ten volá řešič soustav operátorem \, který je v Pascalu nedostupný). Hvězdičku si naprogramujeme sami pro varianty matice krát matice a matice krát vektor, pro lomítko si pomůžeme procedurou gaussjNumerických receptů, popsanou v kapitole 2.1 knihy ve variantě pro C i pro Fortran. (Knihu pro Pascal nevystavili veřejně, ale to ani nevadí.) Procedura provádí Gaussovu eliminaci čtvercové matice uložené ve statickém poli pro několik pravých stran, uložených po sloupcích rovněž ve statické matici. Náš modul uNRGauss jejich proceduru obaluje implementací operátoru, který jako levý operand očekává dynamickou čtvercovou matici a jako pravý operand buď dynamický vektor jedné pravé strany nebo matici několika pravých stran. Pro vektor je v modulu připraven datový typ tVector, pro matici tMatrix, a vůbec – rozhraní modulu vypadá takto:

// Lineární algebra
// operátor * pro maticové násobení: matice*matice, matice*vektor
// operátor / pro řešení soustavy lineárních rovnic: A/b, A/B

unit uNRGauss;

interface

const
  np=1000; // max. počet rovnic
  mp=1000; // max. počet pravých stran
type
  tMatrix = array of array of real;
  tVector = array of real;
  operator * (const A : tMatrix; const b : tVector) : tVector;
  operator * (const A : tMatrix; const B : tMatrix) : tMatrix;
  operator / (const A : tMatrix; const b : tVector) : tVector;
  operator / (const A : tMatrix; const B : tMatrix) : tMatrix;

implementation
// ...

Vypíšeme i část s implementací operátorů:

operator * (const A : tMatrix; const b : tVector) : tVector;
var
  s : real;
  i,k : integer;
begin
  SetLength(result,Length(A));
  for i:=0 to High(A) do begin
    s:=0;
    for k:=0 to High(A[0]) do s:=s+A[i,k]*b[k];
    result[i]:=s;
  end;
end;

operator * (const A : tMatrix; const B : tMatrix) : tMatrix;
var
  s : real;
  i,j,k : integer;
begin
  SetLength(result,Length(A),Length(B[0]));
  for i:=0 to High(A) do
  for j:=0 to High(B[0]) do begin
    s:=0;
    for k:=0 to High(A[0]) do s:=s+A[i,k]*B[k,j];
    result[i,j]:=s;
  end;
end;

operator / (const A : tMatrix; const b : tVector) : tVector;
var
  AA : RealArrayNPbyNP;
  bb : RealArrayNPbyMP;
  n,i,j : integer;
begin
  n:=Length(b);
  SetLength(result,n);
  for i:=1 to n do for j:=1 to n do
    AA[i,j]:=A[i-1,j-1];
  for i:=1 to n do
    bb[i,1]:=b[i-1];
  gaussj(AA,n,bb,1);
  for i:=1 to n do
    result[i-1]:=bb[i,1];
end;

operator / (const A : tMatrix; const B : tMatrix) : tMatrix;
var
  AA : RealArrayNPbyNP;
  BB : RealArrayNPbyMP;
  n,m,i,j : integer;
begin
  n:=Length(B);
  m:=Length(B[0]);
  SetLength(result,n,m);
  for i:=1 to n do for j:=1 to n do
    AA[i,j]:=A[i-1,j-1];
  for i:=1 to n do for j:=1 to m do
    BB[i,j]:=B[i-1,j-1];
  gaussj(AA,n,BB,m);
  for i:=1 to n do for j:=1 to m do
    result[i-1,j-1]:=BB[i,j];
end;

Vidíme, že má-li funkce vracet dynamické pole, je třeba pro něj alokovat paměťový prostor procedurou SetLength(result,...), že se řádně vyptáváme na velikost první dimenze matice dotazem Length(A) a na velikost druhé dimenze dotazem Length(A[0]), že dynamická pole jsou v Pascalu indexována od nuly (zatímco statická pole v proceduře gaussj od 1) a že lehkomyslně nekontrolujeme shodu velikostí dimenzí polí požadovaných pro smysluplnost operací. Celý modul uNRGauss je zde.

A je vhodná chvíle pro test modulu. Program nachystá dynamickou matici A náhodných čísel, vektor b s jednotkovými prvky a operátorem / nalezne řešení x soustavy \(A.x=b\). To se vypíše a pak vynásobí s maticí soustavy, \(A.x\), abychom si ověřili, zda jsme po tom všem zpátky u původního vektoru b.

// Řešení soustavy lineárních algebraických rovnic pomocí operátoru /
program Gauss;
uses uNRGauss;
var
  A : array of array of real;
  b,x : array of real;
  n,i,j : integer;
begin
  n:=100;      // počet rovnic
  if n>np then begin writeln('Error: n > ',np); readln; halt end;
  SetLength(A,n,n); SetLength(b,n); setLength(x,n);
  randomize;
  for i:=1 to n do begin
    for j:=1 to n do A[i-1,j-1]:=random; // náhodná matice
    b[i-1]:=1; // vektor pravé strany
  end;
  x:=A/b;      // operátor pro řešení soustavy algebraických rovnic
  for i:=1 to n do writeln(i:4,x[i-1]:15:10);
  b:=A*x;      // operátor pro maticové násobení (matice*vektor)
  for i:=1 to n do writeln(i:4,b[i-1]:15:10);
  readln;
end.

Domácí úkol: úpravou předchozí ukázky připravte program, který pro čtvercové Hilbertovy matice \(H_{ij}=1/(i+j-1),\) \(i,j=1,\dots,N,\) spočítá inverzní matice a pro každou vypíše součet jejích prvků, \(S(N)=\sum_{i,j=1}^N (H^{-1})_{ij}\). To vše podstupte pro počet řádků \(N\) od 1 do 12 a uhodněte, co je \(S(N)\). Výše to s 8bytovým reálným typem přímočaře nejde, Hilbertovy matice se numericky těžko invertují. Inverzní matici \(A^{-1}\) lze získat řešením soustavy \(AX=I\), kde \(I\) je jednotková matice. Z předchozí četby už víte, že operátor / je nachystán i pro takovou rovnici. Nekazte si příliš brzy radost z poznávání řádkem, kterým by se úkol řešil v Octavu: for n=1:12, disp(sprintf('%5d %20.15f',n,sum(sum(inv(hilb(n)))))), end.

Cvičení č. 8

Následují další numerické metody: budeme prokládat aproximační polynom tabelovanými body. Budeme přitom usilovat o vystižení bodů přesně (polynomická interpolace) a přibližně (polynomická regrese).

Polynomická interpolace

Vstupními daty budou dvojice souřadnic bodů v rovině, \((x_i,f_i),\ i=0,\dots,n\). Cílem polynomické interpolace je sestrojit aproximační polynom proložený vstupními body přesně. Pro \(n+1\) bodů to bude polynom \(n\)-tého stupně,

\[P_n(x)=c_0+c_1x+\dots+c_nx^n\,.\]

Chceme tedy dosáhnout toho, aby \(P_n(x_i)=f_i\), což je vlastně \(n+1\) lineárních algebraických rovnic pro \(n+1\) neznámých koeficientů \(c_0,\dots,c_n\) aproximačního polynomu. Maticově:

\[\begin{split}\left(\begin{array}{cccc} 1 & x_0 & \dots & x_0^n\\ 1 & x_1 & \dots & x_1^n\\ \vdots & & & \vdots\\ 1 & x_n & \dots & x_n^n \end{array}\right) \left(\begin{array}{c} c_0 \\ c_1 \\ \vdots \\ c_n \end{array}\right)= \left(\begin{array}{c} f_0 \\ f_1 \\ \vdots \\ f_n \end{array}\right).\end{split}\]

Pro řešení algebraické soustavy máme od minula přichystaný modul s operátorem /. Algoritmus polynomické interpolace zabalíme do funkce polyfit(x,f), jejíž jméno opisujeme od funkce pro polynomickou fitaci v Octavu. Tomu přizpůsobíme i řazení koeficientů polynomu ve výstupním poli: od \(c_n\)\(c_0\). Počet bodů ve vstupním souboru s daty necháme spočtením řádků souboru zjistit automaticky, pro data pak alokujeme dynamická pole na míru. Po spočtení koeficientů aproximačního polynomu jej vyčíslíme na husté síti překrývající interval s daty, samozřejmě užitím Hornerova schématu (známe z minulého cvičení), funkcí polyval pojmenovanou opět podle vzoru z Octavu. Souřadnice polynomu zapíšeme dvousloupcově do výstupního souboru. Program s polynomickou interpolací a další metodou je přiložen o několik řádků níže.

Polynomická regrese

Vstupní data často nejsou natolik přesná, aby bylo třeba nutit aproximační funkci jimi procházet. Vyřešíme nyní frekventovaný požadavek proložit \(n+1\) body aproximační polynom stupně nižšího než \(n\), tedy vystihnout body jen přibližně. Cílem polynomické regrese je sestrojit v duchu metody nejmenších čtverců takový polynom \(P_m(x)\) stupně \(m\), aby součet druhých mocnin jeho odchylek od vstupních bodů byl minimální,

\[\sum_{k=0}^n\left[P_m(x_k)-f_k\right]^2 = \min\,.\]

V poznámkách k přednášce je tento vztah zderivován vzhledem k \(m+1\) koeficientům aproximačního polynomu a odvozena tak soustava \(m+1\) lineárních algebraických rovnic pro hledané koeficienty:

\[\begin{split}\left(\begin{array}{llll} \sum_k 1 & \sum_k x_k & \dots & \sum_k x_k^m\\ \sum_k x_k & \sum_k x_k^2 & \dots & \sum_k x_k^{m+1}\\ \vdots & & & \vdots\\ \sum_k x_k^m & \sum_k x_k^{m+1} & \dots & \sum_k x_k^{2m} \end{array}\right) \left(\begin{array}{l} c_0 \\ c_1 \\ \vdots \\ c_m \end{array}\right)= \left(\begin{array}{l} \sum_k f_k \\ \sum_k f_kx_k \\ \vdots \\ \sum_k f_kx_k^m \end{array}\right) .\end{split}\]

Algoritmu dáme opět podobu funkce polyfit(x,f,m), tentokrát tříargumentové.

// Polynomická aproximace: polynomická interpolace a polynomická regrese
// Implementace MATLAB/Octave funkcí
// p=polyfit(x,f) pro fitaci dat (x_i,f_i), i=0,..,n, polynomem stupně n
// p=polyfit(x,f,m) pro fitaci dat (x_i,f_i), i=0,..,n, polynomem stupně m (LSQ)
// f=polyval(p,x) pro vyčíslení polynomu p ve skaláru x nebo na poli x (Horner)
// polynom p: 1D vektor počínaje koeficientem u členu s nejvyšší mocninou
// Vstupní data: soubor 'data', výstupní data: soubor 'fit'

program Fit;
uses uNRGauss;
const
  m = 1; // stupeň aproximačního polynomu

// Přesná polynomická fitace (polynomická interpolace)
function polyfit(const x,f : tVector) : tVector; overload;
var
  V : tMatrix;
  c : tVector;
  n,i,j : integer;
begin
  n:=High(x);
  SetLength(result,n+1); SetLength(V,n+1,n+1);
  for i:=0 to n do V[i,0]:=1; // Vandermondova matice
  for j:=1 to n do begin
    for i:=0 to n do V[i,j]:=V[i,j-1]*x[i];
  end;
  c:=V/f; // řešení soustavy rovnic
  for i:=0 to n do result[i]:=c[n-i]; // změna pořadí (flipud)
end;

// Přibližná polynomická fitace (polynomická regrese)
function polyfit(const x,f : tVector; const m : integer) : tVector;
var
  G : tMatrix;
  b,c,xx,sumx,sumf : tVector;
  n,i,j,k : integer;
  s : real;
begin
  SetLength(result,m+1); SetLength(G,m+1,m+1); SetLength(b,m+1);
  n:=High(x); SetLength(xx,n+1); SetLength(sumx,m*2+1); SetLength(sumf,m+1);
  for k:=0 to n do xx[k]:=1;
  for i:=0 to m*2 do begin
    s:=0; for k:=0 to n do s:=s+xx[k]; sumx[i]:=s;
    if i<=m then begin s:=0; for k:=0 to n do s:=s+f[k]*xx[k]; sumf[i]:=s end;
    for k:=0 to n do xx[k]:=xx[k]*x[k];
  end;
  for i:=0 to m do for j:=0 to m do G[i,j]:=sumx[i+j]; // Gramova matice
  b:=sumf; // vektor pravé strany
  c:=G/b; // řešení soustavy rovnic
  for i:=0 to m do result[i]:=c[m-i]; // změna pořadí (flipud)
end;

// Vyčíslení polynomu (Hornerovo schéma) v bodu x
function polyval(const p : tVector; const x : real) : real; overload;
var
  i : integer;
begin
  result:=p[0];
  for i:=1 to High(p) do result:=result*x+p[i];
end;

// Vyčíslení polynomu (Hornerovo schéma) v poli bodů x
function polyval(const p : tVector; const x : tVector) : tVector; overload;
var
  i : integer;
begin
  SetLength(result,Length(x));
  for i:=0 to High(x) do result[i]:=polyval(p,x[i]);
end;

var
  x,f,p,xi,fi : tVector;
  n,i,nmax : integer;
  xmin,xmax : real;
  t : text;
  input : string = 'data';
  output : string = 'fit';
begin
// čtení vstupních dat
  n:=0;
  assign(t,input); reset(t); while not eof(t) do begin readln(t); n:=n+1 end; // eof: end-of-file
  n:=n-1;
  SetLength(x,n+1); SetLength(f,n+1); SetLength(p,m+1);
  reset(t); for i:=0 to n do readln(t,x[i],f[i]); close(t);
  writeln('Soubor ',input,' - precteno radku: ',n);
// výpočet
  if m=n then p:=polyfit(x,f)    // polynomická interpolace
         else p:=polyfit(x,f,m); // polynomická regrese
// tabelování aproximačního polynomu
  nmax:=n*100; xmin:=x[0]-1; xmax:=x[n]+1;
  SetLength(xi,nmax+1); SetLength(fi,nmax+1);
  for i:=0 to nmax do xi[i]:=xmin+(xmax-xmin)*i/nmax;
  fi:=polyval(p,xi);
  assign(t,output); rewrite(t);
  for i:=0 to nmax do writeln(t,xi[i]:10:6,' ',fi[i]:10:6);
  close(t);
  writeln('Soubor ',output,' - zapsano radku: ',nmax);
  readln;
end.

Domácím úkolem je spočítat pro několik (\(n+1\), třeba \(n=3\)) Vámi zvolených bodů koeficienty aproximačních polynomů stupňů \(0,1,\dots,n\), polynomy vyčíslit na husté síti a skreslit je všechny spolu se vstupními body do jednoho obrázku, podobného tomuto:

_images/fit.png

Obrázek těží ze schopnosti gnuplotu (příkaz fit) počítat metodou nejmenších čtverců parametry předepsané funkce aproximující data v souboru a klidně tak odstavit (nebo lépe: ověřit) naše snahy v Pascalu:

# gnuplot: polynomická aproximace
f0(x)=a0                          # konstantní funkce
f1(x)=b1*x+b0                     # polynom 1. stupně
f2(x)=(c2*x+c1)*x+c0              # polynom 2. stupně
f3(x)=((d3*x+d2)*x+d1)*x+d0       # polynom 3. stupně
fit f0(x) 'data' via a0           # metoda nejmenších čtverců
fit f1(x) 'data' via b1,b0
fit f2(x) 'data' via c2,c1,c0
fit f3(x) 'data' via d3,d2,d1,d0
set key top left
set samples 30
plot [0:5][-1:5] 'data' with points linestyle 7 pointsize 3,\
     'fit0' with lines linestyle 1,f0(x) with points pointtype 7 linecolor 1,\
     'fit1' w l ls 2,f1(x) w p pt 7 lc 2,\
     'fit2' w l ls 3,f2(x) w p pt 7 lc 3,\
     'fit3' w l ls 4,f3(x) w p pt 7 lc 4
set terminal pngc enhanced size 1024,768
set output 'fit.png'
replot
set output
unset terminal

Implementace (nelineární) metody nejmenších čtverců v gnuplotu je postavena na obecnějším algoritmu, než jsme předvedli zde, aproximační funkce nemusí být jen polynomem.

Cvičení č. 9

Předvánoční cvičení věnujeme dekódování koledy z Morseovy abecedy a symbolického notového zápisu.

Morseova abeceda

Smyslem je seznámit se s programováním dynamických datových struktur, konkrétně binárního stromu (binary tree). Jeden binární strom dvěma způsoby připravíme, dvakrát celý projdeme a jednou použijeme pro dichotomické hledání. Učebnicovou aplikací je dekódování zpráv v Morseově abecedě – i na Wikipedii je obrázek binárního stromu s Morseovkou a podobný obrázek z gnuplotu (data a skript zde) máme i my:

_images/BTree.png

(Cestu vlevo podstoupíme pro tečku, cestu vpravo pro čárku: pro Morseův kód -.-. tak putováním od kořenu stromu vpravo-vlevo-vpravo-vlevo nalezneme ASCII znak C.) Elegantní dekódování pomocí binárního stromu pak porovnáme s prostým dekódováním pomocí sekvenčního procházení seznamu.

Prvním krokem je příprava datového typu pro uzel binárního stromu. Má obsahovat datovou položku a dva ukazatele na další uzly, bude to tedy záznam:

type
  pBTree = ^tBTree;
  tBTree = record
    A : string;     // datová položka
    pT,pC : pBTree; // ukazatele v tečkovém a čárkovém směru
  end;

Deklarovat typ pBTree pro ukazatel na typ tBTree, jehož deklarace je až za deklarací typu pBTree, je možné, jinak by tyto datové struktury ani programovat nešlo. Pro datovou položku potřebujeme delší než 1znakový typ char (kvůli spřežce Ch), použijeme string. Založení binárního stromu pak můžeme provést třeba takto, uzel po uzlu:

var
  root,p : pBTree;
begin
// alokace a inicializace kořenu
  New(root); root^.A:='@'; root^.pT:=nil; root^.pC:=nil;
// alokace a inicializace uzlů první hladiny
  New(root^.pT); with root^.pT^ do begin A:='E'; pT:=nil; pC:=nil end;
  New(p); with p^ do begin A:='T'; pT:=nil; pC:=nil end; root^.pC:=p;
// alokace a inicializace uzlů druhé hladiny
  New(p); with p^ do begin A:='I'; pT:=nil; pC:=nil end; root^.pT^.pT:=p;
  New(p); with p^ do begin A:='A'; pT:=nil; pC:=nil end; root^.pT^.pC:=p;
  New(p); with p^ do begin A:='N'; pT:=nil; pC:=nil end; root^.pC^.pT:=p;
  New(p); with p^ do begin A:='M'; pT:=nil; pC:=nil end; root^.pC^.pC:=p;
// výpis uzlů průchodem do šířky (po hladinách)
  writeln(root^.A,root^.pT^.A,root^.pC^.A,root^.pT^.pT^.A,root^.pT^.pC^.A,root^.pC^.pT^.A,root^.pC^.pC^.A);

Procedurou New se alokuje prostor pro cíl, na který má ukazovat ukazatel v argumentu, a ukazatel se tam nasměruje. Vhodné je položky v nově vytvořeném uzlu bezprostředně inicializovat, především nulovat ukazatele (přiřazením nil). Přístup k položkám cílového záznamu se děje prostřednictvím dereference ukazatele vyznačené stříškou, root^. K hromadné inicializaci položek záznamu se používá příkaz with. Nový uzel můžeme vytvářet rovnou z rodičovského uzlu, např. New(root^.pT), nebo užitím pomocného ukazatele, New(p), který nakonec na patřičnou pozici v rodičovském uzlu přeneseme, root^.pC:=p. Výše jsme založili dvě hladiny binárního stromu a není pracné sepsat inicializační příkazy pro všech asi 40 uzlů v tomto duchu, je to přehledné copy-and-paste. Nemůžeme si ale odpustit ani přípravu procedury pro inicializaci binárního stromu z dvousloupcové matice obsahující ASCII znaky a jejich Morseovy kódy. Níže náš čeká odkaz na modul pro dekódování Morseovky, který obsahuje obě inicializační procedury.

Kontrolní výpis datových položek ze dvou hladin binárního stromu tak, jak jsme ho provedli výše, je na hranici čitelnosti nebo aspoň trpělivosti čtenáře. Navíc by jistě bylo žádoucí, abychom uměli vypsat obsah všech uzlů stromu, aniž bychom předem znali jeho kompletní strukturu. K tomu slouží dvojice algoritmů, z nichž prvnímu se říká prohledávání binárního stromu do šířky (breadth-first search). Průchod stromu prohledáváním do šířky jsme pro kořen a první dvě hladiny vlastně výše napodobili. Obecně jde o to, že se nenulové ukazatele z každého uzlu počínaje kořenem průběžně umísťují na konec fronty, z jejíhož začátku se postupně odebírají a zpracovávají. Chování ukazatelů ve frontě je nekonfliktní, popisuje se frází First-In-First-Out (FIFO). Konkrétně: umístíme ukazatel na kořen do fronty a pak v cyklu ukazatel z fronty vyzvedneme, data jeho cíle zpracujeme a ukazatele z cíle (jsou-li nenulové) umístíme do fronty, to vše tak dlouho, dokud ve frontě něco je. Výsledný výpis datových položek jde po rozšiřujících se hladinách, do šířky (@ETIANM...). Prohledávání do hloubky (depth-first search) lze popsat stejně, jen se místo fronty užije zásobník, do něhož se vkládá i vybírá na jediném konci neboli vrcholu a který se popisuje termínem Last-In-First-Out (LIFO). Algoritmus lze realizovat ještě snadněji tím, že se práce se zásobníkem zakryje rekurzivním voláním: zpracuj data vstupního uzlu a prohledej levý a pravý podstrom, pokud existují. Výpis dat zkraje sleduje postupné zanořování do hloubky levých podstromů, odtud název algoritmu. Tři procedury pro průchod binárním stromem jsou v modulu odkázaném níže, třetí z nich i zde:

// Průchod binárním stromem do hloubky pomocí rekurze
procedure TraverseDepthFirst(const p : pBTree);
const
  variant = 1; // 1: preorder, 2: inorder, 2: postorder
begin
  if variant=1 then write(p^.A);
  if p^.pT<>nil then TraverseDepthFirst(p^.pT);
  if variant=2 then write(p^.A);
  if p^.pC<>nil then TraverseDepthFirst(p^.pC);
  if variant=3 then write(p^.A);
end;

Vidíme, že algoritmus má tři varianty lišící se pořadím výpisu dat. Např. varianta zvaná preorder vypisuje data z uzlů hned při prvním vstupu do uzlu (@EISH54V3UF...).

Po inicializaci a kontrolních výpisech tak konečně můžeme využít binární strom pro dekódování Morseovky. Čteme po řádcích zakódovanou zprávu, v řádcích lokalizujeme mezerou oddělované Morseovy kódy a dvěma mezerami oddělovaná slova. Každý Morseův kód analyzujeme (funkce DecodeChar) znak po znaku pohyby od kořenu stromu: tečka vede k posunu ve směru ukazatele pT, čárka ve směru pC; datová položka cílového uzlu pak obsahuje hledaný ASCII znak. Nesrozumitelné (nenalezené) kódy ignorujeme skokem z funkce (exit) s prázdnou návratovou hodnotou, skokem z cyklu (break) bychom reagovali vrácením znaku nalezeného poblíž.

// Dekódování Morseova kódu
function DecodeChar(const M : string) : string;
var
  p : pBTree;
  i : integer;
begin
  result:='';
  p:=root; // kořen stromu je deklarován vně funkce
  for i:=1 to Length(M) do begin
    case M[i] of
    '.': begin if p^.pT<>nil then p:=p^.pT else exit end;
    '-': begin if p^.pC<>nil then p:=p^.pC else exit end;
    else exit;
    end;
  end;
  result:=p^.A;
end;

Modul uMorseBTree.pas s popsanými procedurami je spolu s modulem uMorseList.pas, obsahujícím prostou dekódovací funkci sekvenčně procházející seznam Morseových kódů a příslušných ASCII znaků, testovacími programy, drivery MorseDecode a MorseEncode pro dekódování a kódování zpráv v Morseovce a programem MorsePlay pro přepípání zprávy je zde. (O tom, jak se vlastně ve Free Pascalu pípá, se rozepisujeme v následujícím odstavci.) V balíčku je i soubor morse1.in se zprávou v Morseovce, ten dekódujte, ať se něco dozvíte.

Můžete dekódovat i větší zprávu v 10megabytovém souboru morse2.in z doplňkového balíčku zde. Tato zpráva si žádá pozornější práci s řetězci. Generický typ string, který jsme doposud užívali, může být realizován typem ShortString, jehož délka je omezena na 255 znaků, nebo typem AnsiString s neomezenou délkou. Konkrétní typ se vybere podle nastavení překladače (direktivou {$H+} ve zdrojovém souboru, volbou -Sh na příkazovém řádku nebo nastavením Lazaru Projekt–Volby projektu–Analyzování–Použít Ansi řetězce), anebo jej lze v deklaraci typu uvést výslovně. My v naší zprávě řádky delší než 255 znaků máme, na jejich uložení do dlouhých řetězců nám tedy záleží a typ AnsiString v našem programu na patřičných místech použijeme. Jak se ukazuje, dekódování i takto dlouhé zprávy je sekundovou úlohou jak pomocí binárního stromu, tak prostým prohledáváním seznamu Morseových kódů, a jistou úsporu času pomocí binárního stromu si tak ani nestihneme vychutnat.

Temperované ladění

Odstavec zahájíme ukázkou, jak ve Free Pascalu pro Windows pípnout kvartu neboli “Hoří”:

// Pípnutí kvarty ve Windows
program testBeep;
uses Windows;
var
  f,d : integer;
begin
  f:=440;  // frekvence komorního a neboli a1 v Hz
  d:=1500; // doba v milisekundách
  beep(f,d);
  f:=round(f*4/3); // frekvenční interval kvarta
  beep(f,d);
  readln;
end.
_images/piano-550.png

12tónová chromatická řada na klaviatuře. Kresleno v gnuplotu (data a skript zde).

(Nedokončeno. Odkaz: Ladění chromatické řady PDF)

Cvičení č. 10

Obrátíme se k soustavám obyčejných diferenciálních rovnic s počáteční podmínkou. Použijeme Eulerovu metodu 1. řádu přesnosti a Rungovu-Kuttovu metodu 4. řádu přesnosti pro numerické řešení dvou systémů s dostupným analytickým řešením a soustavy plynoucí z Newtonova druhého pohybového zákona – spočteme balistickou křivku.

Rungovy-Kuttovy metody

Úlohou je nalézt numerickou aproximaci funkcí, pro které je předepsána soustava \(R\) obyčejných diferenciálních rovnic prvního řádu,

\[\frac{dy_r(x)}{dx}=f_r(x,y_1(x),\dots,y_R(x)),\ r=1,\dots,R,\]

neboli vektorově

\[\frac{d\mathbf{y}(x)}{dx}=\mathbf{f}(x,\mathbf{y}(x)),\]

doplněná počáteční podmínkou,

\[\mathbf{y}(x_0)=\mathbf{y}_0.\]

Numerickým řešením této tzv. počáteční úlohy na síti \(x_n,\ n=0,1,\dots,\) je pro každé \(r=1,\dots,R\) posloupnost \(y_{nr}\), která lépe či hůře aproximuje přesné řešení \(y_r(x_n)\).

Na přednášce je řeč o explicitní a implicitní Eulerově metodě 1. řádu, dvou Rungových-Kuttových metodách 2. řádu a jedné R.-K. metodě 4. řádu přesnosti. Zde si z těchto metod připravíme první a poslední. Pro explicitní Eulerovu metodu jsme si předvedli vzorec:

\[\mathbf{y}_{n+1}=\mathbf{y}_n+h\mathbf{f}(x_n,\mathbf{y}_n),\quad n=0,1,\dots,\]

vyjadřující posun z bodu \((x_n,\mathbf{y}_n)\) ve směru tečny k funkci \(y(x)\), jejíž směrnice je podle zadání rovna \(\mathbf{f}(x,\mathbf{y}_n)\). To odpovídá nultému a prvnímu členu Taylorova rozvoje \(y(x)\), proto mluvíme o 1. řádu přesnosti. Taylorův rozvoj do 4. řádu vystihuje Rungova-Kuttova metoda (RK4) v následujícím tvaru:

\[\begin{split}\mathbf{y}_{n+1}=\mathbf{y}_n+h\left(\frac16\mathbf{k}_1+\frac13\mathbf{k}_2+\frac13\mathbf{k}_3+\frac16\mathbf{k}_4\right),\quad n=0,1,\dots,\\ \mathbf{k}_1=\mathbf{f}(x_n,\mathbf{y}_n), \\ \mathbf{k}_2=\mathbf{f}(x_n+\frac12h,\mathbf{y}_n+\frac12h\mathbf{k}_1), \\ \mathbf{k}_3=\mathbf{f}(x_n+\frac12h,\mathbf{y}_n+\frac12h\mathbf{k}_2), \\ \mathbf{k}_4=\mathbf{f}(x_n+h,\mathbf{y}_n+h\mathbf{k}_3).\end{split}\]

Pro jeden krok z \(x_n\) do \(x_{n+1}\) vyčísluje metoda pravou stranu soustavy 4krát, krok metody lze interpretovat jako posun ve směru získaném jako vážený aritmetický průměr 4 tečen sestrojených na intervalu \(\langle x_n,x_{n+1}\rangle\).

Oběma metodami nejdříve vyřešíme následující rovnici s počáteční podmínkou:

\[dy(x)/dx=-y(x),\quad y(0)=1,\]

jejímž analytickým řešením je zjevně \(y(x)=e^{-x}\), a soustavu:

\[\begin{split}\begin{array}{ll} dy_1(x)/dx=y_2(x), & y_1(0)=0, \\ dy_2(x)/dx=-y_1(x),& y_2(0)=1, \end{array}\end{split}\]

s analytickým řešením \(y_1(x)=\sin x,\ y_2(x)=\cos x\).

Typické rozhraní řešičů poskytovaných numerickými knihovnami pro počáteční úlohy zahrnuje mezi vstupními argumenty meze nezávisle proměnné xmin a xmax, vektor Y0 s počáteční podmínkou v xmin a také odkaz na funkci, která vyčísluje vektor pravé strany pro argumenty x, Y. My k rozhraní naší procedury ještě přidáme počet nmax ekvidistantních kroků, které má metoda vykonat od xmin k xmax, a jako globální proměnnou zavedeme počet rovnic neq soustavy. Abychom mohli rovnice číslovat od 1 do neq, tlačí nás Pascal k použití statického pole (s dolní mezí nastavitelnou na 1), naopak pro uzly sítě se hodí indexování od 0 do nmax a my si tak můžeme procvičit i dynamická pole (s povinnou 0 v dolní mezi). Výstupní argumenty přinesou síť nezávisle proměnné ve vektoru xa (dynamické pole) a řešení Ya bude maticí (dynamickým vektorem statických vektorů), v první dimenzi indexovanou od 0 do nmax, v druhé od 1 do neq. Numerická řešení \(\mathbf{y_n}\) v jednotlivých uzlech sítě tak budou ležet na řádcích matice Ya.

Zde je program, který řeší naše dvě počáteční úlohy:

// Řešení soustav obyčejných diferenciálních rovnic s počáteční podmínkou
// dvě volitelné úlohy: exp(-x), sin(x)
// dvě volitelné metody: Eulerova, Rungova-Kuttova 4. řádu

program solveODE;

const
  neqmax = 5;    // max. počet rovnic (konstanta pro deklaraci pole)
var
  neq : integer; // aktuální počet rovnic (globální proměnná)

type
  tY = array [1..neqmax] of real; // vektor řešení
  tXa = array of real; // síť pro nezávisle proměnnou
  tYa = array of tY;   // hodnoty řešení na síti
  tF = function (const x : real; const Y : tY) : tY; // pravá strana soustavy

// rovnice s řešením y=exp(-x)
//   y'=-y
function oderhs1(const x : real; const Y : tY) : tY;
begin
  result[1]:=-Y[1];
end;

// soustava rovnic s řešením y=sin(x)
//   y1'=y2; y2'=-y1
function oderhs2(const x : real; const Y : tY) : tY;
begin
  result[1]:=Y[2];
  result[2]:=-Y[1];
end;

// Eulerova metoda 1. řádu
procedure Euler(oderhs : tF; const xmin,xmax : real; const Y0 : tY; const nmax : integer;
  var xa : tXa; var Ya : tYa);
var
  F : tY;
  h : real;
  n,i : integer;
begin
  if xa=nil then SetLength(xa,nmax+1);
  if Ya=nil then SetLength(Ya,nmax+1);
  xa[0]:=xmin;
  Ya[0]:=Y0; // for i:=1 to neq do Ya[0,i]:=Y0[i];
  h:=(xmax-xmin)/nmax;
  for n:=0 to nmax-1 do begin
    F:=oderhs(xa[n],Ya[n]);
    xa[n+1]:=xa[n]+h;
    for i:=1 to neq do Ya[n+1,i]:=Ya[n,i]+h*F[i];
  end;
end;

// Rungova-Kuttova metoda 4. řádu
procedure RK4(oderhs : tF; const xmin,xmax : real; const Y0 : tY; const nmax : integer;
  var xa : tXa; var Ya : tYa);
var
  k1,k2,k3,k4,Y : tY;
  h,x : real;
  n,i : integer;
begin
  if xa=nil then SetLength(xa,nmax+1);
  if Ya=nil then SetLength(Ya,nmax+1);
  xa[0]:=xmin;
  Ya[0]:=Y0; // for i:=1 to neq do Ya[0,i]:=Y0[i];
  h:=(xmax-xmin)/nmax;
  for n:=0 to nmax-1 do begin
    k1:=oderhs(xa[n],Ya[n]);
    xa[n+1]:=xa[n]+h;
    x:=xa[n]+h/2;
    for i:=1 to neq do Y[i]:=Ya[n,i]+h/2*k1[i];
    k2:=oderhs(x,Y);
    for i:=1 to neq do Y[i]:=Ya[n,i]+h/2*k2[i];
    k3:=oderhs(x,Y);
    x:=xa[n+1];
    for i:=1 to neq do Y[i]:=Ya[n,i]+h*k3[i];
    k4:=oderhs(x,Y);
    for i:=1 to neq do Ya[n+1,i]:=Ya[n,i]+h*(k1[i]/6+k2[i]/3+k3[i]/3+k4[i]/6);
  end;
end;

var
  Xa : tXa;
  Ya : tYa;
  Y0 : tY;
  xmin,xmax : real;
  uloha,metoda,nmax,n,i : integer;
  oderhs: tF;
  s : text;
  soubor : string;

begin
  uloha :=2; // typ úlohy:  1 exp(-x), 2 sin(x)
  metoda:=2; // typ metody: 1 Euler, 2 RK4
  soubor:='out2RK';

  case uloha of
  1: begin
       neq:=1;
       oderhs:=@oderhs1;
       xmin:=0;   // meze nezávisle proměnné
       xmax:=2;
       Y0[1]:=1;  // počáteční podmínka
       nmax:=20;  // počet vyčíslení pravé strany
     end;
  2: begin
       neq:=2;
       oderhs:=@oderhs2;
       xmin:=0;
       xmax:=pi*2;
       Y0[1]:=0;
       Y0[2]:=1;
       nmax:=100;
     end;
  end;

  case metoda of
  1: Euler(oderhs,xmin,xmax,Y0,nmax,xa,Ya);
  2: begin
       nmax:=nmax div 4;
       RK4(oderhs,xmin,xmax,Y0,nmax,xa,Ya);
     end;
  end;

  assign(s,soubor);
  rewrite(s);
  for n:=0 to nmax do begin
    write(s,xa[n]:12:5);
    for i:=1 to neq do write(s,' ',Ya[n,i]:12:5);
    writeln(s);
  end;
  close(s);
  writeln('Soubor ',soubor,' ulozen.');
  readln;
end.

Pro exponenciálu na intervalu <0,2> jsme volili 20 kroků Eulerovy metody a 5 kroků RK4 metody, aby počet vyčíslení pravé strany (míra časové náročnosti těchto metod) byl v obou případech totožný. Na obrázku porovnávajícím analytické řešení s oběma numerickými pak vidíme, jak výpočet metodou RK4 lépe přiléhá k analytickému řešení:

_images/ode-exp.png

To platí i pro numerický výpočet goniometrické funkce – všimněme si, že Eulerově metodě zde nestačí pro jednu periodu ani 100 kroků, aby uspokojila oko.

_images/ode-sin.png

Použité skripty pro gnuplot:

set terminal pngc size 800,600
set output 'ode-exp.png'
plot [-.05:2.05][0:1.05] exp(-x) lw 4,\
     'out1E' u 1:2:(.03) with circles lw 2 title 'Euler',\
     'out1RK' u 1:2:(.05) with circles lw 2 title 'RK4'
set output
unset terminal
set terminal pngc size 800,600
set output 'ode-sin.png'
plot [-.25:pi*2+.25][-1.25:1.25] sin(x) lw 4,\
     'out2E' u 1:2:(.06) with circles lw 2 title 'Euler',\
     'out2RK' u 1:2:(.10) with circles lw 2 title 'RK4'
set output
unset terminal

Balistické křivky

Výše připravený program je schopen větších věcí než jen reprodukovat analytické funkce. Vezměme si Newtonův druhý pohybový zákon,

\[Md^2\mathbf{x}(t)/dt^2=\mathbf{F}(t,\mathbf{x}(t)),\]

popisující pohyb \(\mathbf{x}(t)\) tělesa o hmotnosti \(M\) vlivem gravitační síly (při Zemi konstantní a orientované dolů) a odporové síly (působící proti směru pohybu a řekněme úměrné druhé mocnině rychlosti tělesa), \(\mathbf{F}=\mathbf{F}_g+\mathbf{F}_o=-Mg\,\mathbf{e}_z-K|\mathbf{v}|\mathbf{v}\). Víme, že v takovém prostředí se těleso pohybuje po balistické křivce a že analytické řešení tato úloha nemá. Pro nás ideální.

Numerické metody jsme zapsali pro soustavy s prvními derivacemi. Druhý Newtonův zákon obsahuje druhou derivaci posunutí, což však explicitním zahrnutím složek rychlosti snadno upravíme na (větší) soustavu s prvními derivacemi. Problém stačí popsat ve dvou dimenzích, tedy \(\mathbf{x}(t)=(x(t),z(t))\) a \(\mathbf{v}(t)=(v_x(t),v_z(t))\):

\[\begin{split}\frac{d}{dt}\left(\begin{array}{r} x(t) \\ z(t) \\ v_x(t) \\ v_z(t) \\ \end{array}\right) = \left(\begin{array}{l} v_x(t) \\ v_z(t) \\ F_x(t,x,z,v_x,v_z)/M \\ F_z(t,x,z,v_x,v_z)/M \\ \end{array}\right) = \left(\begin{array}{l} v_x(t) \\ v_z(t) \\ -K/M\,v(t)\,v_x(t) \\ -g-K/M\,v(t)\,v_z(t) \\ \end{array}\right)\,.\end{split}\]

Našemu programu doplníme funkci pro pravou stranu soustavy, parametrizovanou pro atletickou kouli a pingpongový míček:

// soustava rovnic pro šikmý vrh v gravitačním poli s odporovou silou
function SikmyVrh(const x : real; const Y : tY) : tY;
const
// atletická koule
// M=7.26;  // hmotnost [kg]
// g=9.81;  // tíhové zrychlení [m/s^2]
// K=0.5*0.48*pi*sqr(0.065)*1.2; // koeficient odporu prostředí [kg/m]
// pingpongový míček
   M=0.0027;// hmotnost [kg]
   g=9.81;  // tíhové zrychlení [m/s^2]
   K=0.5*0.48*pi*sqr(0.020)*1.2; // koeficient odporu prostředí [kg/m]
var
  v : real;
begin
  v:=sqrt(sqr(Y[3])+sqr(Y[4]));
  result[1]:=Y[3];
  result[2]:=Y[4];
  result[3]:=-K/M*v*Y[3];
  result[4]:=-g-K/M*v*Y[4];
end;

Použili jsme vzorec \(K=\frac12C_DS\rho\) (viz Wikipedie), kde \(S\) je průřez tělesa, \(\rho\) hustota vzduchu a \(C_D\) součinitel odporu, pro který jsme z tabulky zde převzali hodnotu 0,48 (stejnou pro kouli i Volkswagen Beetle). Do úseku s přepínačem úloh přidáme požadavek na řešení 2,2sekundového letu koule či míčku mrštěného z výšky 2 m pod úhlem 42 stupňů rychlostí 14 m/s:

3: begin
     neq:=4;
     oderhs:=@SikmyVrh;
     xmin:=0;   // meze nezávisle proměnné
     xmax:=2.2;
     Y0[1]:=0;  // počáteční podmínka
     Y0[2]:=2;
     Y0[3]:=14*cos(pi*42/180);
     Y0[4]:=14*sin(pi*42/180);
     nmax:=1000;// počet vyčíslení pravé strany
   end;

Odpovědí je obrázek, získaný následujícím skriptem (celý program s vizualičními skripty je sbalen zde):

set terminal pngc size 800,600
set output 'ode-ballistic.png'
set xlabel 'x [m]'
set ylabel 'z [m]'
set grid
fac=.25
plot 'out3RK' using 2:3 with lines lw 5 lc 'red' title 'z(x)',\
  '' using 2:3:($4*fac):($5*fac) every 20 with vectors head empty lw 2 lc 'green' title 'v(x)'
set output
unset terminal
_images/ode-ballistic.png

Vidíme, že těleso se po chvilce letu vzhůru vrací k zemi, ovšem po nesymetrické dráze (při nenulovém K), a tušíme, že jeho rychlost se blíží konstantě (díky vyrovnávání sil \(\mathbf{F}_g\) a \(-\mathbf{F}_o\)). Míček jsme vrhli do výšky necelých 4,5 m a dovrhli do vzdálenosti přes 8 m (pro \(z=0\)). Zajímavým tématem je přesnější určení obou údajů, tedy nalezení nulového bodu složky rychlosti \(v_z(t_1)\) a zjištění \(z(t_1)\) a nalezení nulového bodu složky polohy \(z(t_2)\) a zjištění \(x(t_2)\). Hledač kořenů máme zpracován z dřívějška, ten ovšem hledá nulové body funkce spojité proměnné a my zde máme diskrétní body. Bylo by tedy ještě třeba zapracovat vhodnou interpolaci pro složky řešení Ya.

Pokud bychom byli schopni vrhnout atletickou kouli toutéž rychlostí a pod týmž úhlem jako pingpongový míček, měli bychom zanechat programování a věnovat se více vrhání, neboť jsme jen 22 cm pod českým rekordem.

Poznamenejme ještě, že nahradíme-li homogenní gravitační pole obecným gravitačním zákonem (odporovou sílu lze vynechat), můžeme počítat dráhy těles při jejich letu vesmírem. Mohli bychom např. numericky ověřit hodnoty kosmických rychlostí uvedené zde. Poznamenejme také, že úlohy tohoto typu (Hamiltonovské systémy) je lépe řešit pro ně vyvinutými metodami zvanými symplektické integrátory, ovšem s dostatečně malým krokováním zvládne s přijatelnou přesností několik oběhů i Eulerova metoda.

Zbývá zadat domácí úkol. Tím je spočítat a vizualizovat Lorenzův atraktor neboli trajektorii bodu (v 3D stavovém prostoru) popsanou třemi obyčejnými diferenciálními rovnicemi:

\[\begin{split}\frac{d}{dt}\left(\begin{array}{r} x \\ y \\ z \\ \end{array}\right) = \left(\begin{array}{l} \sigma(y-x) \\ \rho x-y-xz \\ xy-\beta z \\ \end{array}\right), \ \sigma=10,\ \beta=8/3,\ \rho=28.\end{split}\]

Počáteční podmínku je třeba zadat, ale může být libovolná (nenulová). Hledaná trajektorie je atraktorem, takže si k ní bod najde cestu odkudkoliv už sám. Je třeba vhodně volit interval integrace a počet kroků, aby se výsledný obrázek podobal obrázku na Wikipedii nebo třeba tady. To je ve vašich rukou. Nezapomeňte, že 3D Lissajousovy obrazce jsme kreslili v ParaView – to by i zde vytvořilo silný dojem, silnější než gnuplot.

Několik pojmenovaných podivných atraktorů zobrazených v ParaView je shrnuto v tomto PDF a jejich kombinovanou ukázkou tuto část uzavřeme:

_images/ode-attractors1.png _images/ode-attractors2.png _images/ode-attractors3.png

Cvičení č. 11

Poslední cvičení je věnováno přípravě a vizualizaci dat ve třech dimenzích. Půjde o banán, srdce, sférické harmoniky a atomové orbitaly. Pro kreslení užijeme gnuplot, Octave a ParaView.

Parametrické funkce

Nejprve se budeme zabývat situací, kdy je objekt (funkce) ve třech dimenzích (3D) popsán parametricky. Máme na mysli zápis

\[x = x(u,v,w),\ y = y(u,v,w),\ z = z(u,v,w),\]

kde \(u,v,w\) jsou nezávislé parametry a \(x,y,z\) kartézské souřadnice bodů náležejících k 3D objektu. Při dvou nezávislých parametrech půjde o 2D plochu, při jediném o 1D křivku. V PDF souboru Plochy – potrava pro oko jsou shrnuty parametrické funkce popisující mj. banán, vajíčko, bagel, šneka či kapříka. Sem převezmeme banán:

\[x=a\sin\theta\cos\phi,\ y=b\sin\theta\sin\phi,\ z=b(\cos\theta+c\sin^2\theta\cos^2\phi),\]

kde \(0\le\theta\le\pi\) a \(0\le\phi\le2\pi\) jsou nezávislé parametry a ostatní fixujeme na hodnotách a=1, b=1/5 a c=2.

Pomocí vnořených cyklů systematicky projdeme požadované intervaly nezávislých parametrů, vytvoříme tak síť uzlových bodů a zapíšeme trojice jejich kartézských souřadnic. K uzlům bychom mohli přidávat další (skalární nebo vektorová) data a ty pak – kromě samotného tvaru objektu – také vizualizovat (např. barevně nebo šipkami). Pokud by se tato data dala vyjádřit jen ze souřadnic, mohli bychom jejich výpočet ponechat až na vizualizační program. My budeme chtít obarvit banán realisticky na žluto a na obou koncích do černa. Uhodneme, že pro konce platí \(\sqrt{x^2+y^2+z^2}\gtrsim a\), a ke každému uzlu přidáme velikost jeho průvodiče, kterou navážeme na vhodnou barevnou paletu, takovou, aby pro hodnoty menší než a byla žlutá a pro hodnoty větší než a přešla do černé. Pro gnuplot budeme generovat datový soubor jen se souřadnicemi, velikost průvodiče dopočteme v gnuplotu. Pro ParaView zapíšeme VTK soubor (ve variantě křivočará síť) s hlavičkou, sekcí souřadnic a sekcí dat. Komentované ukázky, jak mají VTK soubory vypadat, obsahuje PDF soubor Poprvé s ParaView, a vykoukat to snad lze i z následujících zdrojových textů.

// Příprava dat pro vizualizaci banánu v gnuplotu a ParaView

program Banan;
const
  iOut=2;              // vystupni format: 1 gnuplot, 2 ParaView
  a=1; b=0.2; c=2;     // parametry bananu
  nr=0;                // 0 pro 1 plochu
  nth=90;              // nth+1: pocet bodu sfericke site ve smeru theta
  nph=180;             // nph+1: pocet bodu sfericke site ve smeru phi
  thmin=0; thmax=pi;   // theta_min a theta_max
  phmin=0; phmax=pi*2; // phi_min a phi_max

var
  ir,ith,iph : integer;
  th,ph,x,y,z : real;
  s : array [0..nph,0..nth,0..nr] of real; // skalarni data v uzlech site
  ofile : text;
  soubor : string;
begin
  case iOut of
  1: soubor:='banan.dat'; // pro gnuplot
  2: soubor:='banan.vtk'; // pro ParaView
  end;
  assign(ofile,soubor);
  rewrite(ofile);

  case iOut of
  2: begin // hlavicka datoveho souboru VTK s krivocarymi souradnicemi
    writeln(ofile,'# vtk DataFile Version 3.0');
    writeln(ofile,'vtk output');
    writeln(ofile,'ASCII');
    writeln(ofile,'DATASET STRUCTURED_GRID');
    writeln(ofile,'DIMENSIONS ',nph+1,' ',nth+1,' ',nr+1);
    writeln(ofile,'POINTS ',(nph+1)*(nth+1)*(nr+1),' float');
    end;
  end;

  ir:=nr;
  for ith:=0 to nth do begin // sekce souradnic
    th:=thmin+ith*(thmax-thmin)/nth;
    for iph:=0 to nph do begin
      ph:=phmin+iph*(phmax-phmin)/nph;
      x:=a*sin(th)*cos(ph);
      y:=b*sin(th)*sin(ph);
      z:=b*(cos(th)+c*sqr(sin(th)*cos(ph)));
      s[iph,ith,ir]:=sqrt(x*x+y*y+z*z);
      writeln(ofile,x:8:4,y:8:4,z:8:4);
    end;
  end;

  case iOut of
  2: begin // sekce dat
    writeln(ofile,'POINT_DATA ',(nph+1)*(nth+1)*(nr+1));
    writeln(ofile,'SCALARS MyData float 1');
    writeln(ofile,'LOOKUP_TABLE default');
    ir:=nr;
    for ith:=0 to nth do begin
      for iph:=0 to nph do begin
        writeln(ofile,s[iph,ith,ir]:8:4);
      end;
    end;
    end;
  end;

  close(ofile);
  writeln('Soubor ',soubor,' ulozen.');
  readln;
end.

Datový soubor banan.dat pro gnuplot obsahuje pouze tři sloupce se souřadnicemi síťových bodů. Velikost průvodiče v bodech sítě vytvoříme až v gnuplotu příkazem splot jako virtuální čtvrtý sloupec, klauzulí palette zohledněný pro obarvení jednotlivých uzlů. Žlutočernou barevnou paletu si připravíme příkazem set palette defined.

# gnuplot: banán
unset border
unset xtics
unset ytics
unset ztics
unset colorbox
unset key
set palette defined (0 'yellow', 0.8 'yellow', 1.1 'black')
set view 75,135
splot [-1:1][-1:1][-1:1] 'banan.dat' using 1:2:3:(sqrt($1**2+$2**2+$3**2)) palette with lines
set terminal pngc size 2500,1500 crop
set output 'banan-gnuplot.png'
replot
set output
unset terminal

Banán v gnuplotu...

_images/banan-gnuplot.png

S přichystaným souborem banan.vtk lze v ParaView už jen klikat myší: File–Open, načíst soubor banan.vtk, Apply, pootočit banánem při stisknutém levém tlačítku myši (a zkusit, co dělá pohyb myši při stisku ostatních tlačítek), nahoře ikonkami Toggle Color Legend Visibility a Show Orientation Axes vypnout barevnou škálu a osový kříž, ikonkou Load a color palette nastavit White Background, vlevo v části Coloring stisknout Edit pro Color Map Editor, v obrázku Mapping Data dvojkliknout na levém kolečku ve vodorovném pruhu barevné škály a změnit barvu na žlutou, střední kolečko posunout někam pod 1,0 a změnit také na žlutou, pravému kolečku nastavit černou, pak ještě vlevo v Lighting pro parádu přidat na Specular a konečně ve File–Save Screenshot nastavit velikost výstupního PNG obrázku a uložit ho.

Banán v ParaView...

_images/banan-paraview.png

Vše související s těmito dvěma banány je sbaleno zde.

Implicitní funkce

Dále se zaměříme na případ, kdy je plocha popsána implicitně rovnicí

\[f(x,y,z)=0.\]

Ze souboru Plochy – potrava pro oko (a odsud) nyní vybereme rovnici popisující 3D srdce:

\[\left(ax^2+by^2+cz^2-1\right)^3-dx^2z^3-ey^2z^3=0,\]

kde sada konstant a=1, b=9/4, c=1, d=1, e=9/80 definuje spíše široké srdce, zatímco sada a=2, b=2, c=1, d=1, e=1/10 popisuje úzký exemplář.

Běžnou cestou k vizualizaci implicitně zadané plochy je tabelování funkce \(f(x,y,z)\) na jednoduché 3D oblasti (např. na kvádru) a zobrazení její (nulové) izoplochy. Tuto schopnost nemá gnuplot, má ji však Octave (funkce isosurface) i ParaView (modul Contour).

Skript pro Octave, který na kartézské síti [x,y,z] funkci \(f(x,y,z)\) tabeluje do pole s a rovnou i kreslí její izoplochu, isosurface(x,y,z,s,0), následuje. Je použitelný i v MATLABu, který – na rozdíl od dnešního Octavu – umí zacházet i s odlesky, což nahrazením aktivního řádku patch řádkem zakomentovaným také vyzkoušíme.

% Octave: 3D srdce
a=1; b=9/4; c=1; d=1; e=9/80;
xmin=-1.2; xmax=1.2; ymin=-.75; ymax=.75; zmin=-1.2; zmax=1.5;
nx=100; ny=100; nz=300;
xvec=linspace(xmin,xmax,nx+1); yvec=linspace(ymin,ymax,ny+1); zvec=linspace(zmin,zmax,nz+1);
[x,y,z]=ndgrid(xvec,yvec,zvec);
s=(a*x.^2+b*y.^2+c*z.^2-1).^3-d*x.^2.*z.^3-e*y.^2.*z.^3;
set(gcf,'Color','white'); clf; axis off; view([-1 5 0.5]); % bílé pozadí, bez os, pozorovací bod
[f,v]=isosurface(x,y,z,s,0);
patch('Faces',f,'Vertices',v,'FaceColor','red','EdgeColor',[.5 .0 .0]);
% MATLAB (with lighting)
% patch('Faces',f,'Vertices',v,'FaceColor','red','EdgeColor','none'); camlight('right'); lighting gouraud;

Srdce v Octavu (3D dojem vytvářejí hrany elementů tvořících izoplochu) a MATLABu (3D dojem díky stínování)...

_images/srdce-octave.png _images/srdce-matlab.png

Datový soubor pro ParaView s funkcí \(f(x,y,z)\) tabelovanou ve VTK formátu (varianta: rovnoměrná síť podle souboru Poprvé s ParaView) připravíme pomocí následujícího programu. Tabelujeme pole s na síti x,y,z, stejně jako v předešlém skriptu. Přípravu izoplochy necháváme na ParaView.

// Příprava dat pro vizualizaci srdce v ParaView

program Srdce;
const
  a=1; b=9/4; c=1; d=1; e=9/80; // parametry sirokeho srdce
//a=2; b=2;   c=1; d=1; e=1/10; // parametry uzkeho srdce
  xmin=-1.2; xmax=1.2;  // meze ve smeru x
  ymin=-0.75;ymax=0.75; // meze ve smeru y
  zmin=-1.2; zmax=1.5;  // meze ve smeru z
  nx=100;               // nx+1: pocet bodu kartezske site ve smeru x
  ny=100;               // pocty bodu ve smeru y
  nz=100;               // pocty bodu ve smeru z

var
  ix,iy,iz : integer;
  x,y,z,rr3,z3 : real;
  s : array [0..nx,0..ny,0..nz] of real; // skalarni data v uzlech site
  ofile : text;
  soubor : string;
begin
  soubor:='srdce.vtk'; // pro ParaView
  assign(ofile,soubor);
  rewrite(ofile);

  // hlavicka datoveho souboru VTK s rovnomernou kartezskou siti
  writeln(ofile,'# vtk DataFile Version 3.0');
  writeln(ofile,'vtk output');
  writeln(ofile,'ASCII');
  writeln(ofile,'DATASET STRUCTURED_POINTS');
  writeln(ofile,'DIMENSIONS ',nx+1,' ',ny+1,' ',nz+1);
  writeln(ofile,'ORIGIN  ',xmin,' ',ymin,' ',zmin);
  writeln(ofile,'SPACING ',(xmax-xmin)/nx,' ',(ymax-ymin)/ny,' ',(zmax-zmin)/nz);

  for ix:=0 to nx do begin // sekce souradnic
    x:=xmin+ix*(xmax-xmin)/nx;
    for iy:=0 to ny do begin
      y:=ymin+iy*(ymax-ymin)/ny;
      for iz:=0 to nz do begin
        z:=zmin+iz*(zmax-zmin)/nz;
        rr3:=a*sqr(x)+b*sqr(y)+c*sqr(z)-1;
        rr3:=rr3*rr3*rr3;
        z3:=z*z*z;
        s[ix,iy,iz]:=rr3-d*sqr(x)*z3-e*sqr(y)*z3;
      end;
    end;
  end;

  writeln(ofile,'POINT_DATA ',(nx+1)*(ny+1)*(nz+1)); // sekce dat
  writeln(ofile,'SCALARS MyScalar float 1');
  writeln(ofile,'LOOKUP_TABLE default');
  for iz:=0 to nz do begin
    for iy:=0 to ny do begin
      for ix:=0 to nx do begin // prvni souradnice se meni nejrychleji
        writeln(ofile,s[ix,iy,iz]:15);
      end;
    end;
  end;

  close(ofile);
  writeln('Soubor ',soubor,' ulozen.');
  readln;
end.

Izoplochy v ParaView: File–Open, načíst soubor srdce.vtk, Apply, vlevo nahoře kliknout na ikonu Contour (nebo ekvivalentně volit Filters–Common–Contour), vlevo v Isosurfaces změnit hodnotu na 0, Apply, nahoře ikonou Set view direction to +Y nebo myší nastavit vhodnější úhel pohledu, vlevo v Coloring–Edit vybrat červenou barvu, odstranit bounding box kliknutím na ikoně oka vlevo u souboru srdce.vtk a stejně jako u banánu vybrat barvu pozadí, nastavit odrazivost a uložit obrázek.

Srdce v ParaView...

_images/srdce-paraview-1.png _images/srdce-paraview-2.png

Můžeme ovšem ještě pokračovat: View–Animation View, v záhlaví okna nastavit hodnoty Start Time a End Time třeba na -10 a 10, No. Frames pak řekněme 101, dále modré plus u Contour1–Isosurfaces, tam pak dvojklik pro Animation Keyframes, nastavit shodné meze pro čas (-10, 10) a hodnoty Value od -.05 do +.05, OK a zelená šipka pro Play. Po File–Save Animation získáme toto čili animaci izoploch blízkých nulové izoploše.

Program, skript, VTK soubor (binární, získaný jinak než přiloženým programem) a obrázky jsou sbaleny zde.

Sférické harmonické funkce

Se sférickými harmonikami \(Y_l^m(\theta,\phi)\) se potkáváme v situacích, kde hraje roli sférická symetrie (mj. v jaderné fyzice, kvantové fyzice, astronomii nebo geofyzice). Sférické harmoniky jsou ortonormální na povrchu koule, funkce definované na kouli lze reprezentovat jako řadu sférických harmonik, funkce \(r^lY_l^m(\theta,\phi)\) a \(r^{-l-1}Y_l^m(\theta,\phi)\) vyhovují Laplaceově rovnici v 3D, \(\Delta f(r,\theta,\phi)=0\). Ať tak či onak, jejich matematická definice zní:

\[Y_l^m(\theta,\phi)=N_l^m P_l^m(\cos\theta)e^{im\phi},\]

kde \(\theta,\ \phi\) jsou úhlové sférické souřadnice, normovací konstanty \(N_l^m=(-1)^m\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}\) a \(P_l^m(\cos\theta)\) značí přidružené Legendrovy funkce. Ty jsou zavedeny (a ortogonální) na intervalu \(\langle-1,1\rangle\) a numericky se vyčíslují pomocí rekurentního vztahu:

\[P_{k+2}^m(x)=\frac{(2k+3)xP_{k+1}^m(x)-(k+m+1)P_k^m(x)}{k-m+2},\ k=m,m+1,\dots,\]

který se startuje pomocí vzorců \(P_m^m(x)=(-1)^m\frac{(2m)!}{2^m m!}(1-x^2)^{\frac m2}\) a \(P_{m+1}^m(x)=(2m+1)xP_m^m(x)\). Tyto tři vztahy zobecňují rekurentní vzorec pro Legendrovy polynomy \(P_l(x)\equiv P_l^0(x)\), s nimiž jsme se setkali na přednášce o numerickém integrování: \(P_{k+2}(x)=\frac{(2k+3)xP_{k+1}(x)-(k+1)P_k(x)}{k+2},\ k=0,1,\dots,\) se startovacími hodnotami \(P_0(x)=1,\ P_1(x)=x\).

// Vyčíslení přidružené Legendrovy funkce P_l^m(x)
function Leg(const l,m : integer; const x : real) : real;
var
  k : integer;
  P1,P0 : real;
begin
  P0:=1;                   // P_0^0
  if m>0 then begin        // P_m^m
    P0:=power((1-x)*(1+x),m*0.5);
    for k:=1 to m do P0:=-P0*(2*k-1);
  end;
  result:=P0;
  P1:=P0*x*(2*m+1);        // P_{m+1}^m
  if l=m+1 then result:=P1;
  for k:=m to l-2 do begin // P_l^m
    result:=(P1*x*(2*k+3)-P0*(k+m+1))/(k-m+2);
    P0:=P1;
    P1:=result;
  end;
end;

// Vyčíslení reálné sférické harmonické funkce Y_l^m(theta,phi)
function reYlm(const l,m : integer; const theta,phi : real) : real;
var
  Nlm,x : real;
  k : integer;
begin
  Nlm:=(2*l+1)/(4*pi); // normovaci konstanta
  for k:=l-m+1 to l+m do Nlm:=Nlm/k;
  Nlm:=sqrt(Nlm); if odd(m) then Nlm:=-Nlm;
  x:=cos(theta);       // argument pridruzene Legendrovy funkce
  result:=Nlm*Leg(l,m,x)*cos(m*phi);
end;

Sférické harmoniky můžeme vizualizovat poctivě na povrchu koule vynesením barvy přidružené hodnotě jejich reálné části, tedy barevnou reprezentací skalárního pole \(\operatorname{Re}(Y_l^m(\theta,\phi))\) na ploše popsané parametricky:

\[x = \sin\theta\cos\phi,\ y = \sin\theta\sin\phi,\ z = \cos\theta.\]

Zonální sférické harmoniky \(Y_l^0(\theta,\phi),\ l=0,1,2,3,\) na kulové ploše...

\[\textstyle Y_0^0=N_0^0,\ Y_1^0=N_1^0\cos\theta,\ Y_2^0=N_2^0\left(\frac32\cos^2\theta-\frac12\right),\ Y_3^0=N_3^0\left(\frac52\cos^3\theta-\frac32\cos\theta\right)\]
_images/y00-bw.png _images/y10-bw.png _images/y20-bw.png _images/y30-bw.png

Mnozí to však považují za nudné a obarvenou kulovou plochu, \(r=1,\) ještě přetvarují podle předpisu \(r=|\operatorname{Re}(Y_l^m(\theta,\phi))|\) nebo \(r=|\operatorname{Re}(Y_l^m(\theta,\phi))|^2\):

\[x = |\operatorname{Re}(Y_l^m(\theta,\phi))|\sin\theta\cos\phi,\ y = |\operatorname{Re}(Y_l^m(\theta,\phi))|\sin\theta\sin\phi,\ z = |\operatorname{Re}(Y_l^m(\theta,\phi))|\cos\theta.\]

Zonální sférické harmoniky \(Y_l^0(\theta,\phi),\ l=0,1,2,3,\) v ParaView (osa z míří vzhůru mírně napravo)...

\[\textstyle Y_0^0=N_0^0,\ Y_1^0=N_1^0\cos\theta,\ Y_2^0=N_2^0\left(\frac32\cos^2\theta-\frac12\right),\ Y_3^0=N_3^0\left(\frac52\cos^3\theta-\frac32\cos\theta\right)\]
_images/y00.png _images/y10.png _images/y20.png _images/y30.png

Reálné sektorové sférické harmoniky \(\operatorname{Re}(Y_l^l(\theta,\phi)),\ l=0,1,2,3\) (osa z míří doprava mírně vpřed)...

\[\textstyle Y_0^0=N_0^0,\ Y_1^1=N_1^1(-\sin\theta)e^{i\phi},\ Y_2^2=N_2^2(3\sin^2\theta)e^{2i\phi},\ Y_3^3=N_3^3(-15\sin^3\theta)e^{3i\phi}\]
_images/y00.png _images/y11.png _images/y22.png _images/y33.png

...a jejich kombinace – proč ji nenazvat kapřík:

_images/yll-kaprik1.png _images/yll-kaprik2.png _images/yll-kaprik3.png

Příprava dat i práce s ParaView je tu shodná jako při malování banánu neboli parametricky zadané funkce (o dva odstavce výše). Program, který tabeluje předepsanou sférickou harmonickou funkci na předepsané síti ve formátu pro gnuplot a ParaView, je sbalen zde. Jím zhotovená data malovaná v ParaView můžeme porovnávat s obrázkem na Wikipedii.

Atomové orbitaly

Atomový orbital znázorňuje pravděpodobnost výskytu elektronu v okolí atomového jádra. Ta se odvozuje od 3D vlnové funkce \(\psi_{nlm}(r,\theta,\phi)\), která má pro vodíkový atom následující tvar (odvození např. zde):

\[\psi_{nlm}(r,\theta,\phi)=\sqrt{\left(\frac{2}{na}\right)^3 \frac{(n-l-1)!}{2n(n+l)!}} \ e^{-\frac{r}{na}} \ \left(\frac{2r}{na}\right)^l \ L_{n-l-1}^{2l+1}\left(\frac{2r}{na}\right) \ Y_l^m(\theta,\phi).\]

Indexy \(n,l,m\) se nazývají hlavní, vedlejší a magnetické kvantové číslo, konstanta \(a\) = 5,292.10-11 m je Bohrův poloměr, \(L_k^\alpha(x)\) značí zobecněné Laguerrovy polynomy a \(Y_l^m(\theta,\phi)\) sférické harmonické funkce, diskutované o odstavec výše. Zobecněné Laguerrovy polynomy, ortogonální s vahou \(x^\alpha e^{-x}\) na intervalu \(\langle0,\infty)\), vyčíslíme pomocí rekurentního vztahu:

\[L_{k+2}^\alpha(x)=\frac{(2k+3+\alpha-x)L_{k+1}^\alpha(x)-(k+1+\alpha)L_k^\alpha(x)}{k+2},\ k=0,1,\dots,\]

startovaného z hodnot \(L_0^\alpha(x)=1\) a \(L_1^\alpha(x)=1+\alpha-x\).

// Vyčíslení Laguerrova polynomu L_l(x)
function Lag(const l : integer; const x : real) : real;
var
  L0,L1 : real;
  k : integer;
begin
  L0:=1;
  L1:=1-x;
  result:=L0;
  if l=1 then result:=L1;
  for k:=0 to l-2 do begin
    result:=(L1*(2*k+3-x)-L0*(k+1))/(k+2);
    L0:=L1;
    L1:=result;
  end;
end;

// Vyčíslení zobecněného Laguerrova polynomu L_l^alpha(x)
function Lag(const l : integer; const alpha : real; const x : real) : real;
var
  L0,L1 : real;
  k : integer;
begin
  L0:=1;                   // L_0^alpha
  L1:=1+alpha-x;           // L_1^alpha
  result:=L0;
  if l=1 then result:=L1;
  for k:=0 to l-2 do begin // L_l^alpha
    result:=(L1*(2*k+3+alpha-x)-L0*(k+1+alpha))/(k+2);
    L0:=L1;
    L1:=result;
  end;
end;

// Vyčíslení reálné části vlnové funkce psi_{nlm}(r,theta,phi)
function psi(const n,l,m : integer; const r,theta,phi : real) : real;
const
  a=5.291772e-11;        // Bohrův poloměr
var
  Nnl,x : real;
  k : integer;
begin
  Nnl:=power(2/(n*a),3); // normovaci konstanta
  for k:=n-l to n+l do Nnl:=Nnl/k;
  Nnl:=sqrt(Nnl/(2*n));
  x:=2*r/(n*a);          // argument Laguerrova polynomu
  result:=Nnl*exp(-x*0.5)*power(x,l)*Lag(n-l-1,2*l+1,x)*reYlm(l,m,theta,phi);
end;

Pro zobrazení atomového orbitalu pro kvantová čísla \(n,l,m\) je třeba vyčíslit na vhodně zvolené oblasti vlnovou funkci \(\psi_{nlm}(r,\theta,\phi)\), přejít k její reálné reprezentaci (vhodnou lineární kombinací komplexních vlnových funkcí nebo jednoduše omezením se na \(\operatorname{Re}(\psi_{nlm})\)), spočítat pravděpodobnost \(P(r,\theta,\phi)\) jako kvadrát reálné vlnové funkce a vybrat její izoplochu, \(P(r,\theta,\phi)=\mathrm{const},\) která ohraničuje oblast s dostatečně vysokou pravděpodobností výskytu elektronu nebo prostě jen dobře vypadá. (Vytvořit izoplochu v ParaView jsme si už zkusili na příkladu 3D srdce.) Na Wikipedii vizualizují orbitaly pomocí samotné vlnové funkce (nikoliv pravděpodobnosti), tak uděláme totéž, abychom mohli srovnávat.

Orbitaly vodíkového atomu: izoplochy \(\psi_{210},\ \psi_{310},\ \psi_{320},\ \psi_{321}\) v ParaView...

_images/orb-210.png _images/orb-310.png _images/orb-320.png _images/orb-321.png

Orbitaly \(\psi_{410},\ \psi_{420},\ \psi_{421},\ \psi_{432}\)...

_images/orb-410.png _images/orb-420.png _images/orb-421.png _images/orb-432.png

Orbitaly \(\psi_{510},\ \psi_{520},\ \psi_{530},\ \psi_{540}\)...

_images/orb-510.png _images/orb-520.png _images/orb-530.png _images/orb-540.png

Program pro přípravu atomových orbitalů sbalený zde se podobá programu pro přípravu srdce neboli izoplochy implicitní funkce (o dva odstavce výše) a čerpá funkce z programu pro přípravu sférických harmonik (o odstavec výše).

Octave

Programovacím jazykem našeho předmětu je Pascal, mohlo by to být i C, C++ nebo Fortran, tam všude lze získat potřebný vhled. Programování v systémech jako jsou MATLAB, Octave či Python s NumPy a SciPy je pro fyzika neméně potřebné, ale je jiné. Začátečníkovi se zde může mnohé ztrácet ze zřetele, zatímco programátor s jistou praxí už je uchvácen, kolik se tu toho dá zařídit tak jednoduše. My jsme si své začátky už skoro odbyli, takže si můžeme dovolit ke klasickému programovacímu jazyku přibrat i některý z těchto nástrojů – volíme Octave.

Na přednášce v části o numerických metodách se k Octavu také obracíme – většinu probíraných metod lze s jeho pomocí předvést na 2-3 řádcích. Zde na cvičení jej také přibereme a tady a teď v něm zpracujeme domácí úkoly. Studenti tím sice ztrácejí možnost úkoly v Octavu (nebo MATLABu) zpracovat sami, na oplátku mohou získat až trojí užitek – otrkat se s Octavem, inspirovat se postupem, porovnat výsledky.

Plná zadání úkolů jsou v poznámkách ke cvičení, zde jen hesly shrnujeme. Několik poznámek k Octavu přikládáme zde.

DÚ č. 1 – součet aritmetické posloupnosti (od 1 do n)

format compact; more off

% matematická optimalizace
n=input('Zadejte prosim cele cislo: ');
n,s=n*(n+1)/2

% Octave way
tic; n,s=sum(1:n); toc

% cyklus s podmínkou
tic; i=0; s=0;
while i<n, i=i+1; s=s+i; end
disp(sprintf('%7d%15.0f',n,s)); toc

% indexovaný cyklus
tic; s=0;
for i=1:n, s=s+i; end
disp(sprintf('%7d%15.0f',n,s)); toc

Příkaz format compact pro odmítnutí prostorem plýtvajícího výpisu je vhodný v MATLABu, v Octavu není třeba; příkaz more off vypne stránkovaný výpis, jehož vedlejším efektem je zde možné zpožďování výpisu. Příkaz input zajistí dotazem vstup hodnoty z klávesnice. Zápis výrazu nebo přiřazovacího příkazu zakončený čárkou nebo koncem řádku vede k výpisu výsledné hodnoty na obrazovku; středník naopak výpis potlačí. Default datový typ je 8bytový real (double), takže není třeba bojovat s integer přetékáním a dělením. Dvojice příkazů tic a toc změří čas uplynulý mezi jejich voláními. Dublet 1:n vytvoří pole délky n a prvky naplní aritmetickou posloupností od 1 do n, stejně jako by provedl triplet 1:1:n, kde na druhé pozici stojí explicitní krok. Funkce sum sečte hodnoty prvků vektoru. Cyklus while s podmínkou na začátku má v Octavu (nikoliv v MATLABu) k sobě do páru i cyklus do until s podmínkou na konci. Indexovaný cyklus postupně prochází všechny hodnoty vektoru, zde opět dubletu 1:n. Funkce sprintf se užívá k preciznímu formátování výstupu, zde celého (decimal) čísla na 7 míst (%7d) a reálného (floating-point) na 15 míst bez desetinné části (%15.0f). První i druhou variantu výpočtu provede Octave bleskově i pro velké hodnoty n (řekněme stovky milionů), na cyklech třetí a čtvrté varianty se zdrží, což u jazyka interpretujícího příkaz po příkazu nepřekvapuje. (MATLAB tyto ukázky provede bleskově všechny, zjevně si překládá větší sousta než Octave.)

DÚ č. 2 – Wallisův součin (aproximace nekonečného součinu)

format long;
nmax=50000000;
n=1:nmax;
2*prod((n*2)./(n*2-1).*(n*2)./(n*2+1))

Příkaz format long přepne na výpis desetinných míst v plné double přesnosti. Do vektoru n se připraví aritmetická posloupnost od 1 do nmax, na ní se vyčíslí členy Wallisova součinu a voláním funkce prod se všechny vynásobí. Pro násobení a dělení polí téhož tvaru prvek po prvku jsou užity operátory prvkového násobení .* a dělení ./, operátor * by se zde neúspěšně pokoušel o skalární součin. Při opomenutí středníku za třetím příkazem by se Octave začal ztěžka chystat k výpisu desítek milionů čísel. (To by se přerušilo stiskem klávesové kombinace Ctrl+C.)

DÚ č. 3 – 3D Lissajousův obrazec

fx = 3; fy = 5; fz = 7;
ax = 1; ay = 1; az = 1;
phx= 0; phy= 0; phz= 0;
nmax=1000;
t=linspace(0,pi*2,nmax+1);
x=ax*sin(fx*t+phx);
y=ay*sin(fy*t+phy);
z=az*sin(fz*t+phz);
plot(x,y); pause(2); plot(x,z); pause(2); plot(y,z); pause(2);
plot3(x,y,z);

Po úvodních inicializacích se funkcí linspace připraví pole t tvořené aritmetickou posloupností s předepsanými mezemi a daným počtem prvků. Následující trojicí příkazů se provede hromadné vyčíslení výrazů pro všechny prvky pole a všechny se přenesou do odpovídajích prvků polí x, y, z na levé straně přiřazovacích příkazů. Příkazy plot vynesou 2D průměty všech bodů, pause posečká uvedený počet sekund a příkaz plot3 vykreslí body ve třech dimenzích.

DÚ č. 4 – Mandelbrotova množina

xmin=-2; xmax=0.6; ymin=-1.2; ymax=1.2;
nx=1024/4; ny=768/4; nmax=100;
x=linspace(xmin,xmax,nx); y=linspace(ymin,ymax,ny); [c1,c2]=meshgrid(x,y);
c=c1+I*c2;
z=zeros(size(c));
niter=zeros(size(c));
for n=1:nmax
  mask=niter==0;
  z(mask)=z(mask).*z(mask)+c(mask);
  niter(niter==0 & abs(z)>2)=n;
end
hold on
for i=0:12, plot(c(niter>i),['.' num2str(mod(i,6))],'markersize',5); end
axis off % axis([xmin xmax ymin ymax]);

Dvojí volání funkce linspace připraví vektory pro pokrytí obdélníka <xmin,xmax>x<ymin,ymax> v komplexní rovině sítí nx x ny bodů. Matice c1 a c2 reálných souřadnic uzlů této sítě připraví dvojicí návratových hodnot funkce meshgrid a matici c komplexních souřadnic uzlů nachystá následující přiřazení zahrnující imaginární jednotku I (též i, j, J, 1i ad.). Pole nul na téže síti vytvoří funkce zeros s argumentem size(c) popisujícím tvar sítě: pole z bude zachycovat vývoj podle vztahu předepsaného pro Mandelbrotovu množinu a v poli niter se uchová index “únikové” iterace. V cyklu limitovaném na nmax iterací se přepočítává pole z (přepočet lze logickým polem mask omezit jen na ty prvky pole z, které zatím “neunikly”) a do prvků pole niter, které v n-té iteraci unikly z kruhu o poloměru 2, se zapisuje index této “únikové” iterace. V následujícím vizualizačním cyklu se vkreslí do jednoho obrázku (příkaz hold on) všechny komplexní body z pole c, které mají nenulový index niter, přičemž pomocí dotazu mod(i,6) na zbytek po dělení 6 se střídá šestice barev dostupných příkazu plot. Formátovací specifikace pro kreslení je zde řetězcem složeným z tečky (pro vyznačení bodu bodem o velikosti/markersize 5) a znaku s číselným vyjádřením barvy (funkce num2str konvertuje číslo na řetězec). Závěrečný příkaz ve dvou variantách upravuje osy právě zhotoveného obrázku.

DÚ č. 5 – prvočíselná funkce a prvočíselná dvojčata (počet prvočísel a prvočíselných dvojčat)

n=100000000;
p=primes(n);
d=p(p(1:end-1)+2==p(2:end));
plot(p,1:numel(p),d,1:numel(d));

Pole p prvočísel menších nebo rovných n obstará funkce primes, počet prvků pole (funkce numel) je pak hodnotou prvočíselné funkce pro dané n. Vlastně pro každé prvočíslo v poli p je jeho pořadové číslo (index) hodnotou prvočíselné funkce, a pro velká n stačí v obrázku vodorovně krokovat jen po prvočíslech, nikoliv po 1. Prvočíselná dvojčata se z pole prvočísel extrahují pomocí logického pole neboli relace vyjadřující souhrnně prvek po prvku definiční vlastnost prvočíselných dvojčat, tedy vzdálenost sousedů rovnou 2. Poslední prvek pole p nepokrytý logickým polem se do pole d nepřenese.

DÚ č. 6 – Newtonova metoda tečen a fraktál

xmin=-2; xmax=2; ymin=-2; ymax=2;
nx=500; ny=500; nitermax=50; eps=1e-10;
f=@(z) z.^3-1; df=@(z) z.^2*3; zroot0=1; zroot1=exp(pi*2/3*I); zroot2=exp(pi*4/3*I);
x=linspace(xmin,xmax,nx); y=linspace(ymin,ymax,ny); [z1,z2]=meshgrid(x,y);
z0=z1+I*z2;
z=z0;
n=zeros(size(z));
for niter=1:nitermax
  z=z-f(z)./df(z);
  n(n==0 & abs(z-zroot0)<eps)=niter;
end
plot(z0(n>0),'.','markersize',2);

Na síti z0 pokrývající obdélník v komplexní rovině (stejně jako v DÚ č. 4) se spustí Newtonova metoda tečen pro funkci komplexního argumentu \(f(z)=z^3-1\) s derivací \(f'(z)=3z^2\) a kořeny \(x_j=e^{(2/3)ij\pi}\), \(j=0,1,2\). Multiplikativní operace v anonymních funkcích přiřazených ukazatelům f a df jsou zapsány pomocí operátoru .^ tak, aby pro maticové argumenty pracovaly prvek po prvku (nikoliv ve smyslu maticového násobení). V cyklu se do prvků pole n zachycuje index niter Newtonovy iterace, pro který se příslušné prvky pole z poprvé přiblíží na vzdálenost eps k vybranému kořenu. Příkaz plot pak vykreslí komplexní body o souřadnicích v poli z0, pro které je příslušný prvek v poli n nenulový.

DÚ č. 7 – součet prvků inverzních Hilbertových matic

for n=1:12; disp(sprintf('%5d %20.15f',n,sum(sum(inv(hilb(n)))))); end

% pracnější alternativa
for n=1:12
  [mi,mj]=meshgrid(1:n,1:n);
  h=1./(mi+mj-1);
  hinv=h\eye(n);
  disp(sprintf('%5d %20.15f',n,sum(sum(hinv))))
end

Hilbertovu matici nxn vytvoří funkce hilb(n), její inverzi spočte funkce inv, řádek sloupcových součtů matice zajistí jedno sum a součet tohoto řádku druhé. To vše se pro výpis zformátuje funkcí sprintf na 5 znaků pro celočíselné n a 20 znaků (15 za desetinnou tečkou) pro reálný součet. V druhé variantě se připraví čtvercové nxn matice mi a mj s prvky osazenými hodnotami jejich řádkových a sloupcových indexů. Jedním přiřazením mezi poli se pak z hodnot indexů zkonstruují hodnoty prvků Hilbertovy matice h – zápis 1./ není double literál 1. a maticový operátor /, ale (i tak double) 1 a prvkový operátor ./. Inverzní matice hinv se nalezne řešením soustavy lineárních algebraických rovnic \(H.H^{-1}=I\) a dvojím voláním sum se sečtou všechny její prvky.

DÚ č. 8 – polynomická aproximace

x=[1 2 3 4]; f=[0 3 4 2]; xi=0:.1:5; hold off; plot(x,f,'o'); hold on;
for m=0:3, p=polyfit(x,f,m); plot(xi,polyval(p,xi)); end

Připraví se datové vektory x a f a síť xi pro vyčíslení aproximačního polynomu, vstupní data se vyznačí kolečky a do téhož obrázku (díky hold on) se přikreslí polynomy p na síti xi. Koeficienty polynomu p stupně m fitujícího body o souřadnicích v polích x a f se spočítají voláním funkce polyfit a hodnoty polynomů p na síti xi se vyčíslí voláním funkce polyval.

Odkazy

Domácí úkoly

student 1/Lazarus 2/Wallis 3/Lissajous 4/Mandelbrot 5/Erato. 6/Newton 7/Hilbert 8/LSQ 9/Morse 10a/Lorenz 10b/dopad
M. Be. ok ok ok ok ok ok ok ok ok ok ok + Z
T. Ha. ok ok ok                
V. Ha. ok ok ok ok ok ok ok ok ok ok ok + Z
D. Ko. ok ok ok       ok        
M. J. Ma. ok ok ok ok ok ok ok ok ok ok Z
A. Kř. ok ok ok ok ok ok ok ok ok ok ok + Z
M. On. ok ok ok ok ok ok ok ok ok ok ok + Z
N. N. An. ok ok ok ok ok ok ok ok ok ok ok + Z
P. Pe. ok ok ok ok ok ok ok ok ok ok ok + Z
P. Pr. ok ok                  
J. Se. ok 1/2                  
F. Sq.                      
J. Ša. ok ok ok ok ok ok ok ok ok ok ok + Z
O. Še. ok ok ok ok ok ok ok ok ok ok Z
D. Šk.                      
P. Šp. ok 1/2 ok       ok        
F. Šv. ok ok ok ok ok ok ok ok ok ok ok + Z
V. Ty. ok ok ok ok ok ok ok ok ok ok ok + Z

Vysvětlivky: ok = OK, 1/2 = doručeno, ale něco tomu chybí, ? = problém s doručením, Z = zápočet, prázdno = nedoručeno