Numerické výpočty II Jiří Zelinka Tento učební text vznikl za přispění Evropského sociálního fondu a státního rozpočtu ČR prostřednictvím Operačního programu Vzdělávání pro konkurenceschopnost v rámci projektu Univerzitní výuka matematiky v měnícím se světě (CZ.1.07/2.2.00/15.0203). 1 Obsah 1 Polynomy 4 1.1 Hornerovo schema a základní úkony s polynomy . . . . . . . . 5 1.2 Kořeny polynomů . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3 Interpolace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.4 Splajny . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.5 Bersteinovy polynomy a Bézierovy křivky . . . . . . . . . . . . 31 2 Derivace 38 2.1 Limity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.2 Symbolické derivování . . . . . . . . . . . . . . . . . . . . . . 41 2.3 Taylorův rozvoj . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.4 Numerické derivování . . . . . . . . . . . . . . . . . . . . . . . 44 3 Integrály 47 3.1 Symbolické integrování . . . . . . . . . . . . . . . . . . . . . . 47 3.2 Numerické integrování . . . . . . . . . . . . . . . . . . . . . . 50 4 Řady 54 4.1 Symbolické součty řad . . . . . . . . . . . . . . . . . . . . . . 54 4.2 Fourierovy řady . . . . . . . . . . . . . . . . . . . . . . . . . . 56 2 Úvod Tento text je zaměřen na základy numerické aproximace, zejména polynomiální a splajnovou interpolaci, kromě toho ovšem se budeme zabývat také diferenciálním a integrálním počtem funkcí jedné proměnné. Bližší informace o teoretickém základu, na němž tento text staví může čtšnář najít zejména ve skriptech [1] a [2], dále také v knihách [3] nebo [4]. 3 Kapitola 1 Polynomy Kropáčku, jak vám vyšly ty kořeny polynomu? Haló, slyšíte? Promiňte, pane profesore, trochu jsem se zamyslel. Na co jste se ptal? Ptám se, jaké máte kořeny. Moravské. Proč? Nebudeme se zda zabývat příliš teorií a základními pojmy jako např. stupeň polynomu apod. Literatura v tomto směru je velmi obsáhlá a čtenáři jistě nebude činit pražádné potíže si nějakou vyhledat. Nás budou zajímat především výpočty. První věc, na kterou bych rád upozornil, je nejednotnost zápisu polynomů v různé literatuře. Zpravidla se setkáme se zápisem P(x) = anxn + an−1xn−1 + · · · + a1x + a0 ale občas narazíme na opačné pořadí indexů, tedy na polynom ve tvaru P(x) = a0xn + a1xn−1 + · · · + an−1x + an. Pro konkrétní polynom, kdy koeficienty nejsou vyjádřeny obecně ale čísly, je to samozřejmě jedno, ale pokud chceme pro polynom použít nějaké tvrzení, je potřeba si dávat pozor, v jakém pořadí ta či ona věta uvádí pořadí koeficientů. V této příručce se budu držet prvního způsobu, tedy indexy u koeficientů budou stejné jako příslušná mocnina x. Pokud bude potřeba zvýraznit stupeň polynomu, budeme používat značení Pn. 4 1.1 Hornerovo schema a základní úkony s po- lynomy Na polynom budeme nahlížet jako na funkci – vrazíte do ní x a vypadne funkční hodnota podle daného vztahu. Základní věc, kterou tedy s polynomem potřebujeme dělat, je výpočet funkční hodnoty. Nejrychlejší metodou k tomu je tzv. Hornerovo schema. Nebudu uvádět jeho podrobné odvození, které by pro čtenáře bylo snadným cvičením. Hornerovo schema je založeno na dělení polynomu P polynomem prvního stupně P1(x) = x − c, tedy P(x) = (x − c) · Q(x) + A, (1.1) kde Q je podíl (polynom stupně o jedna menší než P) a A je zbytek po dělení, který musí mít stupeň menší než dělitel, tedy je to konstanta. Z uvedeného zápisu je jasné, že P(c) = A, tedy tento zbytek po dělení je současně hodnota polynomu P v bodě c. Z (1.1) se porovnáním koeficientů u stejných mocnin dají odvodit vztahy pro koeficienty polynomu Q, Předpokládáme-li polynom Q ve tvaru Q(x) = bn−1xn−1 + · · · + b1x + b0 dostáváme postupně bn−1 = an bn−2 = an−1 + c · bn−1 ... bk−1 = ak + c · bk ... b0 = a1 + c · b1 A = a0 + c · b0. Hornerovo schema zpravidla v literatuře nalezname coby následujíccí ta- bulku: an an−1 an−2 · · · a2 a1 a0 c bn−1 bn−2 bn−3 · · · b1 b0 A Pravidlo pro zapamatování výpočetního algoritmu je: „První číslo opíšeme (bn−1 = an), každé další dostaneme tak, že k číslu nad čarou připočteme c násobek předchozího čísla (bk−1 = ak + c · bk) . 5 Hornerovo schema se při ručních výpočtech nejčastěji používá k testování, zda dané číslo c je kořenem polynomu, v tom případě vyjde A = 0. Ale vidíme, že kromě hodnoty polynomu v bodě c dostáváme i podíl – polynom Q. Někoho možná napadne, co by se stalo, kdybychom Hornerovo schema použili se stejnou hodnotou c znovu, tentokrát na polynom Q, pak znovu na výsledný podíl a tak dál. Pro tento účel si označíme polynom Q jako Q1 a hodnotu A jakožto A0, v dalším kroku dostaneme podíl Q2 a hodnotu A1, a tak dál, takže obecně máme Qk(x) = (x − c) · Qk+1(x) + Ak. Hornerovo schema pak (symbolicky zkráceno) vypadá takto: P c Q1 A0 c Q2 A1 c Q3 A2 ... ... c An Pro polynom P pak dostáváme P(x) = (x − c)Q1(x) + A0 = = (x − c)((x − c)Q2(x) + A1) + A0 = = (x − c)2 Q2(x) + A1(x − c) + A0 = = (x − c)2 ((x − c)Q3(x) + A2) + A1(x − c) + A0 = = (x − c)3 Q3(x) + A2(x − c)2 + A1(x − c) + A0 = · · · = = An(x − c)n + An−1(x − c)n−1 + · · · + A1(x − c) + A0 Hodnoty An, . . . , A0 jsou tedy koeficienty polynomu P posunutého do bodu c. Toto vyjádření se dá také velmi dobře použít pro výpočet derivací polynomu P v bodě c, ale to bychom předbíhali. V Matlabu definujeme polynom pomocí vektoru jeho koeficientů od nejvyšší mocniny, takže příkazem >> P=[1 3 -2 0 5] definujeme polynom P(x) = x4 + 3x3 − 2x2 + 5, jak se dá snadno ověřit převodem polynomu na symbolický objekt: 6 >> poly2sym(P) ans = x^4 + 3*x^3 - 2*x^2 + 5 Pro výpočet hodnoty polynomu v bodě použijeme funkci polyval a můžeme ji aplikovat nejen na jednu hodnotu ale na vektor či matici hodnot, přičemž hodnoty polynomu se počítají pro každou služku zvlášť. Takže graf polynomu na intervalu můžeme získat třeba pomocí příkazů >> P=[1 3 -2 0 5]; >> x=-3:0.01:3; >> y=polyval(P,x); >> plot(x,y) Operace s polynomy Mezi základní úkony, které lze s polynomy provádět, patří sčítání, násobení skalárem, vzájemné násobení a dělení se zbytek. První dvě operace patrně nepotřebují další komentáře, jen je potřeba vědět, že pro sčítání polynomů v Matlabu musí mít příslušné vektory koeficientů stejnou délku, takže je nutné je doplnit v případě potřeby nulami. Pro násobení polynomů se používá funkce conv, kde jako vstupní parametry zadáme polynomy, které chceme násobit. Dělení polynomu se zbytkem se provádí podobně, jako je tomu u celých čísel. Nejznámějším použitím dělení polynomu se zbytkem je Eukleidův algoritmus pro hledání největšího společného dělitele dvou polynomů. Algerotmus je všeobecně známý, proto jej zde nebudeme uvádět. Pro dělení se zbytkem v Matlabu je možné použít funkci deconv, která má dva vstupní a dva výstupní parametry. Na výstupu dostáváme podíl a zbytek po dělení. Například při dělení polynomu x4 + 3x3 − 4x + 1 polynomem x2 + 1 dostáváme >> P=[1 3 0 -4 1]; >> Q=[1 0 1]; >> [D,R]=deconv(P,Q) D = 1 3 -1 7 R = 0 0 0 -7 2 Zbytek po dělení R má formálně stejný stupeň jako dělenec P, aby platila rovnost P=D*Q+R. Tuhle skutečnost je potřeba brát v potaz při některých výpočtech v Matlabu, například právě u zmíněného Eukleidova algoritmu. 1.1.1 Derivace polynomu V této části trochu předběhneme učivo matematického kursu, protože derivace polynomu se nám bude hodit v dalším výkladu a navíc pro polynomy se počítá poměrně jednoduše. Obecně lze říci, že derivace funkce udává, jak rychle se funkce mění. Tedy pokud funkce hodně rychle roste, má derivace velké hodnoty, pokud se funkce na nějakém intervalu nemění (je konstantní), je tam její derivace nulová a tak podobně. Derivace je kladná tam, kde funkce roste, záporná, když funkce klesá, v bodech, kde má funkce lokální minimum nebo maximum, je derivace nulová. Pro polynom P(x) = anxn + an−1xn−1 + · · · + a1x + a0 se jeho derivace vypočítá podle vztahu P (x) = n · anxn−1 + (n − 1) · an−1xn−2 + · · · + 2 · a2x + a1 V Matlabu se derivace vypočítá příkazem polyder, takže pokud si chceme zobrazit polynom spolu s jeho derivací, můžeme to udělat například takto: >> P=[3 -4 -12 0 5]; >> DP=polyder(P) DP = 12 -12 -24 0 >> x=-2:0.01:3; >> y=polyval(P,x); >> dy=polyval(DP,x); >> z=zeros(size(x)); >> plot(x,y,x,dy,x,z,’k--’) 8 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3 −100 −50 0 50 100 150 Modře je zobrazen polynom, zeleně jeho derivace. 1.2 Kořeny polynomů Najít kořen daného polynomu je jednou ze základních matematických úloh. K ověřování, zda dané číslo je nebo není kořenem polynomu zpravidla používáme Hornerovo schema, existují ale další nástroje, které nám mohou hledání kořenů usnadnit. Například pokud máme polynom s celočíselnými koeficienty a zajímají nás racionální kořeny, tedy kořeny ve tvaru p q , kde p je celé číslo, q je přirozené číslo a p a q jsou nesoudělná, tak se dá lehce zjistit, že koeficient u nejvyšší mocniny polynomu (tedy an) musí být dělitelný číslem q, zatímco absolutní člen (t.j. a0) musí být dělitelný p. Další pomůckou může být stanovení oblasti, kde všechny kořeny leží. K tomu nám může pomoci například tvrzení, které lze najít i s důkazem v [1]. Nechť A = max (|a0|, . . . , |an−1|) , B = max (|a1|, . . . , |an|) , kde ak, k = 0, 1, . . . , n, a0an = 0, jsou koeficienty polynomu P, Pak pro 9 všechny kořeny xk, k = 0, 1, . . . , n, polynomu P platí 1 1 + B |a0| ≤ |xk| ≤ 1 + A |an| . (1.2) Všechny kořeny tedy v komplexním oboru leží v mezikruží, jehož hranice jsou dány uvedenou nerovností. Uvedené hranice, zejména pak horní hranice, jsou poměrně hodně hrubě odhadnuty. Například pro polynom, jehož kořeny jsou čísla od jedné do šesti, dostaneme následující hranice: >> P=poly(1:6) P = Columns 1 through 6 1 -21 175 -735 1624 -1764 Column 7 720 >> n=length(P); >> A=max(abs(P(2:n))) A = 1764 >> B=max(abs(P(1:n-1))) B = 1764 >> H=1+A/P(1) % horni hranice H = 1765 >> D=1/(1+B/P(1)) % dolni hranice D = 5.6657e-04 Existují další kriteria, která v některých případech mohou poskytnout lepší odhad velikosti absolutní hodnoty kořenů, uveďme aspoň některá: Pro všechny kořeny xk, k = 0, 1, . . . , n, polynomu P platí |xk| ≤ max    1, n−1 j=0 aj an    |xk| ≤ 2 max an−1 an , an−2 an , . . . , n a0 an 10 |xk| ≤ max a0 an , 1 + a1 an , . . . , 1 + an−1 an . Dalším zjednodušením při hledání kořenů polynomu může být odstranění násobných kořenů. Obecně platí, že každý polynom lze zapsat ve tvaru P(x) = an(x − x1)n1 · · · (x − xk)nk , kde x1,. . . ,xk jsou vzájemně různá komplexní čísla, n1,. . . ,nk jsou jejich násobnosti a n1 + . . . + nk = n = st(P). Platí také, že pokud je kořen xi násobnosti ni > 1, pak tento kořen je i kořenem derivace polynomu P s násobností ni − 1. Odtud plyne, že snížit násobnost všech kořenů na 1 lze tak, že určíme největší společný dělitel D polynomu P a jeho derivace P . Kořeny tohoto největšího společného dělitele jsou právě ty kořeny původního polynomu, které mají násobnost větší než jedna a jejich násobnost v děliteli D je o jedna menší než v polynomu P. Vydělením polynomu P dělitelem D získáme tedy polynom, který má stejné kořeny jako P, ale všechny jsou násobnosti 1. Z numerického hlediska je potřeba poznamenat, že tento postup v praxi funguje dobře pro polynomy se celočíselnými koeficienty, které se pohybují v „rozumných mezích. Ale i v jednoduchých případech je potřeba postupovat obezřetně: >> format long >> P=poly([1 1 1 1]) P = 1 -4 6 -4 1 >> DP=polyder(P) DP = 4 -12 12 -4 >> roots(P) ans = 1.000217162530750 0.999999969335793 + 0.000217131857745i 0.999999969335793 - 0.000217131857745i 0.999782898797672 >> roots(DP) ans = 11 1.000004930958320 + 0.000008540746793i 1.000004930958320 - 0.000008540746793i 0.999990138083364 Vidíme, že pokud bychom posuzovali shodnost kořenů u uvedeného polynomu a jeho derivace, přičemž ani nepožadujeme příliš velkou přesnost, tak nejenže ke shodě nedochází, ale kořeny ani nejsou násobné. Nicméně pokud bychom hledali největší společný dělitel, postup povede ke správnému výsledku: >> [Q,R]=deconv(P,DP) Q = 0.250000000000000 -0.250000000000000 R = 0 0 0 0 0 >> D=Q % D je nejv.spol.delitel D = 0.250000000000000 -0.250000000000000 >> roots(D) ans = 1 Vidíme, že výsledek je přesný, a to i navzdory numerickým nepřesnostem během výpočtu. Ukažme si ještě jeden příklad, kde budou dva násobné kořeny blízko sebe: >> P=poly([0.9 0.9 1.1 1.1 1.1]) P = 1.0000 -5.1000 10.3800 -10.5380 5.3361 -1.0781 >> DP=polyder(P) DP = 5.0000 -20.4000 31.1400 -21.0760 5.3361 >> [Q1,R1]=deconv(P,DP) Q1 = 0.2000 -0.2040 R1 = 0 0 -0.0096 0.0298 -0.0306 0.0105 >> R1(1:2)=[] % odstranime pocatecni nuly R1 = 12 -0.0096 0.0298 -0.0306 0.0105 >> [Q2,R2]=deconv(DP,R1) Q2 = -520.8333 510.4167 R2 = 1.0e-11 * 0 0 -0.1766 0.2515 -0.0769 V tomto místě je nutné usoudit, že zbytek po dělení je ve skutečnosti nulový, tudíž nějvětší společný dělitel je předchozí zbytek, tedy polynom R1. >> D=R1 D = -0.0096 0.0298 -0.0306 0.0105 >> P0=deconv(P,D) P0 = -104.1667 208.3333 -103.1250 >> roots(P0) ans = 1.1000 0.9000 >> format long >> roots(P0) ans = 1.100000000001784 0.899999999998279 Vidíme, že výsledek je opět poměrně přesný. Kořeny původního polynomu Matlab spočítá s daleko větší chybou: >> roots(P) ans = 1.100021372535598 + 0.000037011183255i 1.100021372535598 - 0.000037011183255i 1.099957254925198 0.900000434848828 0.899999565154781 13 1.2.1 Separace reálných kořenů Ukážeme si algoritmus, s jehož pomocí je možné přesně určit počet reálných polynomů v daném intervalu, pokud předpokládáme, že všechny reálné kořeny jsou jednoduché. Tento algoritmus je založen na konstrukci tzv. Sturmovy posloupnosti, která se definuje následovně: Posloupnost reálných polynomů P = P0, P1, . . . , Pm se nazývá Sturmovou posloupností příslušnou polynomu P, jestliže 1. Všechny reálné kořeny polynomu P0 jsou jednoduché. 2. Je-li x0 reálný kořen polynomu P0, pak signP1(x0) = −signP0(x0). 3. Jestliže x0 je reálný kořen polynomu Pi, platí Pi+1(x0)Pi−1(x0) < 0, pro i = 1, 2, . . . , m − 1. 4. Poslední polynom Pm nemá reálné kořeny. Sturmovu posloupnost lze zkonstruovat poměrně snadno, Pokud předpokládáme, že výchozí polynom P = P0 nemá reálné kořeny, stačí položit P1 = −P0, abychom zaručili splnění druhé vlastnosti. Třetí vlastnost z definice dostaneme, pokud polynom Pi+1 vezmeme jako záporně vzatý zbytek po dělení Pi−1 : Pi,tedy Pi−1 = Pi · Q − Pi+1 Pro kořen x0 polynomu Pi tedy máme Pi−1(x0) = −Pi+1(x0), přičemž daný výraz nemůže být nulový, jinak bychom indukcí dostali, že x0 je kořenem P0 i P1, což je spor s tím, že kořeny P0 jsou jednoduché. Konstrukce založená na dělení se zbytkem zaručuje, že stupně polynomů v posloupnosti postupně klesají. Poslední polynom Pm je zpravidla konstantní, ale je možné konstrukci ukončit i dříve, pokud víme, že polynom nemá reálné kořeny, např. pro polynom x2 + 1. Počet kořenů v daném intervalu pak lze zjistit pomocí Sturmovy věty: 14 Počet reálných kořenů polynomu P v intervalu [a, b] je roven W(b) − W(a), kde W(x) je počet znaménkových změn ve Sturmově posloupnosti P0(x), . . . , Pm(x) v bodě x (z níž jsou vyškrtnuty nuly). Předpokládám, že pojem počet znaménkových změn je natolik intuitivně jasný, že nepotřebuje definici. Důkaz Sturmovy věty je založen na faktu, že k navýšení či snížení počtu znaménkových změn může dojít jen při přechodu přes kořen některého z polynomů ve Sturmově posloupnosti. Ze druhé vlastnosti v definici Sturmovy posloupnosti skutečně plyne, že pokud x roste, pak při přechodu přes kořen P = P0 jsou dvě možnosti: buď přecházíme ze záporných hodnot polynomu do kladných, v tom případě polynom P0 roste, jeho derivace je tedy v okolí kořene kladná, tudíž P1 je tam záporný a přibude jedna znaménková změna. Nebo P0 v kořenu klesá, takže P1 je kladný a také přibude jedna znaménková změna. x − h x x + h P0 − 0 + P1 − − − W(x) 0 0 1 x − h x x + h P0 + 0 − P1 + + + W(x) 0 0 1 Při přechodu přes kořen polynomu Pi pro i > 0 jsou celkem čtyři možnosti, při žádné z nich však podle třetí vlastnosti z definice nedochází ke změně počtu znaménkových změn: x − h x x + h Pi−1 − − − Pi − 0 + Pi+1 + + + W(x) 1 1 1 x − h x x + h Pi−1 + + + Pi − 0 + Pi+1 − − − W(x) 1 1 1 x − h x x + h Pi−1 − − − Pi + 0 − Pi+1 + + + W(x) 1 1 1 x − h x x + h Pi−1 + + + Pi + 0 − Pi+1 − − − W(x) 1 1 1 K navýšení počtu znaménkových změn tak dochází jenom při přechodu přes kořen polynomu P0, odkud bezprostředně plyne tvrzení věty. Příklad 1. Zjistíme počet kladných a záporných kořenů polynomu x4 −4x2 + 8x − 2. Nejprve sestrojím Sturmovu posloupnost: >> P0=[1 -4 0 8 -2]; >> P1=-polyder(P0) 15 P1 = -4 12 0 -8 Protože nás zajímají jen znaménka a jejich změny, je možné polynom vydělit nebo vynásobit kladným číslem. To se hodí hlavně při ručním počítání, když se chceme vyhnout složitějším zlomkům. Pak dál pokračujeme v konstrukci Sturmovy posloupnosti: provedeme dělení se zbytkem, odstraníme nulové koeficienty a polynom opět můžeme vydělit kladným číslem: >> P1=P1/4 P1 = -1 3 0 -2 >> [Q,R]=deconv(P0,P1) Q = -1 1 R = 0 0 -3 6 0 >> P2=-R/3; >> P2(1:2)=[] P2 = 1 -2 0 Celou činnost opakujeme, dokud nedostaneme konstantní polynom: >> [Q,R]=deconv(P1,P2) Q = -1 1 R = 0 0 2 -2 >> P3=-R/2; >> P3(1:2)=[] P3 = -1 1 >> [Q,R]=deconv(P2,P3) Q = -1 1 R = 0 0 -1 >> P4=1; 16 Podle předchozího textu je horní hranice pro absolutní hodnotu všech kořenů rovna 9, takže všechny kořeny leží v intervalu [−9, 9]. Při ručním počítání je jednodušší určit znaménko limity hodnoty v ±∞ podle parity nejvyšší mocniny a znaménka jejího koeficientu. Počet znaménkových změn v Matlabu určíme pomocí příkazů: >> x=-9; >> [polyval(P0,x),polyval(P1,x),polyval(P2,x),... polyval(P3,x),polyval(P4,x)] ans = 9403 970 99 10 1 >> x=0; >> [polyval(P0,x),polyval(P1,x),polyval(P2,x),... polyval(P3,x),polyval(P4,x)] ans = -2 -2 0 1 1 >> x=9; >> [polyval(P0,x),polyval(P1,x),polyval(P2,x),... polyval(P3,x),polyval(P4,x)] ans = 3715 -488 63 -8 1 Při ručním počítání stačí určovat znaménka a výsledky můžeme přehledně umístit do tabulky: x P0(x) P1(x) P2(x) P3(x) P4(x) W(x) −∞ + + + + + 0 0 − − 0 + + 1 ∞ + − + . + 4 Celkem tedy máme 4 kořeny, z toho je jeden záporný a tři kladné. 1.3 Interpolace Dobrý den, studenti, dnes nás čeká interpolace. Co myslíte, Kropáčku, že to je? To je něco v parlamentu, ne? To jste si spletl s interpelací. Aha, no vlastně. A není to pátrání Interpolu? 17 Problém interpolace patří do teorie aproximace. Zpravidla se snažíme aproximovat nějakou funkci, z níž známe jenom hodnoty v diskrétních bodech, někdy je dokonce můžeme mít nepřesné. Lagrangeova interpolace Předpokládejme, že jsou dány body xi, i = 0, 1, . . . , n, xi = xk pro i = k a hodnoty funkce f v těchto bodech: f(xi) = fi, i = 0, 1, . . . , n. Hledáme polynom Pn stupně nejvýše n takový, že Pn(xi) = fi, i = 0, 1, . . . , n. Body xi, i = 0, 1, . . . , n, xi = xk pro i = k, budeme nazývat uzly, polynom Pn interpolační polynom. Platí následující tvrzení, které se dá poměrně snadno dokázat: Pro (n + 1) daných dvojic čísel (xi, fi), i = 0, 1, . . . , n, xi = xk pro i = k, existuje nejvýše jeden polynom Pn stupně menšího nebo rovného n takový, že Pn(xi) = fi, i = 0, 1, . . . , n. Pro důkaz existence interpolačního polynomu použijeme Lagrangeovu konstrukci, podle níž se interpolační polynom nazývá Lagrangeův: Položme li(x) = (x − x0) . . . (x − xi−1)(x − xi+1) . . . (x − xn) (xi − x0) . . . (xi − xi−1)(xi − xi+1) . . . (xi − xn) = n j=0 j=i (x − xj) (xi − xj) . Je zřejmé, že li je polynom stupně n a li(xj) = 0 pro i = j 1 pro i = j. Definujme nyní Lagrangeův interpolační polynom Pn vztahem: Pn(x) = l0(x)f0 + l1(x)f1 + . . . + ln(x)fn = n i=0 li(x)fi. 18 Snadno se ověří, že tento polynom splňuje dané interpolační podmínky a je polynomem stupně nejvýše n. Polynomy li se nazývají Lagrangeovy fundamentální polynomy a tvoří bázi prostoru polynomů stupně nejvýše n. Podívejme se, jaké je výpočetní složitost Konstrukce Lagrangeova interpolačního polynomu. Při konstrukci jednoho fundamentálního polynomu začínáme polynomem prvního stupně x − x0, který postupně násobíme členy x − xj, stupeň polynomu se postupně zvětšuje a pro jeho násobení potřebujeme řádově tolik operaci, jaký je jeho stupeň. Celkový počet operací pro konstrukci li tedy dostaneme jako součet konečné aritmetické řady, což je řádově n2 operací. Abychom sestrojili všechny fundamentální polynomy potřebujeme tedy řádově n3 operací. Při ruční konstrukci interpolačního polynomu ovšem zjistíme, že spoustu operací provádím opakovaně, takže při vhodném postupu by bylo možné konstrukci urychlit. A skutečně, optimální konstrukce Lagrangeova interpolačního polynomu je založena na funkci ωn+1(x) = (x − x0) . . . (x − xn) = n i=0(x − xi). Pro určení funkce ωn+1 potřebujeme řádově také n2 operací, ale tuto funkci konstruujeme jen jedenkrát. Čitatel fundamentálního polynomu li získáme dělením ωn+1(x) : (x − xi), k čemuž se dá využít Hornerovo schema, které je co se týče počtu operací lineární. Dále se dá odvodit, že jmenovatel fundamentálního polynomu li je roven derivaci funkce ωn+1 v bodě xi. Pro derivování polynomu je potřeba také řádově n operací a pro výpočet hodnoty derivace v bodě jakbysmet. Celkově tedy pro jeden fundamentální polynom potřebujeme jen lineární množství operací, pro všechny fundamentální polynomy je to pak řádově n2 operací. Matlabovská funkce pro konstrukci Lagrangeova interpolačního polynomu pak může vypadat následovně: function [lip,lfp] = laintpol2(uzly,ff) % function [lip,lfp] = laintpol(uzly,ff) % Lagrangeuv interpolacni polynom % uzly,ff - radkove vektory % lip - Lagr. interpol. polynom % lfp - matice fundamentalnich polynomu om=poly(uzly); % funkce omega omd=polyder(om); % omega’ n=length(uzly)-1; 19 lip=zeros(1,n+1); lfp=[]; for ii=0:n i1=ii+1; xi=uzly(i1); citatel=deconv(om,[1,-xi]); jmenovatel=polyval(omd,xi); li=1/jmenovatel*citatel; % fundamentalni polynom lfp=[lfp;li]; end % konec cyklu lip=ff*lfp; end % konec funkce Jako nepovinný druhý výstupní parametr funkce vrací fundamentální polynomy řádcích matice. Pro vyzkoušení programu zkusíme interpolovat hodnoty funkce sin. Nejprve definujeme uzly a určíme funkční hodnoty: >> uzly=[1 2 3 4 6]; >> ff=sin(uzly); >> [lip,lfp] = laintpol2(uzly,ff); Protože máme spočtené i Lagrangeovy fundamentální polynomy, můžeme si je prohlédnout: >> l1=lfp(1,:); >> l2=lfp(1,:); >> l2=lfp(2,:); >> l3=lfp(3,:); >> l4=lfp(4,:); >> l5=lfp(5,:); >> xx=0:0.01:7; >> L1=polyval(l1,xx); >> L2=polyval(l2,xx); >> L3=polyval(l3,xx); >> L4=polyval(l4,xx); >> L5=polyval(l5,xx); >> plot(xx,[L1;L2;L3;L4;L5]) 20 >> grid on >> axis([0,7,-5,5]) 0 1 2 3 4 5 6 7 −5 −4 −3 −2 −1 0 1 2 3 4 5 Nakonec si necháme vykreslit i interpolační polynom spolu s funkcí sin pro porovnání: >> Px=polyval(lip,xx); >> fx=sin(xx); >> plot(xx,Px,xx,fx,’--’,uzly,ff,’*’) 21 0 1 2 3 4 5 6 7 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 Vidíme, že polynom funkci aproximuje poměrně dobře, chyba je větší v oblasti, kde jsou uzly dále od sebe. Chyba se samozřejmě zvětšuje vně intervalu obsahující uzly. Hlavní nevýhodou Lagrangeovy konstrukce je ten fakt, že v případě, když potřebujeme přidat uzel, musíme přepočítat všechny fundamentální polynomy. Proto je v některých případech výhodnější Newtonova konstrukce, kterou si teď odvodíme. Je několik způsobů, jak to udělat, my použijeme tzv. iterovanou interpolaci. Iterovaná a Newtonova interpolace Označme Pi,i+1,...,i+k interpolační polynom na uzlech xi,. . . ,xi+k, tedy polynom stupně k. Pro k = 0 dostáváme polynom nultého stupně Pi, tedy se jedná o konstantní polynom, který je všude roven hodnotě fi. Polynom Pi,i+1,...,i+k sestrojíme z polynomů nižších stupňů následujícím způsobem: Nechť Pi,...,i+k−1 a Pi+1,...,i+k jsou interpolační polynomy na příslušných uzlech, oba stupně k − 1. Pak polynom Pi,i+1,...,i+k získáme ze vztahu Pi,...,i+k(x) = 1 xi+k − xi Pi+1,...,i+k(x) Pi,...,i+k−1(x) x − xi+k x − xi = = 1 xi+k − xi [Pi+1,...,i+k(x) · (x − xi) − Pi,...,i+k−1(x) · (x − xi+k)]. Dosazením jednotlivých uzlů do této formule snadno zjistíme, že výsledný polynom skutečně splňuje interpolační podmínky. 22 Uvedený vztah trochu upravíme: Pi,...,i+k(x) = 1 xi+k − xi [Pi,...,i+k−1(x) · ((xi+k − xi) − (x − xi)) + + Pi+1,...,i+k(x) · (x − xi)] = Pi,...,i+k−1(x) + Pi+1,...,i+k(x) − Pi,...,i+k−1(x) xi+k − xi (x − xi) Rozdíl polynomů v čitateli zlomku je samozřejmě polynom stupně k − 1, označme jej například R. Protože jak polynom Pi,...,i+k−1 tak i Pi+1,...,i+k mají stejné funkční hodnoty v uzlech xi+1,. . . ,xi+k−1, jsou tyto uzly kořeny polynomu R. Protože těchto uzlů je k − 1 lze polynom R zapsat ve tvaru R(x) = a · (x − xi+1) · · · (x − xi+k−1), kde a je koeficient u nejvyšší mocniny x, tedy u xk−1 . Tento koeficient je ale roven rozdílu koeficientů u nejvyšší mocniny polynomů Pi+1,...,i+k a Pi,...,i+k−1. Celkem dostáváme Pi,...,i+k(x) = Pi,...,i+k−1(x) + a xi+k − xi (x − xi) · · · (x − xi+k−1). (1.3) Interpolační polynom na uzlech xi,. . . ,xi+k tedy dostaneme z interpolačního polynomu na uzlech xi,. . . ,xi+k−1 přidáním dalšího členu, který je roven součinu kořenových činitelů k−1 j=0 (x−xi+j) vynásobenému koeficientem u nejvyšší mocniny. Označme tento koeficient jako f[xi, . . . , xi+k]. Je jasné, že f[xi] = fi a dále ze vztahu (1.3) dostáváme f[xi, . . . , xi+k] = f[xi+1, . . . , xi+k] − f[xi, . . . , xi+k−1] xi+k − xi . Způsob výpočtu těchto koeficientů jim dává také název - označujeme je jako poměrné diference. Jejich výpočet pro všechny uzly lze provést pomocí ná- 23 sledujícího schematu: x0 f[x0] = f0 > f[x0, x1] x1 f[x1] = f1 > f[x0, x1, x2] > f[x1, x2] ... x2 f[x2] = f2 f[x0, . . . , xn] ... ... ... ... ... > f[xn−2, xn−1, xn] > f[xn−1, xn] xn f[xn] = fn Postupným použitím vztahu (1.3) pak dostáváme P0,...,n(x) = f[x0] + f[x0, x1](x − x0)+ + f[x0, x1, x2](x − x0)(x − x1)+ + . . . + f[x0, . . . , xn](x − x0) · · · (x − xn−1). (1.4) Tato konstrukce se nazývá Newtonův interpolační polynom. Výsledek samozřejmě musí být totožný s Lagrangeovým interpolačním polynomem. Hlavní výhodou Newtonovy konstrukce ale je snadné přidávání dalších uzlů – stačí vypočítat další řádek v tabulce pro výpočet poměrných diferencí a k předchozímu interpolačnímu polynomu přidat součin příslušných kořenových činitelů násobený výslednou poměrnou diferencí. Pro interpolaci v Matlabu je možné použít také jeho vlastní funkci polyfit, která je ale trochu obecnější a umožňuje najít také aproximaci dat polynomem nižšího stupně pomocí metody nejmenších čtverců. Intrpolovat lze samozřejmě take v Sage. Abychom dostali jako výsledek racionální polynom, použijeme pro příklad jiná data než výše uvedená: sage: P = PolynomialRing(QQ, ’x’) sage: P.lagrange_polynomial([(0,1),(2,2),(3,-2),(-4,9)]) -23/84*x^3 - 11/84*x^2 + 13/7*x + 1 Dají se zde určit i poměrné diference pro Newtonův interpolační polynom: sage: P.divided_difference([(0,1),(2,2),(3,-2),(-4,9)]) [1, 1/2, -3/2, -23/84] 24 Předpokádám, že čtenáři nečiní žádné potíže si ověřit, že z vypočítaných nodnot pomocí Newtonovy konstrukce získáme stejný polynom jako v prvním případě. 1.4 Splajny Pane profesore, tady v těch materiálech máte pravopisnou chybu. Opravdum, Kropáčku? Ukažte. Tady máte „splyne s měkkým i. Vypadá to, Kropáčku, že kvůli vám budu muset začít používat český přepis slova spline. Teď zase trochu předběhneme. budeme totiž opět potřebovat derivace, o kterých budeme podrobněji mluvit později, ovšem výklad o splajnech dobře navazuje na interpolaci, takže jsem se rozhodl jej zařadit už sem. Snad to nebude příliš vadit, protože tato část bude spíše jen doplňková a informativní. V češtině už se pro funkce, o nichž teď budeme mluvit, zaběhl počeštěný výraz splajn. Ve starší literatuře je možné narazit na anglický termín spline, výslovnost je ovšem stejná jako u počeštěného výrazu. Splajny patří mezi interpolační techniky, protože většinou procházejí přesně zadanými hodnotami. Existují ale i tzv. vyhlazovací splajny, o těch tady v tuto chvíli taktně pomlčíme. Splajn se zpravidla definuje jako po částech polynomiální funkce, která je spojitá do určitého řádu. Nejjednodušším příkladem splajnu je spojitá po částech lineární funkce, která prochází zadanými body v rovině. Přesněji řečeno předpokládejme opět, že jsou dány body xi, i = 0, 1, . . . , n, xi = xk pro i = k a hodnoty funkce f v těchto bodech: f(xi) = fi, i = 0, 1, . . . , n. Lineární spojitý splajn je funkce s, která je spojitá, lineární na každém intervalu [xi, xi+1], i = 0, . . . , n − 1 a platí s(xi) = fi, i = 0, . . . , n. Na následujícím obrázku můžeme vidět příklad takové funkce: 25 1 2 3 4 5 6 7 0 0.5 1 1.5 2 2.5 3 Nejčastěji se setkáme s kubickými splajny. Jedná se o funkce, které jsou po částech polynomy třetího stupně a jsou spojité do řádu dva, jsou tedy spojité a mají spojité i první a druhé derivace, Uvedeme radši opět přesnou definici: Nechť jsou dány body xi, i = 0, 1, . . . , n, xi = xk pro i = k (zvané uzly) a hodnoty funkce f v těchto bodech: f(xi) = fi, i = 0, 1, . . . , n. Kubický splajn je funkce s, která je spojitá včetně první a druhé derivace a je kubickým polynomem na každém intervalu [xi, xi+1], i = 0, . . . , n − 1 a platí s(xi) = fi, i = 0, . . . , n. Podmínky spojitosti musíme uplatnit v uzlech, jelikož polynomy mají spojité derivace všech řádů, takže uvnitř intervalů není se spojitostí problém. Máme celkem n intervalů, na každém potřebujeme sestrojit polynom 3. stupně, který má čtyři neznámé koeficienty, dohromady tedy máme 4n neznámých koeficientů. Ovšem pro každý z n polynomů máme zadána dvě funkční hodnoty (jednu v každém krajním bodě intervalu [xi, xi+1]), tím se nám počet neznámých koeficientů redukuje na polovinu a současně zajistíme spojitost kubického splajnu. Ve vnitřních bodech (x1, . . . , xn−1) dále máme podmínky na spojitost první a druhé derivace, čímž se nám počet neznámých koeficientů zredukuje ne 2. Další podmínky už ale nemáme, takže kubický splan není funkčními hodnotami jednoznačně zadán a je potřeba dvě podmínky přidat,aby bylo možné jej sestrojit. Nejčastěji se zadávají hodnoty první derivace v krajních bodech x0 a sn, pak hovoříme o úplném kubickém splajnu, nebo hodnoty druhé derivace v těchto bodech. Pokud mají být tyto hodnoty nulové, splajn se nazývá 26 přirozený kubický splan. Nebudeme tady uvádět přesnou konstrukci kubických splajnů, zájemce najde podrobné informace třeba v [2]. Ukážeme si, jak zkonstruovat splajn za pomocí software. V Matlabu se pro konstrukci splajnu dá použít základní funkce spline. V nejjednodušší verzi počítá koeficienty polynomů, jako výsledek pak vrací složitější strukturu, která obsahuje víc informací. Pro data z předchozího obrázku dostaneme: >> x=[1 2 3 4 6 7]; >> f=[3 1 0 2 1 1]; >> plot(x,f,’-*’) >> print -depsc splajn1.eps >> s=spline(x,f) s = form: ’pp’ breaks: [1 2 3 4 6 7] coefs: [5x4 double] pieces: 5 order: 4 dim: 1 Označení ’pp’ ukazuje, že se jedná o po částech polynomiální funkci, Dále následují uzly, tedy hraniční body pro jednotlivé části polynomu. Pak můžeme zjisti koeficienty na jednotlivých intervalech, počet intervalů a stupeň polynomu (ve skutečnosti spíše počet koeficientů na každém intervalu). Význam poslední položky by měl jasný. Pokud bychom tedy chtěli znát koeficienty splajnu na jednotlivých intervalech stačí zadat >> s.coefs ans = 0.6978 -1.5935 -1.1044 3.0000 0.6978 0.5000 -2.1978 1.0000 -1.4891 2.5935 0.8956 0 0.4081 -1.8738 1.6153 2.0000 0.4081 0.5748 -0.9829 1.0000 Nás budou ovšem spíš zajímat hodnoty splanu, abychom si jej mohli případně nakreslit. K tomu lze použít funkci ppval: 27 >> xx=1:0.01:7; >> ss=ppval(s,xx); >> plot(xx,ss,x,f,’r*’) 1 2 3 4 5 6 7 −0.5 0 0.5 1 1.5 2 2.5 3 V případě, že nás nezajímají koeficienty na jednotlivých intervalech, ale jen výsledné funkční hodnoty, stačí zadat >> ss=spline(x,f,xx); Zvídavého čtenáře teď možná napadne, že jsme při konstrukci zadávali jen funkční hodnoty a žádné okrajové podmínky. Právem se tedy může ptát, jaké okrajové podmínky byly použity. Tuto otázku jistě uspokojivě zodpoví dokumentace Matlabu. V Sage se splajn definuje pro body v rovině, jimiž má procházet, takže pro stejná data jako výše můžeme použít následující postup: sage: values = [[1,3],[2,1],[3,0],[4,2],[6,1],[7,1]] sage: S = spline(values) sage: S [[1, 3], [2, 1], [3, 0], [4, 2], [6, 1], [7, 1]] Vidíme, že samotný obsah proměnné S nám moc informací nedá. Příkazem 28 sage: type(S) sage.gsl.interpolation.Spline aspoň zjistíme, o jaký datový typ se jedná, ale koeficienty na jednotlivých intervalech nezjistíme. A pokud vím, u této základní funkce se tyto koeficienty přímo zjistit nedají. Pokud by je člověk opravdu potřeboval, musí si stáhnout další knihovny, které obsahují poněkud více sofistikované funkce pro práci se splajny. Pro naše účely ovšem základní funkce zcela postačuje, například hodnotu polynomu v bodě zjistíme snadno: sage: S(1) 3.0 sage: S(1.5) 1.9917763157894737 sage: S(sqrt(2)) 2.1640477491461865 Nechat si splajn vykreslit (případně i s uzly a příslušnými hodnotami) taky není problém: plot(S,(1,7))+list_plot(values,color=’red’) 1 2 3 4 5 6 7 0 0.5 1 1.5 2 2.5 3 29 V Sage je velmi jednoduché změnit hodnotu v některém uzlu: sage: S[0] [1, 3] sage: S[0]=[1,2] sage: S [[1, 2], [2, 1], [3, 0], [4, 2], [6, 1], [7, 1]] sage: plot(S,(1,7))+list_plot(values,color=’red’) 1 2 3 4 5 6 7 0 0.5 1 1.5 2 2.5 3 První hodnota samozřejmě neleží na splajnu, jelikož jsem ji nepředefinovali. Splajn ale prochází danými body. Jednoduše se dá taky bod ubrat nebo naopak přidat: sage: del S[1] sage: S [[1, 2], [3, 0], [4, 2], [6, 1], [7, 1]] sage: S.append([5,2]) sage: S [[1, 2], [3, 0], [4, 2], [6, 1], [7, 1], [5, 2]] Když si necháme splajn zobrazit, vidíme, že nazáleží na pořadí uzlů. 30 1 2 3 4 5 6 7 0 0.5 1 1.5 2 1.5 Bersteinovy polynomy a Bézierovy křivky Kropáčku, co víte o Bersteinovi? West Side Story je fakt skvělý muzikál. Jako film dostal deset oskarů. Aha, tak příště zkuste do vyhledávače přidat křestní jméno Sergej. Sergej Bernstein byl ruský matematik, který se během svého dlouhého života dosáhl významných výsledků v několika matematických odvětvích. Jeho jméno by se mělo vyslovovat rusky, tedy „bernštejn a někdy se takto i v češtině píše. Často se ovšem setkáme s výslovností německou („bernštajn ) nebo s anglickou („bernstain ). Nás bude zajímat jeho příspěvek k teorií aproximací. Bernsteinovy polynomy totiž mají tu důležitou vlastnost, že s libovolnou přesností aproximují spojitou funkci na intervalu [0, 1] Jednoduchou substitucí pak jimi můžeme aproximovat spojitou funkci na libovolném uzavřeném intervalu. Jejich konstrukce je navíc velmi jednoduchá a tedy i snadno naprogramovatelná. Nechť tedy n je pevně dané přirozené číslo a k je přirozená číslo nebo nula, k ≤ n. Bersteinův bázový polynom bn,k definujeme vztahem bn,k(x) = n k xk (1 − x)n−k . 31 Dále nechť f je funkce definovaná na intervalu [0, 1], položme fk = f(k n ) pro k = 0, . . . , n. Bernsteinův polynom stupně n funkce f definujeme výrazem Bn,f (x) = n k=0 fk · bn,k(x). (1.5) Snadno se ověří, že n k=0 fk · bn,k(x) = 1, odkud pak bezprostředně plyne, že pro funkci f konstantní na [0, 1] platí f(x) = Bn,f (x) pro každé x ∈ [0, 1]. Totéž platí i v případě, že f je polynomem prvního stupně, jak se dá též poměrně snadno spočítat. Pro polynomy vyšších stupňů to už ale neplatí, jak dokládá následující obrázek, kde čárkovaně je zobrazena funkce x2 , křivky, které se k ní blíží jsou postupně Bernsteinovy polynomy druhého, třetího a čtvrtého stupně. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Tento obrázek byl získán pomocí matlabovského programu pro výpočet Bersteinova polynomu, přesněji řečeno jeho koeficientů. Jako vedlejší produkt dostáváme i Bernsteinovy bázové polynomy uspořádané v matici. Při jejich konstrukci se využívá toho, že daný bázový polynom má kořeny 0 a 1 násobností k a n − k. function [BP,B] = bern_pol(fk) % function BP = bern_pol(fk) 32 % Bernsteinuv polynom pro hodnoty fk % B - matice Bern. bazovych polynomu n=length(fk)-1; B=[]; % matice bazovych Bernsteinovych polynomu for k=0:n bk=nchoosek(n,k)*poly( [zeros(1,k),ones(1,n-k)]); bk=bk*(-1)^(n-k); B=[B;bk]; end BP=fk*B; end Uvedený obrázek pak dostaneme například pomocí následujícího postupu: >> f=inline(’x.^2’); >> n=2; >> k=0:n; >> xk=k/n; >> f2=f(xk); >> B2=bern_pol(f2); >> n=3; >> k=0:n; >> xk=k/n; >> f3=f(xk); >> B3=bern_pol(f3); >> n=4; >> k=0:n; >> xk=k/n; >> f4=f(xk); >> B4=bern_pol(f4); >> xx=0:0.01:1; >> Bx2=polyval(B2,xx); >> Bx3=polyval(B3,xx); >> Bx4=polyval(B4,xx); >> plot(xx,[Bx2;Bx3;Bx4],xx,xx.^2,’--’) Nejdůležitější vlastností Bernsteinových polynomů, jak už bylo nazna- 33 čeno, je, že pro funkcí f spojitou na [0.1] konverguje posloupnost Bn,f stejnoměrně k f Kromě teoretického významu se Bernsteinovy polynomy používají i prakticky, totiž ke konstrukci Bězierových křivek. „Legenda praví, že tyto křivky objevili nezávisle na sobě dva pracovníci konstrukčních kancelářích dvou renomovaných francouzských automobilek. První objevitel ovšem přišel o prvenství tím, že kvůli konkurenčnímu boji jeho objev podléhal utajení. V dnešní době jsou Bézierovy křivky běžnou součástí téměř každého kreslícího software, často bez uvedení, o jaký typ křivek se jedná. Nejčastěji můžeme narazit na kubické Bézierovy křivky, jejichž konstrukci si zde ukážeme. Zobecnění pro vyšší stupně je zcela přímé a čtenář by je jistě bez problémů zvládl. Bézierovy křivky jsou zadány parametricky, konstrukce je navíc v podstatě totožná v rovině i v trojrozměrném prostoru. Mějme tedy 4 body v rovině, označme je P0,P1,P2,P3, Pi = [xi, yi], tyto body se nazývají kontrolní nebo řídící. Bézierova křivka je množina všech bodů Q = [x, y) pro něž platí x = x(t) = x0 · b3,0(t) + x1 · b3,1(t) + x2 · b3,2(t) + x3 · b3,3(t) = x0(1 − t)3 + 3x1t(1 − t)2 + 3x2t2 (1 − t) + x3t3 y = y(t) = y0 · b3,0(t) + y1 · b3,1(t) + y2 · b3,2(t) + y3 · b3,3(t) = y0(1 − t)3 + 3y1t(1 − t)2 + 3y2t2 (1 − t) + y3t3 , kde proměnná t probíhá interval [0, 1] Pokud bychom uvažovali body v prostoru, t.j. Pi = [xi, yi, zi], přibyla by nám rovnice pro třetí souřadnici bodu Q: z = z(t) = z0 · b3,0(t) + z1 · b3,1(t) + z2 · b3,2(t) + z3 · b3,3(t) = z0(1 − t)3 + 3z1t(1 − t)2 + 3z2t2 (1 − t) + z3t3 Souřadnice bodu Q na Bézierově křivce jsou tedy hodnoty Bersteinova polynomu pro příslušné souřadnice kontrolních bodů. Stručně se body na Bézierově křivce dají vyjádřit jako Q = Q(t) = P0(1 − t)3 + 3P1t(1 − t)2 + 3P2t2 (1 − t) + P3t3 Z vyjádření navíc plyne, že pro t = 0 dostaneme výchozí kontrolní bod P0, pro t = 1 koncový kontrolní bod P3. Sestrojit Bézierovu křivku v Matlabu je velmi jednoduché, není ani potřeba vytvářet pro tento účel zvláštní funkci: 34 >> P0=[0;1]; P1=[2;3];P2=[4;2];P3=[4;0]; >> t=0:0.01:1; >> Q=P0*((1-t).^3)+3*P1*(t.*(1-t).^2)+... 3*P2*(t.^2.*(1-t))+P3*(t.^3); +>> plot(Q(1,:),Q(2,:)) Pokud chceme zobrazit i kontrolní body, nejjednodušší je poskládat je do matice. kterou zobrazíme po řádcích: >> hold on >> PP=[P0,P1,P2,P3]; >> plot(PP(1,:),PP(2,:),’r*’) 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 Pro konstrukci Bézierových křivek se dá použít také algoritmus de Casteljau pojmenovaný po jednou z objevitelů Bézierových křivek. Spočívá v rozdělení spojnic mezi kontrolními body v daném poměru. Tím získáme tři body, jejichž spojnice opět rozdělíme ve stejném poměru, Stejně tak rozdělíme i spojnici takto získaných dvou bodů a tím získáme bod, který leží na Bézierově křivce. Následující obrázek ukazuje tento proces pro dělící poměr 1:1 a 1:2, to znamená, že všechny spojnice dělíme na poloviny, respektive na třetinu a dvě třetiny. 35 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 To, že uvedený postup skutečně dává body na Bézierově křivce, snadno dokážeme následujícím výpočtem. Dělící poměr pro rozdělení spojnice dvou bodů je určen hodnotou α ∈ (0, 1), např. bod C na spojnici bodů A a B se dá vyjádřit jako C = α · A + (1 − α) · B, přičemž výsledný bod je blízko bodu A pro α blízké jedné. Délky údeček AC a CB jsou pak v poměru (1 − t) : t. Položme tedy R0 = α · P0 + (1 − α) · P1, R1 = α · P1 + (1 − α) · P2, R2 = α · P2 + (1 − α) · P3 čímž získáme první trojici bodů. Pokračujme dále: S0 = α · R0 + (1 − α) · R1, S1 = α · R1 + (1 − α) · R2 a konečně T = α · S0 + (1 − α) · S1. Zpětným dosazováním dostáváme T = α · S0 + (1 − α) · S1 = α · (α · R0 + (1 − α) · R1) + (1 − α) · (α · R1 + (1 − α) · R2) 36 = α2 · R0 + 2α(1 − α) · R1 + (1 − α)2 · R2 = α2 · (α · P0 + (1 − α) · P1) + 2α(1 − α) · (α · P1 + (1 − α) · P2) + +(1 − α)2 · (α · P2 + (1 − α) · P3) = α3 · P0 + 3α2 (1 − α) · P1 + 3α(1 − α)2 · P2 + (1 − α)3 · P3. Vidíme, že pro t = 1 − α jsem dostali přesně vyjádření bodu na Bézierově křivce. 37 Kapitola 2 Derivace Pane profesore, chápu to správně, že když si jako funkci vyjádříme dráhu tělesa v v závislosti na čase, tak její derivace podle času je rychlost tělesa? Ano, Kropáčku, to je běžný matematický model, který ve většině případů dává uspokojivé výsledky. A jak zjistíme rychlost času, pane profesore? To je velmi dobrá otázka, Kropáčku. Jak jste k ní dospěl? No, některé přednášky, tím samozřejmě nemyslím vaše, se vlečou tak, že hodina trvá celé dopoledne, zatímco při posledním testu jsem si sotva stačil opsat zadání a hodina byla pryč. 2.1 Limity Ještě před samotnými derivacemi se aspoň krátce zmíníme o limitách. V Sage to není problém, stačí zadat například sage: n=var(’n’) sage: limit((1+1/n)^n,n=oo) e sage: limit((1-1/n)^n,n=oo) e^(-1) sage: limit((1+10/n)^n,n=oo) e^10 sage: x=var(’x’) 38 sage: limit((1+x/n)^n,n=oo) e^x sage: limit(sin(x)/x,x=0) 1 V Matlabu symbolické limity nenajdeme, nicméně některé limity pro n → ∞ se dají spočítat nebo odhadnout tak, že daný výraz se pokoušíme spočítat pro hodně velká čísla. Má to ale svoje úskalí. Pokud bychom se pomocí výše uvedené limity pokoušeli určit Eulerovo číslo s přesností na 6 desetinných míst, dostáváme: >> format long >> n=1; >> e0=(1+1/n)^n e0 = 2 >> n=10; >> e1=(1+1/n)^n e1 = 2.593742460100002 >> while abs(e1-e0)>0.000001, e0=e1; n=10*n;... e1=(1+1/n)^n, end e1 = 2.704813829421528 e1 = 2.716923932235594 e1 = 2.718145926824926 e1 = 2.718268237192297 e1 = 2.718280469095753 e1 = 2.718281694132082 e1 = 2.718281798347358 Pokud bychom ale chtěli tímto způsobem určit Eulerovo číslo s přesností na 15 desetinných čísel, což je zhruba přesnost, s jakou jsou běžně v počítači 39 uložena čísla řádově v jednotkách, vyjde >> while abs(e1-e0)>0.000000000000001, e0=e1; n=10*n;... e1=(1+1/n)^n, end e1 = 2.704813829421528 e1 = 2.716923932235594 e1 = 2.718145926824926 e1 = 2.718268237192297 e1 = 2.718280469095753 e1 = 2.718281694132082 e1 = 2.718281798347358 e1 = 2.718282052011560 e1 = 2.718282053234788 e1 = 2.718282053357110 e1 = 2.718523496037238 e1 = 2.716110034086901 e1 = 2.716110034087023 e1 = 3.035035206549262 e1 = 1 e1 = 1 40 Pro velká n je totiž hodnota 1/n menší než přesnost, s jakou je uložena jednička, takže při sčítání 1 + 1/n dostáváme výsledek 1. Pokud bychom chtěli přesto s pomocí Matlabu Eulerovo číslo určit jako limitu nějakého nekonečného procesu, můžeme využít vztah e = ∞ n=0 1 1! . Součet nekonečné řady je limita částečných součtů, které můžeme snadno vyjádřit: >> n=0;a=1;s=1; >> while a>0.000000000000001, n=n+1;a=a/n;s=s+a; end >> s s = 2.718281828459046 >> exp(1) ans = 2.718281828459046 >> n n = 18 Ve výpisu vidíme, že dosažená hodnota je v rámci přesnosti stejná jako hodnota, se kterou Matlab pracuje. Na určení stačilo sečíst 18 členů řady, protože 1/18! je řádově 10−16 , tedy další členy by se ve výsledku již neprojevily. 2.2 Symbolické derivování Symbolické derivování je poměrně snadné v Matlabu i v Sage. V Matlabu si nejdříve definujeme symbolický objekt a pak s ním můžeme provádět symbolické operace, tedy i derivování: >> f1=sym(’sin(x)*cos(x)’) f1 = cos(x)*sin(x) >> diff(f1) ans = 41 cos(x)^2 - sin(x)^2 >> f2=sym(’log(t^2)’); >> diff(f2) ans = 2/t Vidíme, že pokud zadaná funkce obsahuje jen jednu proměnnou, derivuje se automaticky podle ní a výsledek je současně upraven. Pokud máme funkci více proměnných, je potřeba si zvolit, podle které chceme derivovat, pokud to není zrovna x. >> f3=sym(’cos(x^2)*sin(y)’); >> diff(f3) ans = -2*x*sin(x^2)*sin(y) >> diff(f3,’y’) ans = cos(x^2)*cos(y) Je taky možné počítat vyšší derivace: >> diff(f2,3) ans = 4/t^3 >> diff(f3,’y’,2) ans = -cos(x^2)*sin(y) Derivování v Sage je podobně jednoduché, jen nejdříve musíme definovat proměnné. Máme dokonce na výběr, který příkaz pro derivování použijeme a samozřejmostí je počítání vyšších derivací. sage: x=var(’x’) sage: t=var(’t’) sage: diff(exp(x)*cos(x),x) -e^x*sin(x) + e^x*cos(x) sage: derivative(exp(x)*cos(x),x) -e^x*sin(x) + e^x*cos(x) 42 sage: diff(log(t)*exp(x),t) e^x/t sage: diff(log(t)*exp(x),x) e^x*log(t) sage: diff(log(t)*exp(x),t,3) 2*e^x/t^3 2.3 Taylorův rozvoj Určit Taylorův rozvoj funkce, když dokážeme počítat derivace, je už snadné. Spočítali bychom derivace funkce v daném bodě, určili tak příslušné koeficienty, věřím, že by s s tímto problémem čtenář bez problémů poradil. V Matlabu ale existuje hezký prográmek taylortool, ve kterém si můžete zadat, po jakou funkci chcete Taylorův rozvoj spočítat, ve kterém bodě a do jakého stupně a jako výsledek uvidíte něco jako na dalším obrázku. −6 −4 −2 0 2 4 6 −10 −5 0 5 10 Taylor Series Approximation TN (x) =− x7 /720 + x5 /24 − x3 /2 + x Jak by se dalo čekat. určit Taylorův rozvoj funkce v Sage je jednoduché. Stačí definovat proměnnou, funkci a pak použít funkci taylor: sage: x=var(’x’) sage: f1=sqrt(x+1) sage: f1.taylor(x,0,6) 43 -21/1024*x^6 + 7/256*x^5 - 5/128*x^4 + 1/16*x^3 - 1/8*x^2 + 1/2*x + 1 sage: f2=sqrt(x) sage: f2.taylor(x,1,6) -21/1024*(x - 1)^6 + 7/256*(x - 1)^5 - 5/128*(x - 1)^4 + 1/16*(x - 1)^3 - 1/8*(x - 1)^2 + 1/2*x + 1/2 Jde to i bez definování funkce: sage: taylor(sin(x)*cos(x),x,pi,6) -pi - 2/15*(pi - x)^5 + 2/3*(pi - x)^3 + x Možná by stálo za to trochu prozkoumat, podle čeho se určuje pořadí členů na výstupu. 2.4 Numerické derivování Může se stát, že známe hodnoty funkce jen v diskrétních uzlech, nicméně potřebujeme spočítat alespoň přibližně hodnotu její derivace v nějakém bodě. Zpravidla se jedná o některý z uzlů, ale nemusí to tak být vždy. V takové situaci máme dvě možnosti: sestrojíme interpolační polynom a ten zderivujeme, nebo si vypomůžeme Taylorovým rozvojem. Oba přístupy dávají stejné výsledky, u Taylorova rozvoje dostaneme navíc odhad chyby, výsledky dosažené derivací interpolačního polynomu zase dostaneme bez velkého přemýšlení. Začneme tedy s interpolačním polynomem. Zde se hodí Lagrangeův tvar Pn(x) = n i=0 fili(x), kde li jsou Lagrangeovy fundamentální polynomy. Derivováním této formule dostaneme f (x) ≈ Pn(x) = n i=0 fili(x). Podrobnosti o přesnosti a dalších vlastnostech této formule se čtenář může dočíst v [1], my se zde blíže podíváme jen na některé speciální případy. Odvodíme formuli pro přibližný výpočet derivace v prostředním ze tří ekvidistantních uzlů. Kvůli symetrii formulí si je označíme trochu jinak než obvykle: x−1, x0, x1, xi = x0 + ih, i = ±1. 44 Nejprve sestrojíme Lagrangeův interpolační polynom: xi x−1 x0 x1 fi f−1 f0 f1 P2(x) = f−1 (x − x0)(x − x1) (x−1 − x0)(x−1 − x1) f−1 + f0 (x − x−1)(x − x1) (x0 − x−1)(x0 − x1) + + f1 (x − x−1)(x − x0) (x1 − x−1)(x1 − x0) , Výraz zderivujeme P2(x) = f−1 2x − x0 − x1 2h2 − f0 2x − x−1 − x1 h2 + f1 2x − x−1 − x0 2h2 , pro xi = x0 + ih, i = ±1. Pokud dosadíme x = x0 dostaneme f (x0) ≈ P2(x0) = 1 2h (f1 − f−1). (2.1) Všimněme si ještě geometrického významu této formule. Výraz (f1 −f−1)/2h je směrnice sečny, která je určena body (x−1, f−1) a (x1, f1). Podobně můžeme odvodit i formule pro výpočet derivace v dalších uzlových bodech x−1, x1: f (x−1) ≈ 1 2h (−3f−1 + 4f0 − f1) (2.2) f (x1) ≈ 1 2h (f−1 − 4f0 + 3f1) (2.3) Tyto formule se nazývají tříbodové. Podívejme se ještě, jak se dají stejné formule odvodit z Taylorova rozvoje. Pro uzly x±1 = x0 ± h dostáváme f(x1) = f(x0) + f (x0)h + f (x0) 2 h2 + f (x0) 6 h3 + f(4) (x0) 24 h4 + · · · f(x−1) = f(x0) − f (x0)h + f (x0) 2 h2 − f (x0) 6 h3 + f(4) (x0) 24 h4 − · · · Odečtením obou rovností a vydělením 2h dostáváme f (x0) = 1 2h (f1 − f−1) − f (x0) 3 h2 − · · · , (2.4) 45 takže máme stejnou formuli jako v předchozím případě, navíc vidíme, že chyba je řádově h2 , takže pokud například hodnotu h zmenšíme na polovinu, chyba klesne přibližně na čtvrtinu. Pokud bychom pro přibližné vyjádření derivace v bodě x0 použili jen jeden Taylorův rozvoj, dostaneme f (x0) = 1 h (f1 − f0) − f (x0) 2 h − · · · , tedy výraz s chybou o řád horší než vztah předchozí. Uvedené Taylorovy rozvoje ale můžeme i sečíst, v tom případě dostaneme vztah pro přibližný výpočet druhé derivace: f (x0) = 1 h2 (f−1 − 2f0 + f1) − f(4) (x0) 12 h2 − · · · . (2.5) Vidíme, že chyba uvedeného vztahu je řádově opět rovna h2 . 46 Kapitola 3 Integrály Pane profesore, nevíte prosím, na jakých matematických principech je postavena důchodová reforma? To netuším, Kropáčku, ale domnívám se, že to nebude nic složitého. Proč se na to ptáte? No, četl jsem včera na jednom internetovém serveru, že „Zřízení sociální pojišťovny je integrální součástí důchodové reformy. Ale není mi jasné, jaké funkce se vlastně integruje. 3.1 Symbolické integrování Co se týče počítání symbolického počítání neurčitých i určitých integrálů, zvládají to oba námi používané softwarové prostředky bez velkých problémů. Třeba v Matlabu pro výpočet neurčitého nebo určitého integrálu používáme funkci int, rozdíl je jenom v zadaných parametrech: >> f=sym(’sin(x)*cos(2*x)’); >> int(f) ans = cos(x) - (2*cos(x)^3)/3 >> int(f,0,pi) ans = -2/3 47 Je možné i integrovat podle jedné z více proměnných: >> g=sym(’exp(-t)*sin(2*x*t)’); >> int(g,’t’) ans = -(sin(2*t*x) + 2*x*cos(2*t*x))/(exp(t)*(4*x^2 + 1)) >> int(g,’x’,0,pi/2) ans = -(cos(pi*t) - 1)/(2*t*exp(t)) Dají se dokonce vyjádřit integrály z funkcí, jejichž primitivní funkce sice existují, ale nedají se vyjádřit v rozumném konečném tvaru: >> G=sym(’exp(-x^2)’); >> int(G) ans = (pi^(1/2)*erf(x))/2 >> I=int(G,0,1) I = (pi^(1/2)*erf(1))/2 >> eval(I) ans = 0.7468 >> int(G,0,inf) ans = pi^(1/2)/2 Výrazem erf v Matlabu získáme primitivní funkci k e−x2 , jak zjistíme z nápovědy, a Matlab umí spočítat hodnoty této funkce: >> help erf erf Error function. Y = erf(X) is the error function for each element of X. X must be real. The error function is defined as: erf(x) = 2/sqrt(pi) * integral from 0 to x of exp(-t^2) dt. 48 See also erfc, erfcx, erfinv, erfcinv. Overloaded methods: sym/erf Reference page in Help browser doc erf >> erf(1) ans = 0.8427 V Sage je to velmi podobné: sage: x=var(’x’) sage: integral(sin(x)*x,x) -x*cos(x) + sin(x) sage: integral(sin(x)*x,x,0,pi) pi sage: integral(sin(x)*x,x,0,pi/2) 1 sage: t=var(’t’) sage: integral(sin(x),t) t*sin(x) sage: integral(1/sqrt(x),x,0,1) 2 sage: integral(1/x^2,x,0,1) Traceback (click to the left of this block for traceback) ... ValueError: Integral is divergent. I primitivní funkce k e−x2 se vyjadřuje podobně: sage: integral(exp(-x^2),x) 1/2*sqrt(pi)*erf(x) sage: integral(exp(-x^2),x,0,oo) 1/2*sqrt(pi) sage: I0=integral(exp(-x^2),x,0,1);I0 49 1/2*sqrt(pi)*erf(1) sage: I0.N() 0.746824132812427 3.2 Numerické integrování Při numerickém integrování se podobně jako při derivování nesnažíme určit primitivní funkci z funkčních hodnot v zadaných diskrétních uzlech, ale snažíme se na základě těchto dat určit přibližně určitý integrál přes nějaký interval. Základem získaných formulí jsou opět interpolační polynomy. Při integrování Lagrangeova tvaru dostáváme: b a f(x)dx ≈ b a Pn(x)dx = n i=0 fi b a li(x)dx = n i=0 fiAi pro Ai = b a li(x)dx. Podobný typ výrazů dostáváme také při odvozování Riemannova integrálu, kde aproximaci určitého integrálu dostáváme jako součet funkčních hodnot násobených koeficienty zahrnujícími délky subintervalů při dělení původního intervalu na menší dílky. Těmto výrazům, tedy b a f(x)dx ≈ n i=0 fiAi říkáme kvadraturní formule. Kromě odhadu chyby nás u kvadraturních formulí zajímá také její stupeň přesnosti, který udává, pro polynomy jakých stupňů je tato formule přesná. Přesněji, pokud je stupeň přesnosti formule roven n pak je tato formule je přesná pro polynomy do stupně n včetně. Dá se snadno ukázat (viz [1]), že pro n + 1 uzlů je stupeň přesnosti roven maximálně 2n + 1. Na druhou stranu je zřejmé, že pokud kvadraturní formuli získáme postupem ukázaným výše, tedy integrací interpolačního polynomu, musí být stupeň přesnosti roven alespoň n. Konstrukce kvadraturních formulí je možné ještě poněkud zobecnit tím, že do integrálu přidáme ještě tzv. vahovou funkci. Do ní můžeme zahrnout například singularity nebo společnou část pro nějakou třídu funkcí. V tom 50 případě kvadraturní formule nahrazují integrál n i=0 fiAi ≈ b a w(x(·f(x)dx. Koeficienty kvadraturní formule pak získáme ze vztahu Ai = b a w(x) · li(x)dx. Výhodou takového postupu je, že vahová funkce není zahrnuta ve funkčních hodnotách funkce f. V tomto textu se ale budeme zabývat pouze jednoduchou situací s vahovou funkcí rovnou jedné. Odvodíme si některé nejpoužívanější formule. Situací s jedním uzlem se zabývat nebudeme, tak je až příliš jednoduchá. Pro dva uzly většinou máme a = x0 < x1 = b. V tomto případě má interpolační polynom tvar například P1(x) = f(a) + f(b) − f(a) b − a (x − a). Jeho integrováním přes interval [a, b] vyjde kvadraturní formule b a f(x)dx ≈ f(a) + f(b) 2 (b − a), což je v podstatě plocha lichoběžníka (postaveného na bok),jehož strany mají délky rovny funkčním hodnotám funkce f v krajních bodech intervalu a jehož výška je rovna b − a. Proto se taky této formuli říká lichoběžníkové pravidlo. 51 a b 0 V případě, že pudeme funkci aproximovat polynomem druhého stupně, tedy parabolou, dostaneme pro uzly x0 = a, x1 = a+b 2 , x2 = b formuli b a f(x)dx ≈ b − a 6 f(a) + 4f( a + b 2 ) + f(b) , která se nazývá Simpsonovo pravidlo. a (a+b)/2 b 0 52 Z těchto dvou formulí se odvozují patrně nejpoužívanější kvadraturní formule. Jejich myšlenka je jednoduchá: pokud máme více uzlů, pak nám rozdělují interval [a, b] na více subintervalů, na každém tomto subintervalu použijeme lichoběžníkové pravidlo nebo Simpsonovo pravidlo a výsledek pak sečteme. Tím dostaneme složené lichoběžníkové nebo složené Simponovo pravidlo. U složeného lichoběžníkovécho pravidla v podstatě počítáme integrál z aproximace funkce lomenou čárou, respektive lineárním splajnem. U Složeného Simpsonova pravidla zase vyžadujeme, aby celkový počet uzlů byl lichý, na což by jistě čtenář přišel sám. pro ekvidistantní uzly x0, . . . , xn, kde vzdálenost sousedních dvou je rovna h dostáváme složené lichoběžníkové pravidlo ve tvaru b a f(x)dx ≈ h 2 (f0 + 2f1 + 2f2 + · · · + 2fn−1 + fn) , kde fi = f(xi). −1 −0.5 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x y Složené Simponovo pravidlo pak je b a f(x)dx ≈ h 3 (f0 + 4f1 + 2f2 + 4f3 + 2f4 + · · · + 2fn−2 + 4fn−1 + fn) . 53 Kapitola 4 Řady Kryšpín Kropáček: Řada Při sčítání řady dočkal jsem se zrady ze strany své spolužačky Novákové Lady. Nevím si s tím rady, skoro šilhám hlady, nejvíc mě však rozptylují její krásné vnady. Řad už jsme částečné dotkli, když jsme se zmínili o Taylorově rozvoji v souvislosti s derivacemi. Teď si zkusíme ukázat, jak se dají počítat některé součty řad a nesmíme zapomenout taktéž na řady Fourierovy. 4.1 Symbolické součty řad Jak se dá čekat, v Matlabu je se symbolickým sčítáním řad potíž. Můžeme zde maximálně odhadnout číselný součet řady, jestliže připočítáváme k součtu členy, které jsou už z numerického hlediska zanedbatelné, jak jsme to ukázali při výpočtu Eulerova čísla. Proto se v dalším výkladu soustředíme na Sage. Nejprve si zkusíme sečíst jednoduchou geometrickou řadu: var(’x’,’k’,’n’) sum(x^k,k,0,oo) Traceback (click to the left of this block for traceback) ... Is abs(x)-1 positive, negative, or zero? 54 Získali jsme pouze chybové hlášení. Chybí totiž bližší informace o proměnné x a víme, že uvedená řada má součet jen v případě, že |x| < 1. Proto tento požadavek zadáme do Sage: sage: assume(abs(x)<1) sage: sum(x^k,k,0,oo) -1/(x - 1) Vidíme, že tohle funguje a můžeme zkoušet další součty: sage: sum(k*x^k,k,1,oo) x/(x^2 - 2*x + 1) sage: sum(1/k^4, k, 1, oo) 1/90*pi^4 sage: sum(x^k/factorial(k), k, 0, oo) e^x Zatím to vypadá všechno bezproblémově. V některých situacích si ovšem ani Sage neporadí. Ukážeme si to na poměrně jednoduchém příkladu. Nejprve si spočítáme Taylorův rozvoj funkce √ 1 − x. sage: f=sqrt(1-x) sage: f.taylor(x,0,6) -21/1024*x^6 - 7/256*x^5 - 5/128*x^4 - 1/16*x^3 - 1/8*x^2 - 1/2*x + 1 Čísla, která se zde objevují, jsou binomické koeficienty, které jsou pro reálný argument a definované vztahem a k = a(a − 1) · · · (a − k + 1) k! V našem případě jsou to koeficienty pro a = 1/2 (až na znaménko), jak snadno ověříme: [binomial(1/2,k) for k in range(7)] [1, 1/2, -1/8, 1/16, -5/128, 7/256, -21/1024] Zkusíme tedy řadu zpětně sečíst. Problémem bude trochu první člen, ale ten by se měl projevit jen posunem výsledku o konstantu: 55 sum(x^k*(-1)^k*binomial(1/2,k),k,0,oo) Traceback (click to the left of this block for traceback) ... TypeError: Either m or x-m must be an integer Vidíme, že na tomto Sage ztroskotal. 4.2 Fourierovy řady Spočítat koeficienty Fourierovy řady na základě toho, co už známe, jistě nebude pro čtenáře problém. Například pokud bychom chtěli v Matlabu určit prvních pár členů Fourierovy řady pro funkci x2 , můžeme postupovat třeba takto: >> f=sym(’x^2’); >> a0=int(f,-pi,pi)/(2*pi) a0 = pi^2/3 >> a1=int(’cos(x)’*f,-pi,pi)/pi a1 = -4 >> a2=int(’cos(2*x)’*f,-pi,pi)/pi a2 = 1 >> a3=int(’cos(3*x)’*f,-pi,pi)/pi a3 = -4/9 Koeficienty u sinových členů jsou nulové, neboť x2 je sudá funkce. Sečteme tedy první čtyři členy řady a výsledek porovnáme s původní funkcí: >> x=-pi:0.01:pi; >> S=eval(a0+a1*cos(x)+a2*cos(2*x)+a3*cos(3*x)); >> plot(x,S,x,x.^2,’--’) Pro výpočet hodnot součtů v bodech daných vektorem x musíme použít funkci eval, protože koeficienty a0,. . . ,a3 jsou symbolické objekty. 56 −4 −3 −2 −1 0 1 2 3 4 −2 0 2 4 6 8 10 Pokud chceme určit Fourierovu řadu v Sage, je potřeba si funkci definovat jen na příslušném intervalu, tedy např. na [−π, π]: sage: x=var(’x’) sage: f=Piecewise([[(-pi,pi),x^2]]) Stejné koeficienty jako v Matlabu určíme pomocí příkazů sage: f.fourier_series_cosine_coefficient(0,pi) 2/3*pi^2 sage: f.fourier_series_cosine_coefficient(1,pi) -4 sage: f.fourier_series_cosine_coefficient(2,pi) 1 sage: f.fourier_series_cosine_coefficient(3,pi) -4/9 První vstupní parametr říká, který koeficient určujeme, druhý je polovina intervalu, na němž řadu počítáme. Částečný součet řady ale můžeme získat i přímo a porovnat ho s původní funkcí: sage: g=f.fourier_series_partial_sum(3,pi) sage: pl=plot([g,x^2],-pi,pi) 57 -3 -2 -1 1 2 3 2 4 6 8 10 58 Literatura [1] Horová, I.; Zelinka, J.: Numerické metody. Brno: Masarykova univerzita, druhé vydání, 2008, ISBN 978-80-210-3317-7. [2] Kobza, J.: Splajny. Univerzita Palackého, 1993, ISBN 9788070672655. URL http://books.google.cz/books?id=-OnXYgEACAAJ [3] Ralston, A.: Základy numerické matematiky. Praha: Academia, druhé vydání, 1978. [4] Schumaker, L.: Spline Functions: Basic Theory. Cambridge University Press, ISBN 9781139463430. URL http://books.google.cz/books?id=2uZLawUhXfgC 59