Univerzita Karlova v Praze Matematicko-fyzik´aln´ı fakulta
´ RSK ˇ ´ PRACE ´ BAKALA A
Ludˇek Gregor Kvadraturn´ı formule pro funkce s vysokou oscilac´ı Katedra numerick´e matematiky
Vedouc´ı bakal´aˇrsk´e pr´ace: Doc. RNDr. Josef Kofroˇ n, CSc. Studijn´ı program: Numerick´a matematika
Prohlaˇsuji, ˇze jsem svou bakal´aˇrskou pr´aci napsal samostatnˇe a v´yhradnˇe s pouˇzit´ım citovan´ych pramen˚ u. Souhlas´ım se zap˚ ujˇcov´an´ım pr´ace a jej´ım zveˇrejˇ nov´an´ım. V Praze dne
Ludˇek Gregor
2
Obsah 1 Pouˇ zit´ e znaˇ cen´ı
5
´ 2 Uvod
6
3 Metody numerick´ e integrace pro obecnou osciluj´ıc´ı funkci 3.1 Pouˇzit´ı vhodn´e v´ahov´e funkce . . . . . . . . . . . . . . . . . 3.2 Vyuˇzit´ı integrace per partes . . . . . . . . . . . . . . . . . . 3.3 Numerick´a integrace mezi nulov´ymi body . . . . . . . . . . . 3.4 Analytick´a integrace osciluj´ıc´ı ˇc´asti . . . . . . . . . . . . . . 3.5 Eulerova transformace . . . . . . . . . . . . . . . . . . . . .
8 8 9 11 13 13
4 Metody numerick´ e integrace speci´ aln´ıch typ˚ u funkc´ı 4.1 Priceova metoda . . . . . . . . . . . . . . . . . . . . . . . 4.2 Filonova metoda . . . . . . . . . . . . . . . . . . . . . . . 4.3 Pantisova metoda . . . . . . . . . . . . . . . . . . . . . . . 4.4 Filonova lichobˇeˇzn´ıkov´a metoda . . . . . . . . . . . . . . . 4.5 Filonova spline metoda . . . . . . . . . . . . . . . . . . . . 4.6 Vyuˇzit´ı ortogon´aln´ıch syst´em˚ u polynom˚ u . . . . . . . . . . 4.7 Pˇreveden´ı probl´emu na ˇreˇsen´ı obyˇcejn´e diferenci´aln´ı rovnice
16 16 18 20 21 22 23 25
. . . . . . .
5 Realizace vybran´ ych metod
27
Literatura
30
3
N´azev pr´ace: Kvadraturn´ı formule pro funkce s vysokou oscilac´ı Autor: Ludˇek Gregor Katedra: Katedra numerick´e matematiky Vedouc´ı bakal´aˇrsk´e pr´ace: Doc. RNDr. Josef Kofroˇ n, CSc. e-mail vedouc´ıho:
[email protected] Abstrakt: V pˇredloˇzen´e pr´aci studujeme metody aproximuj´ıc´ı hodnotu urˇcit´eho integr´alu funkc´ı s vysokou oscilac´ı. Vyuˇz´ıv´ame v praxi obvykl´eho tvaru zkouman´ych funkc´ı, vyskytuj´ıc´ı se napˇr´ıklad u Fourierovov´ych ˇrad - kde je integrovan´a funkce souˇcinem rychle osciluj´ıc´ı a obecnˇe neosciluj´ıc´ı funkce. Pˇrirozenou cestou je aproximace neosciluj´ıc´ı funkce tak, abychom dostali souˇcin funkc´ı, jeˇz je snadno analyticky integrovateln´ y. Typickou volbou jsou funkce, kter´e jsou spojit´e, a po ˇc´astech polynomi´aln´ı. D´ale je moˇzn´e u urˇcit´e tˇr´ıdy probl´em˚ u vyuˇz´ıt integrace per partes, a dostat tak aproximaci zaloˇzenou v´yhradnˇe na bodov´ych hodnot´ach odpov´ıdaj´ıc´ıch funkc´ı. Jsou pops´any metody pˇripisovan´e Filonovi, Pantisovi a Priceovi. Velmi zhruba je pops´ana myˇslenka Eulerovy transformace slouˇz´ıc´ı ke sˇc´ıt´an´ı ˇrad. Informativnˇe je zm´ınˇena metoda popisuj´ıc´ı pˇreveden´ı studovan´eho probl´emu na ˇreˇsen´ı odpov´ıdaj´ıc´ı diferenci´aln´ı rovnice. Kl´ıˇcov´a slova: osciluj´ıc´ı, integrace, Filon, kvadratura.
Title: Quadrature formulae for rapidly oscillating functions Author: Ludˇek Gregor Department: Department of Numerical Mathematics Supervisor: Doc. RNDr. Josef Kofroˇ n, CSc. Supervisor’s e-mail address:
[email protected] Abstract: In the present work we study aproximation methods for values of integrals with strongly oscillating integrand. It’s typical that the integrand can be split into two factors such that one is a rapidly oscillating and the other is a smooth function - typical example are Fourier series. A usual way is aproximation of a smooth function to get product, which is simply integrable. A common choice are continuous and piecewice polynomial functions. It is possible to use per-partes for aproximations, which need only several pointvalues. We describe Filon’s, Pantis’ and Price’s methods. Informatively we describe idea of Euler’s transformation for summation of serries and method, which converts our problem into solution of ODE. Keywords: oscillation, integration, Filon, quadrature.
4
Kapitola 1 Pouˇ zit´ e znaˇ cen´ı N, respektive R, respektive C znaˇc´ı mnoˇzinu pˇrirozen´ych, respektive re´aln´ych, respektive komplexn´ıch ˇc´ısel. R znaˇc´ı integr´al v Lebesgueovˇe smyslu. Pojmem funkce m´ın´ım zobrazen´ı z R do R, nebo do C. Symboly C[a, b], respektive C n [a, b], n ∈ N, −∞ < a < b < ∞, oznaˇcuj´ı spojit´e, respektive spojit´e a n-kr´at spojitˇe diferencovateln´e funkce na [a, b]. Je-li f ∈ C n [a, b], pak symbolem f n (a), respektive f n (b) m´ın´ım odpov´ıdaj´ıc´ı jednostrann´e limity. Symbolem Pn [a, b], −∞ < a < b < ∞ znaˇc´ım mnoˇzinu funkc´ı, kter´e jsou na [a, b] polynomy nejv´yˇse n-t´eho stupnˇe. Je-li f funkce definovan´a na (a, b), −∞ ≤ a < b ≤ ∞, a existuj´ı-li koneˇcn´e limity f v a zprava, respektive v b zleva, pak znaˇc´ıme pˇr´ır˚ ustek f na (a, b) b jako: [f ]a = lim f (x) − lim f (x). x→b−
x→a+
Symbolem k >> 0 m´ın´ım pro k re´alnou, ˇze k je “velk´a” s t´ım, ˇze v´yznam tohoto symbolu nen´ı pˇresnˇeji urˇcen. ˇ {Tn }∞ c´ı posloupnost Cebyˇ sevov´ych ortogon´aln´ıch polynom˚ u (1.druhu) n=1 znaˇ definovan´ych vztahem Tn (x) = cos(n arccos x), x ∈ [−1, 1], n ∈ N. Tak, jak je obvykl´e znaˇc´ı Γ Gamma funkci a Jn , n ∈ N, Besselovu funkci 1. druhu.
5
Kapitola 2 ´ Uvod Rychle osciluj´ıc´ı funkc´ı naz´yv´ame takovou funkci, kter´a na sv´em definiˇcn´ım oboru v “mnoha” bodech nab´yv´a sv´eho lok´aln´ıho maxima, pˇr´ıpadnˇe minima. Za pˇredpokladu, ˇze zkouman´a funkce f m´a na sv´em definiˇcn´ım oboru spojitou prvn´ı derivaci lze tedy snadno zformulovat nutnou podm´ınku k tomu, aby se jednalo o rychle osciluj´ıc´ı funkci ve tvaru: card{x ∈ (a, b); f ′ (x) = 0} >> 0, kde [a, b] leˇz´ı v definiˇcn´ım oboru f , kter´a nen´ı na ˇz´adn´em nedegenerovan´em intervalu konstantn´ı. Dalˇs´ı typickou vlastnost´ı rychle osciluj´ıc´ıch funkc´ı je to, ˇze suprema absolutn´ı hodnoty jejich derivac´ı rostou rychleji, neˇz line´arnˇe (napˇr´ıklad pro f (x) = sin nx plat´ı: sup |f (k) (x)|, x ∈ [0, 2π] = nk , n, k ∈ N). Vzhledem k tomu, ˇze u vˇetˇsiny kvadraturn´ıch formul´ı lze odhad jejich chyby vyj´adˇrit jako n´asobek hodnoty derivace interpolovan´e funkce, pˇr´ıpadnˇe jako integr´al z pˇr´ısluˇsn´e derivace, vid´ıme, ˇze u rychle osciluj´ıc´ıch funkc´ı by klasick´e postupy mohly v´est k velmi nepˇresn´ym v´ysledk˚ um. ˇ Casto lze rychle osciluj´ıc´ı funkci zapsat jako souˇcin neosciluj´ıc´ı a rychle osciluj´ıc´ı funkce, tedy f = gK, kde g je neosciluj´ıc´ı funkce, a K je rychle osciluj´ıc´ı “j´adro” funkce f (ˇcasto je K trigonometrick´a funkce, pˇr´ıpadnˇe exponenci´ala s ryze imagin´arn´ım parametrem). Ve vˇetˇsinˇe aplikac´ı jsou funkce g i K hladk´e a diferencovateln´e alespoˇ n do urˇcit´eho ˇr´adu. Jedn´ım ze z´akladn´ıch pˇr´ıklad˚ u, kde se setk´av´ame s integrac´ı rychle osciluj´ıc´ıch funkc´ı jsou Fourierovy ˇrady. Pro pˇripomenut´ı je zde uvedu v jedn´e z nejjednoduˇsˇs´ıch moˇzn´ych forem. Necht’ g ∈ C(R) je periodick´a funkce s periodou 2T , coˇz lze ekvivalentnˇe zapsat jako: g(x) = g(x+2T ), x ∈ R. Automaticky je tedy splnˇena podm´ınka
6
RT
−T
|g(x)|2dx < ∞ a plat´ı n´asleduj´ıc´ı: ∞ a0 X πnx πnx + an cos + bn sin b, 2 T T n=1 Z 1 T πnx = g(x) cos dx, T −T T Z πnx 1 T = g(x) sin dx, T −T T
g(x) = an bn
pˇriˇcemˇz konvergence trigonometrick´ych polynom˚ u je dokonce stejnomˇern´a. Bohuˇzel se ˇcasto jedn´a o funkce, u nichˇz nelze pˇresnou hodnotu jejich integr´alu spoˇc´ıst, anebo je jejich v´ypoˇcet pˇr´ıliˇs obt´ıˇzn´y, a tud´ıˇz v praxi nepouˇz´ıvan´y. C´ılem t´eto pr´ace je shrnout nejjednoduˇsˇs´ı aproximaˇcn´ı metody, uv´est jejich z´akladn´ı vlastnosti a nˇekter´e z nich aplikovat na jednoduch´ych pˇr´ıkladech. V´yklad nebude veden v nejvyˇsˇs´ı moˇzn´e obecnosti, ale v takov´e formˇe, ve kter´e jsou studovan´e metody nejˇcastˇeji uˇz´ıv´any.
7
Kapitola 3 Metody numerick´ e integrace pro obecnou osciluj´ıc´ı funkci 3.1
Pouˇ zit´ı vhodn´ e v´ ahov´ e funkce
Mˇejme d´anu funkci f ∈ C[a, b], −∞ < a < b < ∞ splˇ nuj´ıc´ı n´asleduj´ıc´ı vztah: f (x) = g(x)K(x), x ∈ [a, b], kde K je rychle osciluj´ıc´ı funkce a g je neosciluj´ıc´ı ˇc´ast integrandu. Pˇredpokl´adejme, ˇze funkce K mˇen´ı znam´enko Rb na [a, b], d´ale necht’ lze a g(x)dx snadno spoˇc´ıst, pˇr´ıpadnˇe aproximovat. V pˇr´ıpadˇe, ˇze bychom chtˇeli pouˇz´ıt Gaussovu, nebo jinou kvadraturn´ı formuli, s nez´apornou v´ahovou funkc´ı, ve kter´e by funkce K vystupovala v pˇr´ısluˇsn´e v´ahov´e funkci, je nutno pˇrej´ıt k nov´e v´ahov´e funkci L dan´e vztahem L = K + K0 , kde K0 = sup K(x), x ∈ [a, b], kter´a nemˇen´ı znam´enko. R R R Nyn´ı m˚ uˇzeme vyuˇz´ıt identity ab f (x) = ab [K(x) + K0 ]g(x)dx−K0 ab g(x)dx. Prvn´ı integr´al spoˇcteme pomoc´ı standardn´ı kvadraturn´ı formule s v´ahovou funkc´ı L, o druh´em integr´alu jsme v´yˇse pˇredpokl´adali, ˇze jeho hodnotu dok´aˇzeme alespoˇ n pˇribliˇznˇe naj´ıt. Naneˇstˇest´ı je cena t´eto metody pomˇernˇe vysok´a - nam´ısto n v´ypoˇct˚ u potˇrebujeme 2n, protoˇze uzly jsou obecnˇe jin´e. Tato metoda je t´eˇz zpracovan´a ve ([2], str. 399).
8
3.2
Vyuˇ zit´ı integrace per partes
Vˇ eta 3.1. Mˇejme d´any n ∈ N, −∞ < a < b < ∞, g ∈ C n [a, b], K ∈ C[a, b] a posloupnost {Ki }ni=0 ⊆ C[a, b] splˇ nuj´ıc´ı K0 = K, Ki′ (x) = Ki−1 (x), i = 1, ..., n, x ∈ [a, b]. Pak pro funkci f ∈ C[a, b] splˇ nuj´ıc´ı vztah f (x) = = g(x)K(x), x ∈ [a, b] plat´ı: Z
b a
f (x)dx =
n−1 X
(−1)j [g (j) (x)Kj+1 (x)]ba + (−1)n
j=0
Z
b
g (n) (x)Kn (x)dx. (3.1)
a
D˚ ukaz. Vyuˇzijeme matematick´e indukce: Pro n = 1 ovˇeˇr´ıme znˇen´ı vˇety pˇr´ımo integrac´ı per partes. Mˇejme n ∈ N, a, b, g, K, {Ki}ni=0 odpov´ıdaj´ıc´ı znˇen´ı vˇety a pˇredpokl´adejme,ˇze vˇeta plat´ı pro m = n − 1. S vyuˇzit´ım identity Z
b
a
h
g (n−1) (x)Kn−1 (x)dx = g (n−1) (x)Kn (x)
ib
a
−
Z
b
a
g (n) (x)Kn (x)dx
a d´ıky pˇredpokladu Z
b
a
f (x)dx =
n−2 X
j
(−1) [g
(j)
(x)Kj+1 (x)]ba
n−1
+ (−1)
j=0
Z
b
a
g (n−1) (x)Kn−1 (x) dx
dost´av´ame prost´ym dosazen´ım vzorec(3.1). Vzorec (3.1) m˚ uˇze platit i v pˇr´ıpadˇe, kdy interval (a, b) nen´ı omezen´y (a to i v pˇr´ıpadˇe a = −∞, b = ∞). V tomto pˇr´ıpadˇe m´ısto funkc´ı z C[a, b] uvaˇzujeme funkce spojit´e na (a, b). K tomu, aby platil vzorec (3.1) pro g maj´ıc´ı spojitou n-tou derivaci na intervalu (a, b) (tedy nejen pro g ∈ C n [a, b]) a posloupnost {Ki }ni=0 funkc´ı spojit´ych a definovan´ych na (a, b) splˇ nuj´ıc´ı pˇredpoklady vˇety na (a, b), staˇc´ı pro K = K0 , je-li integr´al na prav´e stranˇe (3.1) koneˇcn´y (tedy absolutnˇe konvergentn´ı), a existuj´ı-li koneˇcn´e limity lim g (j) (x), lim Kj+1 (x), lim g (j)(x) a lim Kj+1 (x) pro j = 0, ..., n − 1. x→a+
x→a+
x→b−
x→b−
D˚ ukaz se snadno provede pomoc´ı limitn´ıho pˇrechodu. Nyn´ı si uk´aˇzeme aplikaci vzorce (3.1) na jednoduch´ych pˇr´ıkladech.
9
Volba −∞ < a < b = ∞, g ∈ C n [a, c], a < c < ∞, K = eisx , s ∈ R vede ke vzorci: Z
∞ a
eisx g(x)dx = −eisa
n−1 X
ij−1 s−j−1g (j) (a) + (−is)−n
j=0
Z
∞
a
eisx g (n) (x)dx, (3.2)
kter´y plat´ı bez dalˇs´ıch pˇredpoklad˚ u, pokud m´a prav´a strana smysl. n Volba −∞ < a < b < ∞, g ∈ C [a, b], K = eisx , s ∈ R d´av´a vzorec: Z
b
a
eisx g(x)dx = eisb
n−1 X
ij−1 s−j−1g (j) (b) − eisa
j=0
+ (−is)−n
Z
n−1 X
ij−1 s−j−1g (j) (a) +
j=0
b
a
eisx g (n) (x)dx. (3.3)
Dle ([1], str. 149) lze z´ıskat zobecnˇen´ı pro integr´aly obsahuj´ıc´ı singularitu v x = a nebo v x = b, jeˇz je zformulov´ano v n´asleduj´ıc´ı vˇetˇe. Vˇ eta 3.2. Mˇejme n ∈ N, −∞ < a < b < ∞, g ∈ C n [a, b], 0 < µ ≤ 1, 0 < λ ≤ 1. Pak pro funkci f ∈ C[a, b] definovanou vztahem f (x) = g(x)eixs (x − a)λ−1 (b − x)µ−1 , x ∈ [a, b] plat´ı: Z
b
a
f (x)dx =
n−1 X j=0
−
Γ(j + µ) πi(j−µ)/2 −j−µ isb dj e s e [(b − a)λ−1 g(b)] j! dbj
Γ(j + λ) πi(j+λ−2)/2 −j−λ isa dj e s e [(b − a)µ−1 g(a)] + En , j! daj
En = O(s−n ). R
Vid´ıme, ˇze ve speci´aln´ım pˇr´ıpadˇe λ = µ = 0 plat´ı: En = ab eisx g (n) (x)dx = o(s−n ). Pouˇzit´e symboly o, O nejsou ch´ap´any ve sv´em obvykl´em v´yznamu, ale sp´ıˇse jako voln´e vyj´adˇren´ı faktu, ˇze pro g neosciluj´ıc´ı En sn pro n, pro nˇeˇz m´a En smysl, v´yraznˇe neroste v obecn´em pˇr´ıpadˇe. Analogicky v´yraz En = o(s−n ) vyjadˇruje, ˇze pro g ∈ C ∞ [a, b], g neosciluj´ıc´ı, tedy speci´alnˇe takovou, ˇze jej´ı 10
derivace lze stejnomˇernˇe omezit, existuje n0 ∈ N, ˇze En1sn >> 0 pro n > n0 . Nyn´ı aplikujeme tuto metodu pro aproximaci n´asleduj´ıc´ıho integr´alu. In =
Z
0
π
2
ex cos nx dx.
(3.4)
Opakovanou aplikac´ı integrace per partes dostaneme In =
2 (−1)n 2 2πeπ + 3 2 n n
Z
π
0
2
(xex )′′ sin nx dx.
(3.5)
Nyn´ı se nab´ız´ı za aproximaci In povaˇzovat prvn´ı sˇc´ıtanec vzorce (3.5), kter´y oznaˇc´ıme jako Jn . R Uv´aˇz´ıme-li, ˇze lim 0π h(x) sin(nx)dx = 0 pro kaˇzdou h ∈ C[0, π], t´ım sp´ıˇse 2
n→∞
pro h = (xex )′′ , zjist´ıme, ˇze pro chybu aproximace (kterou znaˇc´ıme En ) plat´ı: En = In − Jn = o(n3 ). (3.6) Pro velk´a n tedy m´ame velice jednoduch´y vzorec, kter´y dobˇre aproximuje In .
3.3
Numerick´ a integrace mezi nulov´ ymi body R
V tomto odstavci se budeme zab´yvat pˇribliˇzn´ym v´ypoˇctem 01 f (x)dx, kde funkci f , respektive jej´ı restrikci na interval [0,1] uvaˇzujeme z mnoˇziny M = {f ∈ C[0, 1] : card({x ∈ [0, 1] : f (x) = 0}) < ∞}. Tedy f na intervalu [0,1] nab´yv´a hodnoty nula pouze v koneˇcnˇe mnoha bodech. Mˇejme nyn´ı d´anu pevnou f ∈ M. Pak zˇrejmˇe existuje n ∈ N a koneˇcn´a posloupnost 0 ≤ t1 < t2 < ... < tn ≤ 1, kde n = card({x ∈ [0, 1] : f (x) = 0}, tedy ti jsou pr´avˇe vˇsechny nulov´e body f v [0,1]. Vyuˇzijeme rovnosti Z1 0
f (x)dx =
t n Zi+1 X
f (x)dx,
(3.7)
i=0 ti
kde t0 = 0 a tn+1 = 1. Rt Pro v´ypoˇcet tii+1 f (x) dx m˚ uˇze b´yt vhodn´e pouˇz´ıt nˇejakou kvadraturn´ı formuli (samozˇrejmˇe v Rpˇr´ıpadˇe t1 = 0, respektive tn = 1, nepoˇc´ıt´ame R t1 t f (x)dx, respektive tnn+1 f (x) dx. t0
11
Mˇejme tedy kvadraturn´ı formuli ˇr´adu m + 2 ve tvaru: Z1 0
m
X . g(x)dx = Qm+2 g = A0 g(0) + Ai g(xi ) + Am+1 g(1)
(3.8)
i=1
D´ıky substituci x = R1
y−tj , tj+1 −tj
jeˇz n´am d´av´a identitu
R tj +1 tj
f (y)dy =
= (tj+1 − tj ) 0 f (tj + x(tj+1 − tj ))dx dost´av´am pro j = 1, 2, ..., n − 1 aproximaˇcn´ı vzorec ve tvaru: tZ j +1 tj
n
X . f (y)dy = (tj+1 −tj )(A0 f (tj )+ Ai f ((tj +xi )(tj+1 −tj ))+Am+1 f (tj+1)). i=1
(3.9)
Snadno se m˚ uˇzeme pˇresvˇedˇcit, ˇze vzorec (3.9) plat´ı i pro j = 0 nebo j = n. Dosazen´ım (3.9) do (3.7) a s vyuˇzit´ım vlastnosti f (tj ) = 0, j = 1, 2, ..., n dost´av´ame hledan´y vzorec. Z1 0
f (x)dx = A0 t1 +
m X n X
Ai (tj−1 − tj )f (tj + xi (tj+1 − tj )) + An+1 (1 − tn )f (1)
i=1 j=1
(3.10)
Ve vzorci (3.10) potˇrebujeme celkovˇe mn + 2 funkˇcn´ıch hodnot. V obecn´em pˇr´ıpadˇe n´am nic nezaruˇc´ı, ˇze f (x) nen´ı rychle osciluj´ıc´ı funkce v [ti , ti+1 ]. Je-li napˇr´ıklad f (x) 6= 0 ∀x ∈ [0, 1], pak n´am popsan´a metoda d´av´a Gaussovu kvadraturn´ı metodu ˇr´adu m + 2. R Moˇzn´ym v´ychodiskem je nalezen´ı fuknce g(x) takov´e, ˇze zn´ame 01 g(x)dx a z´aroveˇ n funkce f (x) − g(x) nab´yv´a nuly bˇehem kaˇzd´e sv´e oscilace (jednou oscilac´ı m´ın´ım, volnˇe ˇreˇceno, dostateˇcnˇe velk´y interval, v nˇemˇz funkce f (x)−g(x) nab´yv´a sv´eho lok´aln´ıho maxima i minima v pr´avˇe jednom bodˇe). Pot´e pouˇziji (3.10) na funkci f (x) − g(x) a vyuˇziji identity R1 R1 R1 0 f (x)dx = 0 (f (x) − g(x))dx + 0 g(x)dx. Moˇzn´ym zdrojem informac´ı o t´eto metodˇe je ([1], str. 398-399).
12
3.4
Analytick´ a integrace osciluj´ıc´ı ˇ c´ asti
Zab´yvejme se situac´ı, kdy integrand lze z´ıskat jako souˇcin neosciluj´ıc´ı a rychle R R osciluj´ıc´ı funkce (tedy 01 f (x)dx = 01 g(x)K(x)dx, kde K(x) je osciluj´ıc´ı j´adro integrandu a g(x) je neosciluj´ıc´ı ˇc´ast integrandu). M˚ uˇze b´yt uˇziteˇcn´e aproximovat neosciluj´ıc´ı ˇcinitel pomoc´ı vhodn´ych funkc´ı tak, ˇze souˇcin osciluj´ıc´ı a aproximuj´ıc´ı funkce je analyticky integrovateln´y. Necht’ neosciluj´ıc´ı ˇc´ast integrandu g je aproximov´aR na vhodnou line´arn´ı kombinac´ı fukc´ı gi (x), i = 1, ..., n takov´ych, ˇze bi := 01 gi (x)K(x)dx lze snadno . P spoˇc´ıst, pˇr´ıpadnˇe aproximovat. Necht’ tedy g(x) = ni=1 ai gi (x). Potom: Z1 0
f (x)dx =
Z1 0
n
. X g(x)K(x)dx = ai i=1
Z1
gi (x)K(x)dx =
0
n X
ai bi .
(3.11)
i=1
Funkce gi (x) se ˇcasto vol´ı jako spline kˇrivky. Aplikace t´eto metody detailnˇe pˇredstav´ım v kapitole vˇenovan´e integraci speci´aln´ıch typ˚ u funkc´ı, je t´eˇz zpracov´ana v ([2], str. 400-401).
3.5
Eulerova transformace
Pro diference u ˇrad pouˇz´ıv´am znaˇcen´ı: ∆un = un+1 − un , ∆2 un = un+2 − 2un+1 + un , ∆k+1 un = ∆(∆k un ); n, k ∈ N0 . Pro identitu: Iun = un . Pro posunut´ı: Eun = un+1 ; E k un = un+k . Vˇ eta 3.3. Necht’ posloupnost (un )∞ aporn´ a a d´ ale necht’ ˇrada n=1 je nez´ P∞ P∞ 1 n n ı. Pak tak´e ˇrada n=0 2n+1 ∆ u0 je konvergentn´ı, n=0 (−1) un je konvergentn´ a konverguje ke stejn´emu souˇctu. D˚ ukaz. Pˇredpokl´adejme, ˇze (un )∞ nuje pˇredpoklady vˇety, tedy speci´alnˇe, n=1 splˇ ˇze je omezen´a. Matematickou indukc´ı lze snadno pro libovoln´a pˇrirozen´a ˇc´ısla k a n dok´azat n´asleduj´ıc´ı: |∆k un | ≤ sup{um , m ∈ N0 }. (3.12) Vˇsimnˇeme si, ˇze plat´ı: (I − E + E 2 − · · ·)(I + E) = I = (I + E)(I − E + E 2 − · · ·) 1 1 1 1 1 (I + ∆)(I − ∆ + ∆2 − · · ·) = I = (I − ∆ + · · ·)(I + ∆), 2 2 4 2 2 13
kde jsm´e v druh´e rovnosti vyuˇzili (3.12). Dost´av´ame tedy n´asleduj´ıc´ı vztahy: (I + E)−1 = (I − E + E 2 − · · ·) 1 1 1 (I + ∆)−1 = (I − ∆ + ∆2 − · · ·). 2 2 4 (3.13) Nyn´ı nen´ı tˇeˇzk´e pomoc´ı element´arn´ıch u ´ prav pˇr´ısluˇsn´ych oper´ator˚ u a zˇrejm´e identity E = ∆ + I odvodit znˇen´ı v´yˇse uveden´e vˇety: ∞ X
∞ X
(−1)n un =
(−1)n E n u0 = (I − E + E 2 − E 3 + · · ·)u0 =
n=0
n=0
1 1 = (I + E)−1 u0 = (2I + ∆)−1 u0 = (I + ∆)−1 u0 = 2 2 ∞ X 1 1 1 1 = (u0 − ∆u0 + ∆2 u0 − · · ·) = ∆n u0 . (3.14) n+1 2 2 4 2 n=0
Vˇ eta 3.4. Pro (un )∞ ı: ∆n u0 = n=1 plat´ m ∈ N.
Pn
n−m m=0 (−1)
n m
um , pro libovoln´e
D˚ ukaz. Vyuˇzijeme matematick´e indukce: Pro n = 1 vˇeta zˇrejmˇe plat´ı. Necht’ vˇeta plat´ı pro n − 1. Pak dost´av´ame: ∆n u0 =
n−1 X
(−1)n−1−m
m=0
= (−1)n =
n X
m=0
!
!
n−1 X n u0 + (−1)n−m 0 m=1 n−m
(−1)
!
n−1 n−1 um+1 − (−1)n−1−m um = m m
!
!
n−1 n−1 + m m−1
!!
um + un
!
n = n
n um . m
V mnoha pˇr´ıpadech ˇrada na konci (3.14) konverguje podstatnˇe rychleji, neˇz p˚ uvodn´ı zadan´a ˇrada. Samozˇrejmˇe lze zkonstruovat ˇrady, u kter´ych Eulerova transformace d´av´a pomalejˇs´ı konvergenci, neˇz klasick´e sˇc´ıt´an´ı, nicm´enˇe se jedn´a sp´ıˇse o teoretick´e pˇr´ıklady, neˇz o ˇrady, s nimiˇz se setk´av´ame 14
v praxi. Pˇri poˇc´ıt´an´ı integr´al˚ u funkc´ı, jejichˇz obor hodnot nen´ı nez´aporn´y, pˇr´ıpadnˇe nekladn´y, m˚ uˇze b´yt v´yhodn´e spoˇc´ıtat jednotliv´e kladn´e a z´aporn´e “pˇr´ır˚ ustky” zvl´aˇst’, a potom seˇc´ıst pˇr´ısluˇsnou ˇradu. Ve mnoha pˇr´ıpadech se m˚ uˇze st´at, ˇze ˇclen˚ u ˇrady je nekoneˇcnˇe mnoho, pˇr´ıpadnˇe je jich koneˇcnˇe mnoho, nicm´enˇe pˇr´ıliˇs na to, aby se realizace pˇresn´eho v´ypoˇctu dala prov´est bez velk´e z´atˇeˇze na v´ypoˇcetn´ı prostˇredky, pˇr´ıpadnˇe na ˇcas. Pˇr´ıkladem integr´alu, u nˇejˇz by se dala uplatnit Eulerova transformace je: I :=
Z∞ 0
sin x dx, x R
π(n+1) sin x kde bychom postupnˇe aproximovali hodnoty integr´al˚ u In = πn dx, x n=0,1,2,... a aproximaci v´ysledn´e hodnoty integr´alu bychom dostali jako . P I= ∞ n=0 In .
15
Kapitola 4 Metody numerick´ e integrace speci´ aln´ıch typ˚ u funkc´ı 4.1
Priceova metoda
Hledejme aproximaci x = kt dost´av´ame: Z
Rπ 0
g(t) sin kt dt, kde k ∈ N, g ∈ C 4 [0, π]. Poloˇz´ıme-li
Z
k x 1X sin x dx = k k i=1
Z
x sin x dx. k 0 (i−1)π 0 (4.1) Rπ 4 Pˇrirozenˇe se nab´ız´ı hledat aproximaci 0 gi (x) sin x dx pro gi ∈ C [0, π], i = 1, ..., k urˇcenou vztahem: π
g(t) sin kt dt =
1 k
kπ
g
iπ
g
!
x + (i − 1)π . gi (x) = g k Pot´e staˇc´ı vyuˇz´ıt identity Z
π 0
k Z π 1X g(t) sin kt dt = gi (x) sin(x + (i − 1)π) dx = k i=0 0
Z π k 1X i−1 = (−1) gi(x) sin(x) dx. k i=0 0
(4.2)
Pro pˇrehlednost zavedeme znaˇcen´ı f = gi . Necht’ xi , i=0,...,4 je ekvidistantn´ı dˇelen´ı intervalu [0, π] a fi znaˇc´ı f (xi ). 16
Hledejme nyn´ı aproximaˇcn´ı formuli ve tvaru: S(f ) =
Z
0
R
π
4 . X f (x) sin x dx = yn fn , yn ∈ R.
(4.3)
n=0
R
Vzhledem k tomu, ˇze 0π f (π − x) sin x dx = 0π f (x) sin x dx je pˇrirozen´e poˇzadovat, aby platilo S(f ) = S(h), kde h(x) = f (π − x), x ∈ [0, π], coˇz lze ekvivalentnˇe zapsat jako: y4 = y0 ; y3 = y2 , S(f ) lze tedy vyj´adˇrit ve tvaru: S(f ) = y0 (f0 + f4 ) + y1 (f1 + f3 ) + y2 f2 .
(4.4)
Definujme chybu aproximace jako: E(f ) = S(f ) −
Z
0
π
f (x) sin x dx,
→ d´ale oznaˇcme − y = (y0 , y1 , y2 )T . Podm´ınky E(xa ) = 0, a = 0, 1, 2, 3 → n´am jednoznaˇcnˇe urˇc´ı − y . Prost´ym dosazen´ım dostaneme n´asleduj´ıc´ı soustavu line´arn´ıch rovnic:
2 π 2 π π3
2 π 5 2 π 8 7 3 π 16
T 2 y 1 π π 2 y2 = 2 1 2 π π −4 4 y 3 1 3 π π 3 − 6π 8
1
(4.5)
Element´arn´ımi u ´ pravami dostaneme jednoznaˇcn´e ˇreˇsen´ı (4.5): 8 32 12 32 2 16 − 2 ; y1 = − + 2 ; y2 = − 2. π π π π π π Nyn´ı definujme Taylor˚ uv polynom 3. ˇr´adu: y0 = 1 +
(4.6)
f ′′ (0) 2 f (3) (0) 3 P (x) = f (0) + f (0)x + x + x , x ∈ [0, π]. 2 6 ′
Potom pro kaˇzd´e x ∈ [0, π] existuje ξx , ˇze f (x) − P (x) = Z(x) = Snadno nahl´edneme (d´ıky E(P ) = 0 a linearitˇe E, ˇze E(f ) = E(Z(x)) = S(Z) − 17
Z
0
π
E(x) sin x dx.
f (4) (ξx ) 4 x. 24
Odhadem chyby aproximace pˇri oznaˇcen´ı M = sup{|f (4) (x); x ∈ [0, π]} je |E(f )| ≤ cM, c ∈ R. R
V pˇredchoz´ım jsme pouˇzili identitu 0π sin x dx = 2. Dle ([1], str. 150) existuje ξ ∈ [0, π], ˇze E(f ) = −0.01002f (4) (ξ). S vyuˇzit´ım vztah˚ u (4.2), (4.4), (4.6) dostaneme hledanou aproximaˇcn´ı formuli: Z π k 1X . g(t) sin kt dt = Sk (g) = (−1)i−1 S(gi). (4.7) k i=0 0
Chybu aproximace m˚ uˇzeme vyj´adˇrit jako: Z
0
4.2
π
g(t) sin kt dt − Sk (g) ≤ c sup{|g (4) (x); x ∈ [0, π]}
Filonova metoda
Idea t´eto metody je zaloˇzena na myˇslence, jeˇz je naznaˇcena v kapitole 3.2. R Poˇc´ıtejme tedy aproximaci 01 g(x) sin(nπx) dx, n ∈ N, kde g nen´ı osciluj´ıc´ı funkce. 1 , Vezmeme m ∈ N dostateˇcnˇe velk´e, a zavedeme n´asleduj´ıc´ı znaˇcen´ı h = 2m i xi = 2m , i = 0, ..., 2m a rozdˇel´ıme interval [0, 1] na m podinterval˚ u d´elky 2h. Ty budeme znaˇcit Jk = [xk−1 , xk+1 ], k = 1, 3, 5..., 2m − 1. Pro tato k aproximujeme g na intervalu Jk ˇc´ast´ı paraboly, jeˇz v bodech xk−1 , xk a xk+1 nab´yv´a stejn´ych hodnot, jako aproximovan´a funkce. . Tedy g(x) = ak + bk x + ck x2 na Jk . T´ım z´ısk´av´ame n´asleduj´ıc´ı podm´ınky na koeficienty - pro pˇrehlednost pouˇziji v n´asleduj´ıc´ım znaˇcen´ı g(xk ) = gk : ak + bk (xk − h) + ck (xk − h)2 = gk−1 ak + bk xi + ci x2i = gk ak + bk (xk + h) + ck (xk + h)2 = gk+1. (4.8) je soustava tˇr´ı line´arn´ıch rovnic o tˇrech nezn´am´ych. Mˇejme pevn´e k ∈ {1, 3, 5, ..., 2m − 1}. Oznaˇc´ıme-li
1 xk − h (xk − h)2 A = 1 xk x2k 1 xk + h (xk + h)2 18
(4.8)
a
g − → k−1 b = gk , gk+1
lze snadno ovˇeˇrit, ˇze A je regul´arn´ı matice, a soustavu (4.8) lze ekvivalentnˇe − → → → zapsat ve tvaru A− x = b , kde − x = (ak , bk , ck )T . Jej´ı jednoznaˇcnˇe urˇcen´e ˇreˇsen´ı je: 1 [xk xk−1 fk+1 − 2xk−1 xk+1 fk + xi xk+1 fk−1 ] 2h2 −1 = [(xk + xk−1 )fk+1 − 4xi fk + (xk + xk+1 )fk−1 ] 2h2 1 = [fk+1 − 2fk + fk−1]. 2h2
ak = bk ck
(4.9)
R
Oznaˇcme Lp,k = Jk xp sin(nπx) dx, kde p = 0, 1, 2; k = 1, 3, ..., 2m − 1, z´ısk´ame pro k = 1, 3, ..., 2m − 1 n´asleduj´ıc´ı vzorec: Z
(ak + bk x + ck x2 ) sin(nπx) dx = ak L0,k + bk L1,k + ck L2,k .
(4.10)
Jk
S vyuˇzit´ım identit cos(a + b) − cos(a − b) cos(a + b) + cos(a − b) sin(a + b) − sin(a − b) sin(a + b) + sin(a − b)
= = = =
−2 sin(a) sin(b), 2 cos(a) cos(b), 2 cos(a) sin(b), 2 sin(a) cos(b),
dost´av´ame pˇri znaˇcen´ı nπi , 2m nπ b = , 2m
a =
n´asleduj´ıc´ı: L0,k =
Z
xk+1
xk−1
sin(nπx)dx = −
1 2 (cos(a + b) − cos(a − b)) = sin a sin b, nπ nπ 19
Z
xi+1 1 1 x cos nπx − sin nπx = nπ nπ xk−1 xi−1 i 1 2 = sin a sin b + cos a cos b + 2 2 cos a sin b, mnπ m nπ xi+1 2 2 1 2 = − x cos(nπx) + 2 2 x sin(nπx) + 3 3 cos(nπx) = nπ nπ nπ xi−1 i 2n2 π 2 − 4 = sin a sin b − cos a cos b + n3 π 3 nπm2 4 i 1 + 2 2 cos a sin b + sin a sin b . n π 2m 2m
L1,k =
L2,k
xk+1
x sin(nπx)dx = −
V´ysledn´y aproximaˇcn´ı vzorec lze tedy zapsat ve tvaru: Z1 0
m
. X f (x) sin(nπx)dx = (ak L0,k + bk L1,k + ck L2,k ), k = 2i − 1.
(4.11)
i=1
Hodnoty Lp,k je moˇzno aproximovat pomoc´ı koneˇcn´ych sum, coˇz ale nic nemˇen´ı na podstatˇe metody. Jedn´a se historicky o jednu z nejstarˇs´ıch metod slouˇz´ıc´ıch k aproximaci integr´al˚ u rychle osciluj´ıc´ıch funkc´ı a m˚ uˇzeme se s n´ı setkat jako s Filonovou, respektive s Filonovou-Simpsonovou metodou.
4.3
Pantisova metoda
Pˇredchoz´ı kapitola spolu s vyuˇzit´ım line´arn´ı substituce d´av´a n´avod, jak pro funkci Rg dostateˇcnˇe hladkou definovanou na intervalu [0, α] hledat aproxusobem rozˇs´ıˇril Filonovu imaci ab g(x) sin(nπx) dx. Pantis pˇrirozen´ym zp˚ metodu pro integraci na intervalu [0, ∞) rozdˇelen´ım tohoto intervalu na ˇc´asti J1 = [0, α], J2 = [α, ∞], α >> 0, aplikac´ı Filonovy metody na J1 a vyuˇzit´ım asymptotick´eho rozvoje g na J2 . M´a-li g spojit´e derivace aˇz do ˇr´adu n + 1 pro n ∈ N, kter´e v nekoneˇcnu konverguj´ı k 0, potom vyuˇzit´ım integrace per partes dost´av´ame: Z
J2
g(x) sin(nπx)dx =
n X i=0
+
Z
"
∞
0
20
cos(i) (nπx) g (x) (nπ)2i+1 (i)
g (m+1) (x)
#
x=α
(n)
cos (nπx) dx. (nπ)2n+1
(4.12)
Vzorec (4.12) plat´ı pro g takovou, ˇze integr´al na prav´e stranˇe je koneˇcn´y. Tato metoda m˚ uˇze d´av´avat dobr´e v´ysledky pro n >> 0, a m´a smysl ji pouˇz´ıvat pouze v pˇr´ıpadˇe, kdy nen´ı obt´ıˇzn´e nal´ezt pˇr´ısluˇsn´e derivace funkce g. Struˇcnˇe je tato metoda pops´ana t´eˇz v ([3], str.209).
4.4
Filonova lichobˇ eˇ zn´ıkov´ a metoda
R1 ˇ s´ıme u Reˇ ´ lohu: nal´ezt aproximaci −1 f (x)dx, kde f (x) = g(x)eikx , x ∈ [−1, 1], kde k ∈ R je konstanta, k >> 0. Pokud by jsme aplikovali sloˇzen´e lichobˇeˇzn´ıkov´e pravidlo na funkci f , dostali bychom pro N ∈ N, h = 1/N aproximaˇcn´ı vzorec ve tvaru: Z
1
−1
N X . f (x)dx = h wn f (nh),
(4.13)
n=−N
kde w−N = wN = 1/2, wn = 1, n = −N + 1, ..., N − 1. Vzhledem ke speci´aln´ı struktuˇre probl´emu je vˇsak vhodn´e volbu koeficient˚ u wn modifikovat n´asleduj´ıc´ım zp˚ usobem. Rozdˇelme interval [−1, 1] na 2N subinterval˚ u d´elky h, kter´e budeme znaˇcit [xn , xn+1 ] = [nh, (n + 1)h] = Jn , n = −N, ..., N − 1. Nyn´ı sestrojme funkci h, ˇze h je spojit´a, je line´arn´ı na Jn , n = −N, ..., N −1 a nav´ıc h(xn ) = g(xn ) pro vˇsechna n. Tyto podm´ınky funkci h jednoznaˇcnˇe urˇcuj´ı na intervalu [−1, 1]. Vid´ıme, ˇze h m˚ uˇzeme zapsat ve tvaru: g(xn+1) − g(xn )) , x ∈ Jn . (4.14) h Snadno nahl´edneme, ˇze funkce h je definovan´a korektnˇe, a to i v bodech xn , ve kter´ych je spojit´a zleva i zprava a nab´yv´a hodnoty g(xn ). Element´arn´ımi u ´ pravami dost´av´ame hledan´y aproximaˇcn´ı vzorec ve tvaru: h(x) = g(xn ) + (x − xn )
Z
1
. f (x)dx =
−1
Z
1
ikx
h(x)e
−1
w−N wn wN
=h
N X
wn g(xn ),
n=−N
(1 + ikh − eikh ) = , k 2 h2 sin(kh/2)2 = , n = −N + 1, ..., −N, (kh/2)2 (1 − ikh − e−ikh ) = . k 2 h2 21
(4.15)
Je vidˇet, ˇze pro h → 0 koeficienty ve vzorci (4.15) konverguj´ı ke konstantn´ım koeficient˚ um ve (4.13). Nicm´enˇe vzhledem k tomu, ˇze k >> 0 bude tato konvergence pomal´a. Protoˇze line´arn´ı interpolace neosciluj´ıc´ı funkce je obecnˇe pˇresnˇejˇs´ı neˇz interpolace funkce, kter´a rychle osciluje, popsan´a metoda by mˇela d´avat mnohem lepˇs´ı v´ysledky, neˇz klasick´e rovnobˇeˇzn´ıkov´e pravidlo, a to pˇresto, ˇze potˇrebuje stejn´y poˇcet funkˇcn´ıch hodnot. Tato metoda je t´eˇz informativnˇe zm´ınˇena v ([1], str.155)
4.5
Filonova spline metoda R
Mˇejme d´anu pevnou g ∈ C 3 [a, b] a hledejme aproximaci ab g(x) sin kx dx, kde k >> 0. Necht’ n ∈ N, a = x0 ≤ x1 < . . . < xn−1 < xn = b je dˇelen´ı intervalu [a, b]. Nyn´ı aproximujeme g pomoc´ı spline kˇrivky S splˇ nuj´ıc´ı n´asleduj´ıc´ı podm´ınky: (1) (2) (3) (4)
S|[xi,xi+1 ] ∈ P3 [xi , xi+1 ]. S ∈ C 2 [a, b]. S(xi ) = g(xi ), i = 1, . . . , n. S ′ (a+) = g ′ (a+); S ′(b−) = g ′(b′ )
Konstrukci funkce S, jeˇz je podm´ınkami (1)-(4) jednoznaˇcnˇe urˇcena, je moˇzn´e nal´ezt napˇr´ıklad v ([1], str. 62-67), odkud pˇreb´ır´am i pouˇzit´e znaˇcen´ı: S ′′ (xj ) = Mj , j = 0, 1, ..., n; hj = xj − xj−1 . Vyuˇzijme integrace per partes: Z
1 1 1 S(x) sin kxdx = [− S(x) cos kx − 2 S ′ (x) sin kx + 3 S ′′ (x) cos kx]ba + k k k b 1 Z b ′′′ S (x) cos kx dx. (4.16) − 3 k a Z konstrukce S(x) vid´ıme, ˇze S ′′′ (x) je po ˇc´astech konstantn´ı, a na intervalu (xj−1 ; xj ) nab´yv´a hodnoty a
Nj =
Mj − Mj−1 hj 22
. Mus´ı tedy nutnˇe platit: Z
b a
n 1X S (x) cos kxdx = Nj (cos kxj − cos kxj−1 ). k j=1 ′′′
(4.17)
Vztahy (4.16) a (4.17) n´am d´avaj´ı pomˇernˇe jednoduch´y vzorec, pomoc´ı R kter´eho lze pˇresnˇe spoˇc´ıst ab S(x) sin kx dx.
4.6
Vyuˇ zit´ı ortogon´ aln´ıch syst´ em˚ u polynom˚ u P
Necht’ pn (x) = nk=0 ak,n xk , n ∈ N je mnoˇzina polynom˚ u ortogon´aln´ıch vzhledem k nez´aporn´e v´ahov´e funkci w(x) na intervalu [a, b],tedy: Z
b
a
w(x)pn (x)pm (x) = an δn,m , m ∈ N, n ∈ N,
(4.18)
kde an jsou nenulov´e konstanty. Polynomy pn splˇ nuj´ı 3-ˇclennou rekurenci: pn+1 (x) = (An x + Bn )qn (x) + Cn qn−1 (x), n ∈ N, x ∈ [a, b], p1 (x) = (A0 x + B0 )p0 (x), (4.19) An h n An = an+1,n+1 , Cn = An−1 , an,n hn−1 Bn = An (rn+1 − rn ), rn = an−1,n an,n , n > 0, r0 = 0. Nyn´ı volme pn = Tn , tedy w(x) = (1 − x2 )−1/2 . Mˇejme d´anu funkci g, ˇze g(x)/w(x) lze na intervalu [a, b] vyj´adˇrit jako stejnomˇernˇe konvergentn´ı ˇrada ortogon´aln´ıch polynom˚ u pn . Tedy:
g(x)/w(x) =
∞ X
bn pn (x).
n=0
S vyuˇzit´ım platn´ych identit Z
sin mxTn (x) dx = (1 − x2 )1/2 −1 = Z 1 cos mxTn (x) dx = (1 − x2 )1/2 −1 = 1
0, n = 2k (−1)k πJn (m), n = 2k + 1 0, n = 2k + 1 (−1)k πJn (m), n = 2k. 23
(4.20)
Dost´av´ame n´asleduj´ıc´ı vzorce: Z
1
g(x) sin mx dx = π
−1
Z
1
g(x) cos mx dx = π
−1
∞ X
k=0 ∞ X
b2k+1 (−1k )J2k+1 (m) b2k (−1k )J2k (m)
k=0
Nyn´ı zb´yv´a urˇcit koeficienty bm : po pˇren´asoben´ı g Tm a integraci pˇres interval (−1; 1) dostaneme: Z
1
−1
g(x)Tn (x) dx =
∞ X
bn Tn (x)Tm (x)w(x) dx = bn
n=0
Z
1 −1
Tn2 (x) dx.
Z´amˇena ˇrady a integr´alu v pˇredchoz´ım plyne ze stejnomˇern´e konvergence P∞ redchoz´ıho element´arnˇe plyne jednoznaˇcnost stejnomˇernˇe n=0 bn Tn (x). Z pˇ R1 konvergentn´ıho rozvoje, kde pˇri znaˇcen´ı cn = −1 Tn2 (x) dx plat´ı: bn =
R1
−1
g(x)Tn (x) dx , n = 0, 1, ... cn
(4.21)
V t´eto f´azi staˇc´ı dopoˇc´ıtat koeficienty cn . Element´arn´ımi u ´ pravami dostaneme: c0 = 2; cn = π/2, n ∈ N. Popsan´a metoda urˇcov´an´ı koeficient˚ u bn m´a sice dobrou teoretickou hodnotu, nicm´enˇe je v praxi t´emˇeˇr nepouˇziteln´a vzhledem k n´arok˚ um na v´ypoˇcet integr´al˚ u vyskytuj´ıc´ıch se v (4.21). N´avod na prakticky realizovatelnou modifikaci t´eto metody, kde se poˇc´ıtaj´ı aproximace v´yˇse zm´ınˇen´ych koeficient˚ u mnohem efektivnˇeji, lze nal´ezt napˇr´ıklad v ([1], str. 161-163), stejnˇe jako podobnou metodu zaloˇzenou na expanzi g pomoc´ı Legenderov´ych polynom˚ u. Dalˇs´ı ot´azkou je existence koeficient˚ u bn uveden´ych vlastnost´ı. Lze uk´azat, ˇze {bn }∞ ych vlastnost´ı existuje pro g a w spojit´e a maj´ıc´ı koneˇcnou n=0 uveden´ variaci na [−1; 1] . Omez´ıme se na jednoduˇsˇs´ı pˇr´ıpad, kdy g i w jsou takov´e, ˇze g/w je re´alnˇe analytick´a funkce na [a,b], tedy existuje posloupnost {dn }∞ ze g/w je na n=0 ⊆ R, ˇ P Pn n n [−1, 1] stejnomˇern´a limita ˇrady ∞ d x . Definujeme-li D = n n=0 n n=0 dn x , Pn pak nutnˇe existuj´ı koeficienty Dn,m , m ≤ n, ˇze Dn = m=0 Dn,m pm . Vid´ıme, ˇze g/w lze libovolnˇe dobˇre aproximovat pomoc´ı ˇc´asteˇcn´ych souˇct˚ u pouˇzit´ych ortogon´aln´ıch polynom˚ u, coˇz jsme chtˇeli dok´azat.
24
4.7
Pˇ reveden´ı probl´ emu na ˇ reˇsen´ı obyˇ cejn´ e diferenci´ aln´ı rovnice
Mˇejme d´any funkce g, q ∈ C[a, b], −∞ < a < b < ∞, kde g nen´ı osciluj´ıc´ı funkce, a funkce K(x) = eiq(x) , x ∈ [a, b] je rychle osciluj´ıc´ı. Oznaˇcme f (x) = g(x)K(x), x ∈ [a, b] a hledejme hodnotu integr´alu Z
b a
f (x)dx.
(4.22)
Pˇredpokl´adejme, ˇze existuje funkce p ∈ C 1 [a, b], ˇze g(x) = iq ′ (x)p(x) + p′ (x), x ∈ [a, b],
(4.23)
potom d´ıky identitˇe [p(x)eiqx ]′ = (iq ′ (x)p(x) + p′ (x))eiq(x) = g(x)K(x) = f (x) plat´ı
Z
b
a
f (x)dx = p(b)eiq(b) − p(a)eiq(a) .
(4.24)
D´ıvejme se na (4.23) jako na diferenci´aln´ı rovnici pro p(x). Vˇsechna ˇreˇsen´ı pˇr´ısluˇsn´e homogenn´ı rovnice (p′ (x) = −iq ′ (x)p(x)) se daj´ı zapsat ve tvaru pbc (x) = ce−iq(x) , kde c ∈ R. R b hJako partikul´aˇrn´ı ˇreˇsen´ı (4.23) m˚ uˇzeme zvolit p(x) = e−iq(x) [ ax g(t)eiq(t) dt]. Z teorie obyˇcejn´ych diferenci´aln´ıch rovnic v´ıme, ˇze kaˇzd´e ˇreˇsen´ı (4.23) lze zapsat ve tvaru b pec (x) = pbc (x) + p(x) = e−iq(x)
Z
x
a
g(t)eiq(t) dt + c , c ∈ R.
(4.25)
Zd´a se, ˇze n´am pˇredchoz´ı identity nijak nepomohou, protoˇze v´yraz v (4.25) obsahuje integr´al typu (4.22). Nicm´enˇe plat´ı, ˇze kdyˇz f a q ′ nejsou osciluj´ıc´ı funkce, potom existuje neosciluj´ıc´ı ˇreˇsen´ı (4.23). Snaˇzme se nyn´ı naj´ıt vhodnou aproximaci nˇejak´eho ˇreˇsen´ı (4.23). Zvolme n ∈ N a (x0 , x1 , ..., xn ) dˇelen´ı intervalu (a, b), kde x0 = a; xn = b. Pro jednoduchost budeme pˇredpokl´adat, ˇze {xj }nj=0 je ekvidistantn´ı (xj = a + j (b − a), j=0,...,n. D´ale mˇejme vhodnou posloupnost {uk (x)}nk=0 ⊆ C 1 [a, b] n vhodnˇe zvolen´ych b´azov´ych funkc´ı. P Hledejme koeficienty {αk }nk=0 ⊆ C, aby pro funkci pn (x) = nk=0 αk uk (x) platilo: g(xj ) = iq ′ (xj )pn (xj ) + p′n (xj ), j = 0, ..., n. (4.26) 25
Pokusme se nyn´ı (4.26) pˇrepsat do maticov´eho tvaru. Oznaˇcme aj,k = iq ′ (xj )uk (xj )+u′k (xj ), necht’ A = (ai,j )ni,j=0, x = (α0 , ..., αn )T , g = (g(x0 ), ..., g(xn )T . Nyn´ı je (4.26) ekvivalentn´ı u ´ loze Ax = g. Probl´emy nastanou, kdyˇz matice A nebude regul´arn´ı. V pˇr´ıpadˇe, kdy g 6∈ rank(A) u ´ loha (4.26) nebude m´ıt pˇresn´e ˇreˇsen´ı, v opaˇcn´em pˇr´ıpadˇe ˇreˇsen´ı nebude urˇceno jednoznaˇcnˇe. Nab´ızelo by se pouˇz´ıt ˇreˇsen´ı ve smyslu nejmenˇs´ıch ˇctverc˚ u, nicm´enˇe v dalˇs´ım bude pˇredpokl´adat, ˇze matice A je regul´arn´ı. M´ame tedy jednoznaˇcnˇe urˇcen´e koeficienty α0 , ..., αn . Nyn´ı se pˇrirozenˇe nab´ız´ı urˇcit pˇribl´ıˇzen´ı (4.22) jako Z
b
a
n
n
i=0
i=0
X X . f (x)dx = eiq(b) αk uk (b) − eiq(a) αk uk (b).
(4.27)
Samozˇrejmˇe n´am nic nezaruˇc´ı, ˇze s rostouc´ı hodnotou n se nebude naˇse aproximace st´avat osciluj´ıc´ı funkc´ı. Moˇzn´ym zdrojem informac´ı o t´eto metodˇe je ([1], str.164-165).
26
Kapitola 5 Realizace vybran´ ych metod V t´eto kapitole bude 5 vybran´ych metod aplikov´ano na element´arn´ıch u ´ loh´ach, kter´e jsou voleny tak, ˇze zn´ame jejich analytick´e ˇreˇsen´ı. Zdrojov´e k´ody v Borland Pascalu 7.0 a spustiteln´e soubory jsou na pˇriloˇzen´em CD, v souborech osc1.pas, osc21.pas, osc22.pas, osc31.pas, osc33.pas, osc41.pas, osc43.pas, osc51.pas a osc53.pas, respektive .exe, kde prvn´ı ˇc´ıslovka odpov´ıd´a pouˇzit´e metodˇe. V t´eto kapitole bude vˇzdy i pˇredstavovat poˇcet iterac´ı. Jedn´a se o aproximativn´ı ˇreˇsen´ı n´asleduj´ıc´ıch u ´ loh: P∞ (−1)n S1 = n=1 n = − log 2. P (−1)n e S2 = ∞ n=0 en = 1+e . I1 =
I2 = I3 =
Rπ
R π (−10i+1)x x e dx 0 sin(10x)e dx = Im R0 R∞ π (−2i−1)x −x dx 0 sin(2x)e dx = Im 0 e R 10π sin(2x)e−x dx = 25 (1 − e−10π ). 0
27
=
=
−10eπ +10 . 101 2 . 5
1. Eulerova transformace Program ze vstupn´ıho souboru naˇcte ˇradu (un )∞ ı dva ˇcleny odeˇcte n=0 , prvn´ a od tˇret´ıho ˇclenu d´ale aplikuje Eulerovu transformaci. V´ysledn´y vzorec tedy vypad´a takto: i−2
X . Sj = u 0 − u 1 +
1
n+1 n=0 2
∆n u2 , i > 1, j = 0, 1,
kde lze k v´ypoˇctu ∆n u2 pouˇz´ıt (3.12). Jak vid´ıme, popsan´a implementace m´a konstant´ı pamˇet’ovou a kvadratickou ˇcasovou sloˇzitost. K vyj´adˇren´ı ∆n by t´eˇz bylo moˇzno vych´azet pˇr´ımo z rekurzivn´ı definice tohoto funkcion´alu, coˇz by vedlo k line´arn´ı z´avislosti pamˇet’ov´e n´aroˇcnosti na i. Nyn´ı m˚ uˇzeme srovnat chybu klasick´eho sˇc´ıt´an´ı (E1 ) a Eulerovy transformace (E2 ) pˇr´ı sˇc´ıt´an´ı ˇrady S1 : i E1 E2 5 0.11 0.0109 10 0.0525 0.000211 Pˇr´ı sˇc´ıt´an´ı ˇrady S2 byly v´ysledky n´asleduj´ıc´ı: i E1 5 4.93 ∗ 10−3 10 3.32 ∗ 10−5
E2 5.66 ∗ 10−4 3.2 ∗ 10−6
2. Vyuˇ zit´ı integrace per partes Dle vzorce (3.3) plat´ı: . I1 =
2i X
(−1)j/2+1 10−j−1eπ + (−1)j/2 10−j−2)
j=0,2,4,...
Stejnˇe tak dle (3.2): . I2 =
2i X
(−1)j/2 2−j−1
j=0,2,4,...
. Nyn´ı m˚ uˇzeme shrnout chyby jednotliv´ych metod (znaˇc´ıme E1 pro aproximaci I1 a E2 pro I2 ). i E1 4 2.2 ∗ 10−8 16 7.3 ∗ 10−12 28
E2 1.6 ∗ 10−3 9.3 ∗ 10−11
3. Numerick´ a integrace mezi nulov´ ymi body Aproximoval jsem I1 i I3 pomoc´ı sloˇzen´eho lichobˇeˇzn´ıkov´eho pravidla. V´ysledn´e chyby jsou: i E1 E3 10 1.8 ∗ 10−2 4.1 ∗ 10−3 100 1.8 ∗ 10−4 4.1 ∗ 10−5 4. Priceova metoda Pouˇzit´ım substituce a prost´ym dosazen´ım dost´av´ame pˇri aproximaci I1 chybu 1.1 ∗ 10−6 a 1.2 ∗ 10−4 pˇri aproximaci I3 . 5. Filonova lichobˇ eˇ zn´ıkov´ a metoda Vzhledem k tomu, ˇze funkce integrovan´e pˇri v´ypoˇctu integr´al˚ u I1 a I3 nab´yvaj´ı v krajn´ıch bodech integrace hodnoty 0, jedn´a se o kvadraturn´ı metodu s konstantn´ımi koeficienty d´avaj´ıc´ı n´asleduj´ıc´ı v´ysledky: i E1 100 3.7 ∗ 10−2 1000 3.6 ∗ 10−4
E3 2.9 ∗ 10−2 3 ∗ 10−4
Samozˇrejmˇe je moˇzn´e a dokonce i u ´ˇceln´e na tuto metodu aplikovat postup popsan´y v (3.3), kter´y ale na tˇechto konkr´etn´ıch pˇr´ıkladech ˇr´ad chyby pˇr´ıliˇs nezlepˇsuje. Na pˇredchoz´ıch pˇr´ıkladech vid´ıme, ˇze vybran´e metody funguj´ı s rozumnou pˇresnost´ı, a pˇri praktick´e realizaci je moˇzn´e nˇekter´e z nich vyj´adˇrit v mnohem jednoduˇsˇs´ı formˇe, neˇz pˇri obecn´em studiu jejich vlastnost´ı.
29
Literatura [1] Davis, Rabinowitz: Methods of Numerical Integration, Second Edition, Academic Press, 1984. [2] Engels: Numerical Quadrature and Cubature, Academic Press, 1980. [3] Kythe, Sch¨aferkotter: Handbook of computational methods for integration, Chapman & Hall 2005.
30
