Seminar Programovani v systemu Mathematica

Praha, CVUT, zari - rijen 2003

Zakladni stavebni bloky

Atomy

Head["Co je to?"]

Head[jmeno]

Atomy - cisla

{Head[0], Head[1/3], Head[1.0], Head[1], Head[1 + I]}

Pozor na podstatny rozdil mezi "1.0" a "1". Ve vstupnich datech numerickych vypoctu je dobre mit aspon jeden udaj ve formatu "Real".

Slozene vyrazy

h[(Subscript[e, 1]),  ... , (Subscript[e, n])]
h  ...  hlavicka,  (Subscript[e, i])  ...  elementy
Elementu muze byt nulovy pocet:

Date[]
Date[][[
      Range[3]]]
Date[][[{1, 2, 3}]]
Take[Date[], {1, 3}]
Drop[Date[], -3]

Priklad na vytvareni diferenci

Drop[{(Subscript[ee, 1]), (Subscript[ee, 2]), (Subscript[ee, 3])}, 1] - 
      Drop[{(Subscript[ee, 1]), (Subscript[ee, 2]), (Subscript[ee, 3])}, -1]

Priklady vyrazu

f[a, b], g[a][e1, f[3]],   2[4, 5]

Analyza vyrazu - hlavicka

e2 = f[a, b][][g[1], h[x, y]]
Head[e2]

Head[%]

Head[%]

Head[%]

Analyza vyrazu - elementy

e2

{Part[e2, 2], Part[e2, 1]}

{Part[h[x, y], 1], Part[h[x, y], 3]}

h[x, y][[1]]

Part[e2, 2] == e2[[2]]

Dalsi priklad

Take[{1, 2, 3}, {3}]
Part[{1, 2, 3}, 3]
{3}[[1]]
Take[{1, 2, 3}, {3}][[1]]
Part[{1, 2, 3}, {2, 3}]

Operatory a jejich vnitrni representace

Plus[a, b, c]

FullForm[a + b + c]

FullForm[{x, y, z}]

FullForm[f'[0]]

Derivative[1] je funkce Derivative v bode 1, jejimz argumentem je funkce f, takze Derivative[1][f] je jina funkce, zde v bode 0. Vsimneme si, ze derivace radu -1 je integral:

fd[x_] := x
Derivative[-1][fd][x]

Netrivialni priklady

FullForm[1 - 1]

FullForm[HoldForm[1 - 1]]

FullForm[a/b]

Priklad, kdy zabudovana funkce ma atribut HoldForm

Attributes[Plot]

Short[cheb = Table[ChebyshevT[i, x], {i, 10}], 3]

Plot[cheb, {x, -1, 1}]

Plot[Evaluate[cheb], {x, -1, 1}];

Definice funkce

FullForm[HoldForm[f[x] := (x^2)]]

Prepisovaci pravidlo (rewrite rule)

FullForm[
      HoldForm[log[a_*b_] :> log[a] + log[b]]]

Prepisovaci pravidla

Substituce, substitucni operator /.

h[x, y, z] /. z -> α

h[x, y, y] /. z -> α

Zde neni rozdilu mezi okamzitym a odlozenym prirazenim:

{f[2] /. f[x_] -> (x^2), f[2] /. f[x_] :> (x^2)}

?f

Definice

x//f znamena totez co f[x], s vyrazem typu TableForm se neda pocitat:

Clear[f];
f[x_]:=(x^2)
tab1 = 
      Table[{x, f[x]}, {x, 4}] // TableForm

2*tab1

Tato konstrukce celou vec napravi:

(tab1 = Table[{x, f[x]}, {x, 4}]) // TableForm

2*tab1

Priklad - definice logaritmu

Clear[log];
log[a_*b_] := 
      log[a] + log[b]

log[2*x*(y^2)]

Zde je pricina, proc neni rozvinut i log[(y^2)]: porovnavani vyrazů se děje syntakticky, a ne semanticky.

{FullForm[HoldForm[(y^2)]], FullForm[HoldForm[y*y]]}

Pridame tedy dalsi pravidlo:

log[x_^y_] := y*log[x]
log[2*x*(y^2)]

Zatim nase funkce log neni "listovatelna":

Attributes[log]

log[{a*b, c*d, c^d}]

Ale tento atribut muzeme snadno pridat:

SetAttributes[log, Listable];
Attributes[log]
log[{a*b, c*d, c^d}]

Ale pozor!!!

log[2] = 3;

log[2*x*(y^2)]

?log

Unset[log[2]]
?log

Priklad definice funkce se stejnou hlavickou a ruznym poctem argumentu

hu[x_] := 2*x
hu[x_, y_] := x + y
{hu[1], hu[3, 3]}

?hu

S jakymi objekty jsme v teto seanci zatim pracovali?

?Global`*

Programovaci styly

Uloha: vytvorit proceduru pro druhe mocniny listu

list = {1, -1, a, 2.5};

Proceduralni styl

mocproc[list_List] := 
      Module[{n, r, i}, 
n = 
          Length[list]; 
r = 
          Table[0, {n}]; 
Do[
          r[[i]] = 
            list[[i]]^2, {i, 1, 
            n}]; 
Return[r]
]

mocproc[list]

Strukturovany styl

mocstruk[list_List] := 
      Table[list[[i]]^2, {i, 
          Length[list]}]
mocstruk[list]

Mapovani funkce

Ryzi funkce

Function[x, x^2]
%[2]

Ve zkracenem zapisu

#1^2 &
%[2]

Zakladni funkce Map

Clear[f];
Map[f, {a, b, c}]

mocmap[list_List] := 
      Map[Function[x, (x^2)], list]
mocmap[list]

Ryzi funkci je samozrejme mozne pojmenovat:

ryzi = Function[x, x^2]

ryzi[2]
Map[ryzi[#] &, Range[5]]
Map[ryzi, Range[5]]

Automaticke mapovani

Attributes[Power]

list^2

Alternativni vyuziti listovatelnosti (1. element 1. listu umocnen na 1. element 2. listu ...)

list^Range[4]

Priklad efektivniho vyuziti listovatelnosti

Timing[Table[Sin[x], {x, 0, π, π/10000. }];]

Timing[Sin[Range[0, π, π/10000. ]];]

Rekursivni styl

Kazdy list je bud prazdny nebo ma prvni element a zbytek - ten zbytek muze byt prazdny list.

Klicova funkce pro tuto proceduru:

Prepend[{a, b, c}, e]

Predlohy: x_ prave jeden argument, x__ jeden a vice argumentu, x___ nula a vice argumentu

mocrek[{}] = {};
mocrek[{e_, rest___}] := 
      Prepend[mocrek[{rest}], e^2]
mocrek[list]

Prostudujte si fungovani krok po kroku.

Trace[mocrek[list], mocrek]

Dalsi konstrukce

Predlohy, pravidla, definice

(Patterns, rules, definitions)

_  nebo ekvivalentne Blank[] ... predloha, ktere vyhovuje vse (any expression matches it)

e3 = {1, x, f[], f[a], f[a, b], g[e4, e5]}

MatchQ[f[] g[a, b] h[1, 3, 4] 2[3, 4], _]

{MatchQ[h[a, b], _h], MatchQ[h[a, b], _a]}

?*Q

Vyhledavani vyhovujicich si vyrazu

Finding Matching Expressions

?Cases

Cases[e3, _]

Cases[e3, _f]

Cases[e3, _[_, _]]

Cases[e3, _[__]]

Cases[e3, f[___]]

Cases[e3, f[_, _]]

Podminene predlohy (conditional patterns)

Tvar: predloha/; podminka. Nejdrive se vyhledaji vyrazy vyhovujici predloze a pote ty, ktere splnuji podminku.

Cases[e6 = {f[1], f[0], f[-1], f[a]}, 
      f[x_] /; x ≥ 0]

Jiny mozny zapis:

Cases[e6, f[x_ /; x ≥ 0]]

e7 = Table[fff[j], {j, 9}];

Logicke operatory (ekvivalentni zapisy): a&&b  (And)         a||b (Or)

(2 ≤ x ≤ 5 || x > 6) == 
      Or[2 ≤ x ≤ 5, x > 6]

Cases[e7, 
      fff[x_] /; 2 ≤ x ≤ 5 || x > 6]

Predikaty

Predikat je funkce, ktera vraci True nebo False. Nasledujici zapisy jsou ekvivalentni:

{Cases[e6, f[x_ /; NonNegative[x]]], 
      Cases[e6, f[x_?NonNegative]], 
Cases[e6, 
        f[_?NonNegative]]}

Aplikace pravidel

Transformacni pravidla (v prvnim pripade chybi mezera za "/." a tudiz vyraz ".1" je vyhodnocen jako jmenovatel. V druhem pripade jiz mezera je.

Clear[f]

{f[1]/ .1 -> 2, f[1] /.  1 -> 2, 
      f[1] /. {1 -> 2}}

Predlohova promenna v pravidlech

Odlozene prirazeni:

(e8 = {f[0], f[2], g[0]}) /. f[x_] :> (x^2)

Okamzite prirazeni, v tomto pripade ekvivalentni odlozenemu:

e8 /. f[x_] -> (x^2)

Doporucuje se prednostne pouzivat odlozene pravidlo, prirazeni (delayed rule) leva strana  :>  prava strana pred okamzitym  -> . Duvody v nasledujicich dvou prikladech.

⚠ Varovny priklad 1

SeedRandom[13]; {x, x, x} /. 
      x -> Random[]
SeedRandom[13]; {x, x, x} /. x :> Random[]

⚠ Varovny priklad 2

z = 17;
{f[z_] :> (z^3), f[z_] -> (z^3)}

{f[5] /. f[z_] :> (z^3), 
      f[5] /. f[z_] -> (z^3)}

Opakovana substituce (ReplaceRepeated)

Clear[fac];
fac[4] /. {fac[0] -> 1, 
        fac[n_] :> n*fac[n - 1]}

Nahrazovani se deje pouze pouze jednim pruchodem listu pravidel. Nasledujici "//." (ReplaceRepeated) aplikuje pravidla tak dlouho, az se vyraz nemeni.

fac[4] //. {fac[0] -> 1, 
        fac[n_] :> n*fac[n - 1]}

⚠ Varovny priklad

Clear[fact];
TimeConstraint[
      fact[-3] //. {fact[1] -> 1, 
          fact[n_] :> n*fact[n - 1]}, 3]

Duvod spociva v tom, ze prvni pravidlo se neaplikuje nikdy, a tudiz druhe se aplikuje do nekonecna. Odpomoc spociva v pojistce, dotazu, na typ argumentu. V pripade, ze skutecny parametr nevyhovuje, vypocet se neprovede.

fact[-3] //. {fact[1] -> 1, 
        fact[n_Integer?Positive] :> n*fact[n - 1]}

fact[5] //. {fact[1] -> 1, 
        fact[n_Integer?Positive] :> n*fact[n - 1]}

Dalsi priklady definice funkci

Priklad 1

odm[x_] := "x je zaporne" /;  x < 0
odm[x_] := (Sqrt[x ]) /; 
        x ≥ 0
{odm[-3], odm[3]}

Priklad 2

sgn[x_] := 
      Which[
x < 0, -1, 
x == 0, 
        0, 
x > 0, 1]
Plot[sgn[x], {x, -1, 1}];

?sgn

Clear[sgn]   
?sgn

Remove[sgn]   
?sgn

Priklad 3 (Fibonacciova cisla)

Exisuje zabudovana funkce  Fibonacci:

Attributes[Fibonacci]

Fibonacci[Range[7]]

Definice A (naivni, vypocetni narocnost znacna):

Clear[fibA];
fibA[1] = 
        1;
fibA[2] = 
        1;
fibA[n_Integer?Positive] := 
      fibA[n - 1] + fibA[n - 2]

Table[fibA[n], {n, 25}]

?fibA

Definice B (Tato konstrukce umoznuje si zapamatovat jiz vypoctene hodnoty, podstatne rychlejsi):

fibB[1] = 1;
fibB[2] = 1;
fibB[n_Integer?Positive] := fibB[n] = 
        fibB[n - 1] + fibB[n - 2]

?fibB

Table[fibB[n], {n, 5}]

Narocnost na pamet je zanedbatelna. Vyzkousejte sina vice hodnot a pote zkuste

MaxMemoryUsed[]

Vsimneme si nyni na zapamatovanych hodnot:

?fibB

Priklad 4 (Programovani zalozene na pravidlech, vlastni definice linearni funkce)

Clear[f];
f[c_?NumericQ*x_] := c*f[x]
f[x_ + y_] := f[x] + f[y]

f[π + (Sqrt[2]) x + 3*a*(x^2)]

To je v poradku, ale  v poradku neni

f[2*π]

Odstranime tim, ze za pivotni hodnotu budeme uvazovat f[1].

f[c_?NumericQ] /; c =!= 1 := c*f[1]
{f[1], f[2*π]}

Uzitecne funkce Thread a MapThread

Clear[a, b, c];

Thread[h[a + b + c], Plus]
Thread[h[a*b*c], Times]

Thread[{a, b, c} -> {1, 2, 3}]

MapThread[Equal, {{x, y, z}, {1, 2, 3}}]

Clear[z];
MapThread[Equal, {{x, y, z}, {1, 2, 3}}]

Solve[%, {x, y, z}]

Funkcionalni programovani

Opakovani ryzich funkci (pure functions)

Map[Sqrt, {a, b, c}]

Map[Function[x, (x^2)], {a, b, c}]

#^2 &[3]

(#^2 &)'

(#1^2 + #2^2 &)[x, y]

Funkce Map a Apply

Map[f, {a, b, c}]
Apply[f, {a, b, c}]

Prvni argument obou funkci je funkce. Druhy argument je obvykle list, ale neni to pravidlem, viz nasledujici prikaz. Rozdil mezi obema funkcemi je ze Map ponecha hlavicku druheho argumentu a aplikuje prvni argument na elementy druheho argumentu zatimco Apply vymeni hlavicku druheho argumentu.

{Map[f, g[a]], Apply[f, g[a]]}

Map[Function[x, 1/x], e10 = {1, 2, 1/3, I}]

Ponekud mene prehledne, ale funkcne ekvivalentni zapis:

1/# &  /@  e10

⚠ Vyhnete se uziti funkce Sum (Σ) vsude, kde je to jen mozne. V nasledujicim pripade to mozne neni:

{(∑_(i = 1)^∞) (1/(i^2)), (∑_(i = 1)^∞) (1/(i^4))}

Apply[Plus, Range[5]]

Ukazka vyuziti pri vypoctu aritmetickeho prumeru a rozptylu

Clear[prumer];
prumer[data_?VectorQ] := 
      Apply[Plus, data]/Length[data] /; 
        Length[data] > 0
prumer[
      Range[5]]

mean = prumer;

rozptyl[data_?VectorQ] := 
      mean[(((data - mean[data]))^2)]
rozptyl[Range[5]]

Samozrejme funguje i pro symbolicke promenne.

prumer[{a, b, c}]

Vsimnete si zkratkoveho zapisu funkce Apply (@@):

?prumer

Funkce typu Nest a Fold

Nest[f, x0, 6]

NestList[f, x0, 4]

NestList[f, x0, 0]

Vypocet druhe odmocniny ze tri Newtonovou metodou s pocatecni hodnotou 1:

Nest[Function[x, (x + 3/x)/2], 1, 5]

My ale vetsinou pozadujeme numerickou hodnotu. Staci za kterymkoli cislem napsat tecku:

Nest[Function[x, (x + 3/x)/2], 1. , 5]

V predchozim pripade bylo nutne urcit pocet iteraci. Nasledujici funkce umoznuje najit pevny bod zobrazeni. Presnost druheho argumentu udava, kdy se proces ukonci (dve po sobe jdouci iterace se lisi o mene nez je presnost druheho argumentu). Zde na 40 mist.

FixedPoint[Function[x, (x + 3/x)/2], N[1, 40]]

Odvozena funkce NestList vraci vsechny iterace.

anu = NestList[Function[x, x/(1 + s)], c, 4]

To je vlastne diskontovana anuita pri urokove mire s. Soucasna hodnota je tedy:

presentvalue = Apply[Plus, anu]
presentvalue = 
      Apply[Plus, anu] // FullSimplify
FullSimplify[presentvalue = Apply[Plus, anu], 
      c > 0 && s > 0]

Prekvapeni, ze vysledek neni ve tvaru castecneho soouctu geometricke rady, a to ani pri pouziti volby predpokladu v poslednim pripade. Zkusime vyjadreni pomoci diskontniho faktoru a pouziti funkce Sum a dalsi finty:

rule1 = 
      Solve[v == 1/(1 + s), s] // 
        Flatten
Sum[c/(((1 + s))^i) //. rule1, {i, 0, 4}] // 
      FullSimplify
{Sum[(v^i), {i, 0, 4}], Sum[(v^i), {i, 0, n}]}

Sum[v^i, {i, 0, n}] /. n -> 4

Zahada, ze? Kdo znate pricinu, napiste mi.

Fold[g, x0, {a, b, c}]

FoldList[Plus, 0, {a, b, c}]

Priklad nahodne prochazky startujici z bodu 100

SeedRandom[13];
ceny1 = 
        FoldList[Plus, 100, 
          Table[Random[] - 0.5, {5000}]];
ListPlot[ceny1, PlotJoined -> True];

Totez, ale vyuzitim funkce NestList. Vyhodou je, ze nepotrebujeme predem vytvorit tabulku nahodnych cisel.

SeedRandom[13]; 
    ceny2 = NestList[Random[] - 0.5 + # &, 100, 
        5000];
ListPlot[ceny2];

Packages a upozorneni na funkci NMinimize

Vlastni package

Vytvoreni vlastniho package je snadne pri vyziti sablony (kostry) v adresari nize uvedenem, kde jsou i uzitecne informace Romana Maedra - v podstate vychazejici z jeho knih.

"C:\Program Files\Wolfram \
Research\Mathematica\5.0\AddOns\ExtraPackages\ProgrammingInMathematica\S\
keleton.m"

Dulezite je, aby nas package byl v adresari, kde ho Mathematica umi najit, pripadne je nutne popsat uplnou cestu. Seznam adresaru, ktere se prohledavaji pri hledani externiho souboru je uveden v promenne:

$Path

Neni-li nas Package v nekterem z vyse uvedenym adresaru, muzeme prislusny adresar pridat pomoci funkce

?AppendTo

podle vzoru

AppendTo[$Path, "C:\hurt\...."]

Dobre je prostudovat si tez nektery z jednodussich standardnich package (vybrat si funkci, kterou zname), pro mne napr. "C:\Program Files\Wolfram \
Research\Mathematica\5.0\AddOns\StandardPackages\Statistics\DataSmoothing.m\
"

Zabudovana funkce NMinimize

?NMinimize

Asociovana pole

mat=(({{"names", "BC", "EFG", "HIJ", "KLM"}, {"a", 1, 2, 3, 4}, {"b", 5, 6, 7, 8}, {"c", 9, 10, 11, 12}, {"d", 13, 14, 15, 16}}));

Clear[data, firstrow];
firstrow = 
      mat[[1]]
FullForm[firstrow]
fields = Map[data[#] &, firstrow]
matdrop = Drop[mat, 1] // Transpose
MapThread[Set, {fields, matdrop}]

fields

?fields

data["KLM"]

Select[mat, #[[2]] == 1 &]
Select[mat, #[[2]] > 7 &]
Select[
      fields, #[[4]] == 
          16 &]

?? fields

fields // Dimensions

{fields[[1, All]], fields[[All, 1]]}

Kernel estimate

Jadrovy odhad - neparametricky odhad hustoty
x ... hodnota, ve ktere odhad pocitame
k ... jadro
h ... sirka okenka - nutno vyzkouset experimentalne; cim vetsi, tim hladsi krivku dostaneme

kernel[x_,k_,h_,observ_]:=Module[{n=Length[observ],au},
		au=Map[k[(x-#)/h]&,observ];
		Apply[Plus,au]/(n*h)
		             ]

k[x_]:=Exp[-x^2/2]/Sqrt[2*Pi]

observ = Table[Random[] + 2*Random[], {50}];

Plot[kernel[x, k, 0.3, observ], {x, -1, 4}];

Odkazy

Mathematica version 5

\(\*""C:\Program Files\Wolfram Research\Mathematica\5.0\Documentation\English\Tour.nb"\)

Maeder, R.

Programming in Mathematica, Wolfram Educational Group, Version 1 2002

Maeder, R.

Programming in Mathematica, Addison-Wesley, Redwood, 2nd ed. 1990

Maeder, R. (obtizne)

The Mathematica Programmer II, Academic Press, San Diego, 1996

Shaw, W. T. , Tigg, J.

Applied Mathematica, Addison-Wesley, Reading (MA), 1994

Gaylord, R. J., Kamin, S. N., Wellin, P. R.

An Introduction to Programming with Mathematica, Springer Telos, New York 1996