DelodiplomskegaseminarjaMentorica:prof.dr.MarjetkaKnezLjubljana,2021 Numeriˇcnemetodezaintegracijohitrooscilirajoˇcihfunkcij UNIVERZAVLJUBLJANIFAKULTETAZAMATEMATIKOINFIZIKOFinanˇcnamatematika–1.stopnjaAljoˇsaRebolj

UNIVERZA V LJUBLJANI

FAKULTETA ZA MATEMATIKO IN FIZIKO Finanˇ cna matematika – 1. stopnja

Aljoˇsa Rebolj

Numeriˇ cne metode za integracijo hitro oscilirajoˇ cih funkcij

Delo diplomskega seminarja

Mentorica: prof. dr. Marjetka Knez

Ljubljana, 2021

Kazalo

1. Uvod 4

2. Asimptotska metoda 6

2.1. Odsotnost stacionarnih toˇck 6

2.2. Ob prisotnosti stacionarnih toˇck 11

3. Filonove metode 16

3.1. Ideja Filonove metode 17

3.2. Filon-Lagrangeve metode 17

3.3. Remesov postopek 19

3.4. Metoda najmanjˇsih kvadratov 22

4. Zakljuˇcek 25

Slovar strokovnih izrazov 25

Literatura 25

Numeriˇcne metode za integracijo hitro oscilirajoˇcih funkcij

Povzetek

V delu diplomskega seminarja si pogledamo numeriˇcne metode za integracijo hitro oscilirajoˇcih funkcij specifiˇcnega tipa. Podrobneje obravnavamo dva razreda metod, asimptotske metode in Filonove metode. V obeh primerih integrand loˇcimo na del, ki povzroˇca hitro oscilacijo, in del, ki ne oscilira oziroma oscilira poˇcasi. Asimptotske metode so primerne predvsem pri zelo hitrih oscilacijah in se loˇcijo glede na lastnosti hitro oscilirajoˇcega dela. Ideja Filonovih metod je, da pohleven del integranda aproksimiramo s polinomi, iz preostanka pa izpeljemo tako imenovane momente, ki se jih pod doloˇcenimi pogoji da eksaktno izraˇcunati. Vse obravnavane metode implementiramo in preizkusimo na veˇc numeriˇcnih primerih, pri ˇcemer primerjamo absolutno napako glede na razliˇcne hitrosti oscilacije.

Numerical methods for highly oscillatory integrals

Abstract

In this thesis, we look at numerical methods for the integration of highly oscillatory functions of a specific type. We discuss two classes of methods in detail, asymptotic methods and Filon Methods. In both cases, we divide the integrand into a part, which causes high oscillation, and a part that does not oscillate or oscillates slowly.

Asymptotic methods are particularly suitable for very high oscillations and differ according to properties of the part, which causes high oscillation. The idea of Filon methods is that we approximate the part that does not contribute to high oscillation with polynomials, and from the rest of the integrand we derive the so- called moments, which can be calculated analitically under certain conditions. We implement all the discussed methods and test them on several numerical cases, comparing the absolute error with respect to different oscillation speeds.

Math. Subj. Class. (2010): 65D15, 65D30

Kljuˇcne besede: hitro oscilirajoˇce funkcije, asimptotska metoda, Filonove metode Keywords: highly oscillatory functions, asymptotic method, Filon methods

1. Uvod

V delu diplomskega seminarja bodo predstavljene metode za numeriˇcno raˇcunanje integralov hitro oscilirajoˇcih funkcij, to je integralov oblike

(1)

∫︂ b a

f(x)e^iωg(x)dx,

kjer sta funkciji f in g gladki realni funkciji in ω ∈ R tak, da je |ω| ≫1. Integrali oblike (1) se pojavljajo na podroˇcju elektromagnetike, kvantne teorije, akustike, elektrodinamike in to je le vrh ledene gore uporabe. Ko se skuˇsa take integrale raˇcunati numeriˇcno, npr. s sestavljenimi Newton-Cotesovimi pravili, nastanejo ve- like ovire. Pri sestavljenem trapeznem pravilu za integracijo se interval, po katerem se integrira, deli na manjˇse podintervale. Ko se izraˇcuna vrednosti integranda na robovih podintervalov, se izraˇcuna ploˇsˇcino trapezov, ki nastanejo, vsota ploˇsˇcin trapezov pa predstavlja oceno za vrednost integrala. Pri hitro oscilirajoˇcih funkci- jah zavzame funkcija velik interval vrednosti na majhnem intervalu v definicijskem obmoˇcju. Razumljivo je, da ploˇsˇcina nastalih trapezov ne predstavlja dobre ocene za ploˇsˇcino integranda na podintervalih. ˇCe pogledamo od blizu, se nam tudi hitra oscilacija zdi poˇcasna oziroma, ˇce pogledamo od daleˇc, se nam poˇcasna oscilacija zdi hitra. V nekaterih primerih je raˇcunanje vrednosti integrala prezahtevno zaradi na- rave integranda, lahko pa funkcije ne poznamo v celoti in poznamo le nekaj vrednosti v razliˇcnih toˇckah. V takih primerih se uporabi interpolacija oziroma aproksimacija integranda, na kateri se lahko uporabi metode za integracijo. Za raˇcunanje integra- lov oblike (1) potrebujemo drugaˇcen pristop, saj bi za dosego zadovoljive natanˇcnosti vrednosti numeriˇcnega pribliˇzka takega integrala ob uporabi osnovnih metod za in- tegracijo potrebovali izjemno raˇcunsko moˇc in ˇcas, ki ju ponavadi nimamo.

Primer 1.1. Oglejmo si primer, ki bo prikazal naravo funkcij, ki jih obravnavamo.

Izraˇcunati ˇzelimo numeriˇcen pribliˇzek za integral

(2) I =

∫︂ 1

−1

cos(x)e^iωx2dx

pri ω = 30. Na sliki 1 si lahko pogledamo hitro oscilirajoˇco naravo integranda.

Poznamo veliko razliˇcnih metod za aproksimacijo dane funkcije. Najbolj znani aproksimacijski shemi sta zagotovo polinomska interpolacija in aproksimacija po metodi najmanjˇsih kvadratov. Poglejmo si, kako lahko z metodo najmanjˇsih kva- dratov aproksimiramo integrand v (2). Za funkcijski podprostor, kjer bomo iskali aproksimante, si izberimo

S = lin{1, x, x²,cos(x),sin(x),cos²(x),sin²(x),cos(2x),sin(2x),cos(4x),sin(4x)}.

Pri aproksimaciji funkcij po metodi najmanjˇsih kvadratov, ki bo podrobno opisana kasneje, uporabimo diskretni skalarni produkt funkcij f in g na intervalu [−1,1],

⟨f, g⟩=

N

∑︂

i=0

f(x_i)g(x_i),

ki ga v tem primeru raˇcunamo v N = 20 enakomerno razporejenih toˇckah {x_i}²⁰_i=1. Ko pogledamo aproksimacijo na sliki 2, vidimo, da je slaba in ni uporabna za na- daljnje raˇcunanje. Tudi ob veˇcjem in drugaˇcnem naboru funkcij, ki doloˇcajo nek podprostor S_k, ne bi dobili boljˇsih, bolj natanˇcnih rezultatov. Z naraˇsˇcanjem ω se oscilacija poveˇcuje v neskonˇcnost in zato za vsak izbor podprostoraSkobstajaωk, da

je za vsakω ≥ω_knapaka prevelika glede na podane kriterije. Potrebujemo robustno metodo za raˇcunanje vrednosti integrala, ki ne temelji na obiˇcajnih aproksimacijah oziroma interpolacijah integranda, niti na deljenju intervala in raˇcunanju vrednosti integrala na manjˇsih podintervalih.

Slika 1. Graf funkcije h(x) = cos(x)e^iωx2za ω = 30 na intervalu [−1,1].

Slika 2. Funkcija h in njena aproksimacija po metodi najmanjˇsih kvadratov.

♢

V nadaljevanju si bomo pogledali metodo, s katero lahko izraˇcunamo integrale oblike (1) veliko hitreje in natanˇcneje. Ima tudi fascinantno lastnost, da je metoda z veˇcanjem ω, ko funkcija hitreje oscilira, tudi vedno bolj natanˇcna, kar se na prvi pogled zdi neintuitivno. Metoda se imenuje asimptotska metoda in bo predstavljena v naslednjem poglavju. Metoda je med drugim odvisna tudi od funkcijeg, ki nastopa v eksponentu, in od tega, ali ima g kakˇsno stacionarno toˇcko na intervalu [a, b], saj se v tem primeru postopek raˇcunanja spremeni.

2. Asimptotska metoda

Pogledali si bomo metodo raˇcunanja integralov oblike (1), kjer bomo najprej pred- postavili, da g nima stacionarnih toˇck na intervalu [a, b]. V nadaljevanju pa bomo metodo nadgradili, tako da bo primerna tudi za integrale, pri katerih ima funkcijag eno ali veˇc stacionarnih toˇck na intervalu integracije. Pri izpeljavi prve metode bo postalo oˇcitno, zakaj moramo primera loˇciti.

2.1. Odsotnost stacionarnih toˇck. Predpostavimo, da funkcija g nima staciona- rih toˇck na intervalu [a, b].

Definicija 2.1. Naj bo g : R → R gladka funkcija na intervalu [a, b]. Definiramo linearen funkcional

I :C^∞[a, b]→C f ↦→

∫︂ b a

f(x)e^iωg(x)dx, ki dani gladki funkciji priredi integral, ki ga ˇzelimo izraˇcunati.

Izrek 2.2. Naj bosta funkciji f : R → R in g : R → R gladki in naj bo g brez stacionarnih toˇck na intervalu [a, b], a, b∈R, a < b. Za ω → ∞ velja

I[f] =

∫︂ b a

f(x)e^iωg(x)dx=

∞

∑︂

m=0

−1 (−iω)^m+1

[︄e^iωg(x) g^′(x)f_m(x)

]︄⃓

⃓

⃓

⃓

b

a

, kjer so

f₀(x) = f(x), f_m+1 = d dx

(︃f_m(x) g^′(x)

)︃

za m= 0,1, . . . Dokaz. Poglejmo si integral I[f ·(g^′)²] = ∫︁b

a f(x)(g^′(x))²e^iωg(x)dx. Z uporabo inte- gracije po delih dobimo

I[f ·(g^′)²] = 1 iω

(︂

f(x)g^′(x)e^iωg(x))︂

⃓

⃓

b a− 1

iω

∫︂ b a

d

dx(f(x)g^′(x))e^iωg(x)dx

= 1 iω

(︂

f(x)g^′(x)e^iωg(x))︂

⃓

⃓

b a− 1

iωI [︃ d

dx(f g^′) ]︃

. (3)

Tu postane pomembna naˇsa predpostavka, da g nima stacionarnih toˇck na intervalu [a, b]. Tako lahko pogledamo, kaj se zgodi, ko v (3) funkcijof nadomestimo z _(g^f′)².

Dobimo I

[︃ f

(g^′)² ·(g^′)² ]︃

=I[f]

= 1 iω

(︃f(x) g^′(x)e^iωg(x)

)︃ ⃓

⃓

⃓

⃓

b

a

− 1 iωI

[︄ d dx

(︃f g^′

)︃]︄

= 1 iω

(︃f(x) g^′(x)e^iωg(x)

)︃ ⃓

⃓

⃓

⃓

b

a

− 1

iωI[f₁]. Drugemu ˇclenu se moramo ˇse posvetiti. Najprej definirajmo funkcijo

G(x) = 1

iωe^iωg(x). Ker je

d

dxG(x) = g^′(x)e^iωg(x), se ˇclen I[f₁] z integracijo po delih razˇsiri v

I[f₁] =

∫︂ b a

1

g^′(x)f₁(x)g^′(x)e^iωg(x)dx

=G(x) (︃ 1

g^′(x)f₁(x) )︃ ⃓

⃓

⃓

⃓

b

a

−

∫︂ b a

d dx

(︃ 1

g^′(x)f₁(x) )︃

G(x)dx.

Dalje velja I[f₁] =1

iω

(︄e^iωg(x) g^′(x)f₁(x)

)︄⃓

⃓

⃓

⃓

b

a

− 1 iω

∫︂ b a

d dx

(︃ 1

g^′(x)f₁(x) )︃

e^iωg(x)dx

=1 iω

(︄e^iωg(x) g^′(x)f₁(x)

)︄⃓

⃓

⃓

⃓

b

a

− 1

iωI[f₂]. Sledi, da je

I[f] = 1 iω

(︃f(x) g^′(x)e^iωg(x)

)︃ ⃓

⃓

⃓

⃓

b

a

− 1 (iω)²

(︄e^iωg(x) g^′(x) f₁(x)

)︄⃓

⃓

⃓

⃓

b

a

+ 1

(iω)²I[f₂]. S ponovno integracijo po delih izpeljemo

I [f] = 1 iω

(︃f(x) g^′(x)e^iωg(x)

)︃ ⃓

⃓

⃓

⃓

b

a

− 1 (iω)²

(︄e^iωg(x) g^′(x) f1(x)

)︄⃓

⃓

⃓

⃓

b

a

+ 1

(iω)³

(︄e^iωg(x) g^′(x)f₂(x)

)︄⃓

⃓

⃓

⃓

b

a

− 1

(iω)³I[f₃].

Po enakem postopku lahko, s pomoˇcjo rekurzivnega zaporedja funkcijf_m, zapiˇsemo vrednost integrala kot vsoto

I[f] =

∞

∑︂

m=0

−1 (−iω)^m+1

[︄e^iωg(x) g^′(x)f_m(x)

]︄⃓

⃓

⃓

⃓

b

a

. (4)

□

Posledica 2.3. Ce v izpeljani formuliˇ (4) uporabimo prvih p ˇclenov vsote, dobimo aproksimacijo za vrednost integrala, ki jo oznaˇcimo s Q^A_p[f] in jo imenujemo asimp- totska metoda s p ˇcleni:

Q^A_p[f] =

p−1

∑︂

m=0

−1 (−iω)^m+1

[︄e^iωg(x) g^′(x) f_m(x)

]︄⃓

⃓

⃓

⃓

b

a

. (5)

Pri aproksimaciji vrednosti integrala z vsoto prvih p ˇclenov vidimo, da je glede na definicijo Q^A_p[f] napaka vsote enaka

Q^A_p[f]−I[f]∼O(ω^−p−1).

Priˇcakovali bi, da se bo napaka manjˇsala z naraˇsˇcanjem ˇstevila ˇclenov v vsoti (ˇstevila p), a vidimo, da se napaka manjˇsa tudi z naraˇsˇcanjem ω. Tako vidimo, da nam bo asimptotska metoda dala vedno natanˇcnejˇse rezultate takrat, ko bi ostale konvenci- onalne metode bile vedno slabˇse, z naraˇsˇcanjem ω.

Opomba 2.4. Pri dokazu izreka smo uporabili predpostavko, da funkcija g nima stacionarne toˇcke na [a, b] in smo zato lahko pogledali zamenjavo f z _g^f′2. Po pred- postavki o gladkosti vemo, da je na intervalu [a, b] funkcija g^′e^iωg zvezna. Tako po osnovnem izreku o infinitezimalnem raˇcunu sledi, da je

G(x) =

∫︂ x a

g^′(u)e^iωg(u)du odvedljiva funkcija, in velja

d

dxG(x) = g^′(x)e^iωg(x). Primer 2.5. Raˇcunamo vrednost integrala

I =

∫︂ 1

−1

cos(x)e^iωxdx, (6)

pri katerem je f(x) = cos(x) in g(x) =x. S pomoˇcjo vira [5], ki ga bomo uporabili tudi v ostalih primerih, lahko izraˇcunamo toˇcno vrednost integrala, ki je

I = (ω−1) sin (ω+ 1) + sin (ω−1)ω+ sin (ω−1)

ω²−1 ∈R.

Poglejmo si vrednost integrala (6) z uporabo asimptotske metode. Ker je g^′(x) = 1, jo lahko uporabimo, saj g nima stacionarnih toˇck na intervalu integracije. Za j = 0,1, . . . lahko izraˇcunamo kar sploˇsen ˇclen zaporedja funkcij f_m,

f4j(x) = cos(x), f4j+1(x) =−sin(x) f4j+2(x) = −cos(x) f4j+3(x) = sin(x).

Uporabimo asimptotsko metodo s tremi ˇcleni:

Q^A₃ =

2

∑︂

m=0

−1 (−iω)^m+1

[︂

e^iωxf_m(x)]︂⃓

⃓

1

−1

=1 iω

[︂

e^iωxf₀(x) ]︂⃓

⃓

1

−1+ 1 ω²

[︂

e^iωxf₁(x) ]︂⃓

⃓

1

−1− 1 iω³

[︂

e^iωxf₂(x) ]︂⃓

⃓

1

−1

= 1 iω³

[︂

ω²f₀(x)e^iωx+iωf₁(x)e^iωx−f₂(x)e^iωx]︂

⃓

⃓

1

−1

= 1 iω³

[︂

ω²cos(1)e^iω−ω²cos(1)e^−iω−iωsin(1)e^iω−iωsin(1)e^−iω + cos(1)e^iω−cos(1)e^−iω]︂

.

Ce upostevamo ˇse, da jeˇ e^iω = cos(ω) +isin(ω), se formula poenostavi v Q^A₃ = 2

ω³

(︁−ωcos(ω) sin(1) + (1 +ω²) cos(1) sin(ω))︁

.

Vidimo, da je tudi aproksimacija integrala realno ˇstevilo. Na sliki 3 vidimo gibanje napak |Q^A₃ −I| in |Q^A₃ −I| ·ω⁴ glede na vrednost uteˇzi ω. Na sliki 4 je prikazana absolutna napaka pri raˇcunanju integrala, kjer imamo razliˇcno ˇstevilo ˇclenov v vsoti.

Slika 3. Napaka pri raˇcunanju vrednosti integrala (6) z asimptotsko metodo v odvisnosti od ω ∈N.

♢

Slika 4. Napaka pri raˇcunanju integrala (6) z asimptotsko metodo glede na ˇstevilo ˇclenov v vsoti.

Primer 2.6. Poglejmo si vrednost integrala (1) za f(x) = 1,g(x) =x², a =−1 in b = 1,

I =

∫︂ 1

−1

e^iωx2, (7)

katerega toˇcno vrednost lahko izraˇcunamo in je I =

√π(i+ 1) erf(︂√

−iω)︂

√2√︁

|ω| , erf(z) = 2

√π

∫︂ z 0

e^−t2dt.

Z uporabo asimptotske metode (5) s tremi ˇcleni dobimo aproksimacijo integrala Q^A₃ = 1

4ω²ie^iω(︁

3−4ω²+ 2iω)︁

. (8)

Na sliki 5 vidimo gibanje absolutnih napak |Q^A₃ −I|, |Q^A₃ −I| ·ω⁴ in|Q^A₃ −I| ·ω^0.5 glede na vrednosti uteˇziω. Absolutna napaka se, v primerjavi s prejˇsnjim primerom, ne zmanjˇsuje dovolj hitro in zato pomnoˇzena absolutna napaka|Q^A₃ −I| ·ω⁴ naraˇsˇca ˇ

cez vse meje. V tem primeru – pri izboru prvih treh ˇclenov v vsoti – vidimo, da napaka pada s faktorjem ω^−0.5. V naslednjem poglavju si bomo pogledali metodo, ki nam pomaga raˇcunati vrednosti takih integralov. ♢

Slika 5. Napaka pri raˇcunanju integrala (7) z asimptotsko metodo v odvisnosti od ω ∈N.

2.2. Ob prisotnosti stacionarnih toˇck. Recimo, da ima funkcijag konˇcno ˇstevilo stacionarnih toˇck na intervalu [a, b]. Brez ˇskode za sploˇsnost lahko predpostavimo, da ima samo eno stacionarno toˇcko, saj lahko interval [a, b] delimo na manjˇse podin- tervale in na ta naˇcin doseˇzemo, da ima g najveˇc eno stacionarno toˇcko na vsakem podintervalu. Nato raˇcunamo integral na vsakem od podintervalov ter seˇstejemo dobljene vrednosti. Brez ˇskode za sploˇsnost lahko dodatno predpostavimo, da inte- griramo na intervalu [0,1], saj lahko vsak interval [a, b] transformiramo na [0,1].

Izrek 2.7. Naj bosta funkciji f : R → R in g : R → R gladki funkciji in naj ima g eno stacionarno toˇcko ξ stopnje r na intervalu (0,1), to je g(ξ) = g^′(ξ) =· · · = g^(r−1)(ξ) =g^(r)(ξ) = 0in g^(r+1)(ξ)̸= 0. Potem velja

I[f] =

∫︂ 1 0

f(x)e^iωg(x)dx=

r−1

∑︂

j=0

1

j!µ_j(ω, ξ)

∞

∑︂

m=0

1

(−iω)^mρ^(j)_m[f](ξ) (9)

−

∞

∑︂

m=0

(︄ e^iωg(x) (−iω)^(m+1)g^′(x)

(︁

ρ˜︁_m[f](x)−ρ˜︁_m[f](ξ))︁

)︄⃓

⃓

⃓

⃓

1

0

, kjer so

ρ˜︁_k[f](x) :=ρ_k[f](x)−

r−1

∑︂

j=0

ρ_k[f]^(j)(ξ)

j! (x−ξ)^j, k= 0,1,2, . . . ρ0[f](x) :=f(x)

ρ_k+1[f](x) := d dx

(︃

ρ˜︁_k[f](x) g^′(x)

)︃

, k ≥0 µ_j(ω, ξ) :=

∫︂ 1 0

(x−ξ)^je^iωg(x)dx, j ≥0.

Dokaz. V integral dodamo razvoj funkcije f v Taylorjevo vrsto okoli ξ do stopnje r−1:

∫︂ 1 0

f(x)e^iωg(x)dx=

∫︂ 1 0

⎛

⎝f(x) +

r−1

∑︂

j=0

f^(j)(ξ)

j! (x−ξ)^j−

r−1

∑︂

j=0

f^(j)(ξ)

j! (x−ξ)^j

⎞

⎠

·e^iωg(x)dx

=

∫︂ 1 0

r−1

∑︂

j=0

f^(j)(ξ)

j! (x−ξ)^je^iωg(x)dx +

∫︂ 1 0

⎛

⎝f(x)−

r−1

∑︂

j=0

f^(j)(ξ)

j! (x−ξ)^j

⎞

⎠e^iωg(x)dx

=

r−1

∑︂

j=0

f^(j)(ξ) j!

∫︂ 1 0

(x−ξ)^je^iωg(x)dx

+

∫︂ 1 0

⎛

⎝f(x)−

r−1

∑︂

j=0

f^(j)(ξ)

j! (x−ξ)^j

⎞

⎠e^iωg(x)dx.

Vstavimo funkcijo iz izreka:

ρ˜︁₀[f](x) =f(x)−

r−1

∑︂

j=0

f^(j)(ξ)

j! (x−ξ)^j.

Pomembno je pripomniti, da velja ρ˜︁₀[f](ξ) = ρ˜︁^′₀[f](ξ) = · · · = ρ˜︁^(r−1)₀ [f](ξ) = 0.

Torej ima funkcija ρ˜︁₀[f] pri x = ξ niˇclo stopnje r −1. Drugi ˇclen pomnoˇzimo z

g^′(x)iω

g^′(x)iω in dobimo

∫︂ 1 0

f(x)e^iωg(x)dx=

r−1

∑︂

j=0

f^(j)(ξ) j!

∫︂ 1 0

(x−ξ)^je^iωg(x)dx

+ 1 iω

∫︂ 1 0

(︂

f(x)−∑︁r−1 j=0

f^(j)(ξ)

j! (x−ξ)^j )︂

g^′(x) iωg^′(x)e^iωg(x)dx

=

r−1

∑︂

j=0

f^(j)(ξ) j!

∫︂ 1 0

(x−ξ)^je^iωg(x)dx (10)

+ 1 iω

∫︂ 1 0

ρ˜︁0[f](x) g^′(x)

d dx

(︂

e^iωg(x) )︂

dx.

Ker imata takoρ˜︁0[f] kotg^′ prix=ξniˇclo stopnjer−1, je kvocient ^˜︁ρ0_g^[f](x)′(x) funkcija, ki je na [0,1] dobro definirana, to pomeni, da nima polov na [0,1]. Po definiciji v izreku je

ρ0[f](x) =f(x).

Ce sedaj na drugem ˇˇ clenu v (10) uporabimo integracijo po delih, dobimo

∫︂ 1 0

f(x)e^iωg(x)dx=

r−1

∑︂

j=0

ρ0[f]^(j)(ξ) j!

∫︂ 1 0

(x−ξ)^je^iωg(x)dx (11)

+ e^iωg(x) iω

ρ˜︁0[f](x) g^′(x)

⃓

⃓

⃓

⃓

⃓

1

0

− 1 iω

∫︂ 1 0

d dx

(︃

ρ˜︁0[f](x) g^′(x)

)︃

e^iωg(x)dx.

Sedaj si pomagamo s funkcijami ρ˜︁k[f](x) = ρk[f](x)−

r−1

∑︂

j=0

ρ_k[f]^(j)(ξ)

j! (x−ξ)^j, k = 0,1,2, . . .

Funkcije ρ˜︁_k[f] imajo pri x = ξ niˇclo stopnje r −1, zato je kvocient ^ρ˜︁k_g^[f′(x)^](x) dobro definiran za vsak x∈[0,1]. Po definiciji v izreku velja

ρ_k+1[f](x) = d dx

(︃

ρ˜︁_k[f](x) g^′(x)

)︃

. Ce to vstavimo v (11), dobimoˇ

∫︂ 1 0

f(x)e^iωg(x)dx=

r−1

∑︂

j=0

ρ₀[f]^(j)(ξ) j!

∫︂ 1 0

(x−ξ)^je^iωg(x)dx (12)

+e^iωg(x) iω

ρ˜︁₀[f](x) g^′(x)

⃓

⃓

⃓

⃓

⃓

1

0

− 1 iω

∫︂ 1 0

ρ₁[f](x)e^iωg(x)dx.

Opazimo, da je integrand v zadnjem ˇclenu (12) enake oblike kot leva stran enaˇcbe (10), zato ga lahko razvijemo v

− 1 iω

∫︂ 1 0

ρ₁[f](x)e^iωg(x)dx=− 1 iω

r−1

∑︂

j=0

ρ₁[f]^(j)(ξ) j!

∫︂ 1 0

(x−ξ)^je^iωg(x)dx

− 1 (iω)²

∫︂ 1 0

ρ˜︁₁[f](x) g^′(x)

d dx

(︂

e^iωg(x))︂

dx.

Dobljen drugi ˇclen zopet integriramo po delih in dobimo

− 1 iω

∫︂ 1 0

ρ1[f](x)e^iωg(x)dx=− 1 iω

r−1

∑︂

j=0

ρ₁[f]^(j)(ξ) j!

∫︂ 1 0

(x−ξ)^je^iωg(x)dx

− e^iωg(x) (iω)²

ρ˜︁₁[f](x) g^′(x)

⃓

⃓

⃓

⃓

⃓

1

0

+ 1

(iω)²

∫︂ 1 0

ρ˜︁₂[f]e^iωg(x)dx, torej je

∫︂ 1 0

f(x)e^iωg(x)dx=

⎛

⎝

r−1

∑︂

j=0

(︄ρ₀[f]^(j)(ξ)

j! − 1

iω

ρ₁[f]^(j)(ξ) j!

)︄⎞

⎠

∫︂ 1 0

(x−ξ)^je^iωg(x)dx

+ e^iωg(x) g^′(x)

(︃

ρ˜︁₀[f](x)

iω − ρ˜︁₁[f](x) (iω)²

)︃⃓

⃓

⃓

⃓

⃓

1

0

+ 1

(iω)²

∫︂ 1 0

ρ₂[f](x)e^iωg(x)dx.

Celoten postopek ponavljamo in z vsakim korakom pridobimo ˇclen (−1)^m 1

(iω)^m

r−1

∑︂

j=0

ρ_m[f]^(j)(ξ) j!

∫︂ 1 0

(x−ξ)^je^iωg(x)dx, m∈N. Vsaka integracija po delih pridobi ˇclen oblike

(−1)^m e^iωg(x) (iω)^m+1

ρ˜︁_m[f](x) g^′(x)

⃓

⃓

⃓

⃓

⃓

1

0

, m∈N. Po definiciji v izreku je

µ_j(ω, ξ) =

∫︂ 1 0

(x−ξ)^je^iωg(x)dx, j ≥0.

Tako sledi (9). □

Opomba 2.8. Pogojg(ξ) = 0 lahko izpustimo iz pogojev izreka, saj lahko raˇcunamo vrednost integrala

e^iωg(ξ)

∫︂ 1 0

f(x)eiω(g(x)−g(ξ))

dx.

Opomba 2.9. Vrednosti ρ_m[f]^(j)(ξ) je treba raˇcunati z limitami, kjer si pomagamo z l’Hospitalovim pravilom.

Opomba 2.10. Raˇcunanje integrala (1) na intervalu [a, b], ko imageno stacionarno toˇcko na intervalu (a, b), raˇcunamo z uporabo enakosti

∫︂ b a

f(x)e^iωg(x)dx= (b−a)

∫︂ 1 0

f(u)e˜︁ ^iω˜︁g(u)du, pri ˇcemer sta

f(x) =˜︁ f(a+x·(b−a)), ing˜︁(x) = g(a+x·(b−a)).

Posledica 2.11. Ce v formuliˇ (9) uporabimo prvih pˇclenov, dobimo aproksimacijo za vrednost integrala, ki jo oznaˇcimo s Q^A_p[f] in jo imenujemo asimptotska metoda s p ˇcleni, kjer ima g stacionarno toˇcko ξ∈(0,1) stopnje r:

Q^A_p[f] =

r−1

∑︂

j=0

1

j!µ_j(ω, ξ)

p−1

∑︂

m=0

1

(−iω)^mρ_m[f]^(j)(ξ)

−

p−1

∑︂

m=0

(︄ e^iωg(x) (−iω)^m+1g^′(x)

(︁

ρ˜︁m[f](x)−ρ˜︁m[f](ξ))︁

)︄⃓

⃓

⃓

⃓

1

0

.

Poglejmo si, kakˇsna je napaka pri asimptotski metodi ob prisotnosti stacionarnih toˇck.

Trditev 2.12. Za integral oblike (1), kjer staf in g realni gladki funkciji in ima g eno stacionarno toˇckoξ∈(a, b)stopnjer,f pa je neniˇcelna le na majhnem intervalu okrog toˇcke ξ, velja naslednja asimptotiˇcna ocena (za velike ω):

∫︂ b a

f(x)e^iωg(x)dx∼

∞

∑︂

j=0

a_j ω^j+1r+1

.

Pri tem so koeficienti aj odvisni le od konˇcno mnogo odvodov funkcij f in g.

Dokaz trditve 2.12 se nahaja v virih [3] in [2]. Z uporabo izpeljav sledi, da je napaka oblike

Absolutna napaka =

⃓

⃓

⃓

⃓

⃓

∫︂ b a

f(x)e^iωg(x)dx−Q^A_p[f, g]

⃓

⃓

⃓

⃓

⃓

∼O(ω^−p−r+11 ), ko ω → ∞.

Primer 2.13. Preizkusimo izpeljano metodo na primeru, ko je g(x) = (x−0.5)², torej ima funkcija g eno stacionarno toˇcko ξ = 0.5 stopnje r = 1 na [0,1], za f pa izberemo funkcijo cos(x). Tako dobimo

(13) I =

∫︂ 1 0

cos(x)e^{iω(x−0.5)2}dx.

Raˇcunamo vrednost integrala, za katerega lahko izraˇcunamo analitiˇcno reˇsitev.

Toˇcna vrednost integrala v odvisnosti od ω je enaka

I =−√ π

⎛

⎝erf

(︄(ω−1)√

−iω 2ω

)︄

+ erf (︄√

−iω(ω+ 1) 2ω

)︄⎞

⎠·

(i−1) sin(︁_2ω+1

4ω

)︁+ (−i−1) cos(︁_2ω+1

4ω

)︁+ (1−i) sin(︁_2ω−1

4ω

)︁+ (−i−1) cos(︁_2ω−1

4ω

)︁

2⁵²√︁

|ω|

Poglejmo si aproksimacijo integrala z uporabo asimptotske metode ob prisotnosti stacionarnih toˇck za p = 2. Ker ima g stacionarno toˇcko ξ = 0.5 stopnje r = 1, potrebujemo le µ0, ki se ga da v tem primeru analitiˇcno izraˇcunati:

µ₀(ω,0.5) =

∫︂ 1 0

(x−0.5)⁰e^{iω(x−0.5)2}dx=−

√π(i+ 1) erf (︂ iω

2√

−iω

)︂

√2√︁

|ω| .

Pri izbiri p = 2 potrebujemo ˇse ˇclena ρ_k[f] in ρ˜︁_k[f] za k = 0,1. Ker je stacionarna toˇcka stopnje r= 1, se nam predpis za ρ˜︁k[f] poenostavi v

ρ˜︁_k[f](x) =ρ_k[f](x)−ρ_k[f](0.5).

Dobimo

ρ₀[f](x) = cos(x), ρ˜︁₀[f](x) = cos(x)−cos(0.5),

ρ₁[f](x) = −(2x−1) sin(x)−2(cos(x)−cos(0.5))

(2x−1)² .

Z uporabo l’Hospitalovega pravila izraˇcunamo ρ₁[f](0.5) = −¹₄cos(0.5), od koder dobimo

ρ˜︁₁[f](x) = −sin (x) 2x−1 −2(︁

cos (x)−cos (0.5))︁

(2x−1)² +1

4cos(0.5).

Vrednost aproksimacije integrala v odvisnosti od ω je tako Q^A₂ =

0

∑︂

j=0

1

j!µ_j(ω, ξ)

1

∑︂

m=0

1

(−iω)^mρ_m[f]^(j)(ξ)

−

1

∑︂

m=0

(︄ e^iωg(x) (−iω)^m+1g^′(x)

(︁

ρ˜︁_m[f](x)−ρ˜︁_m[f](ξ))︁

)︄⃓

⃓

⃓

⃓

1

0

=µ0(ω,0.5)· (4ω−i) cos(0.5) 4ω

− e^0.25iωcos(0.5)(8 cos(0.5) + 4 sin(0.5)−9−8iωsin²(0.25))

2ω² .

Na sliki 6 imamo primerjavo absolutne napake pri raˇcunanju vrednosti integrala glede na ˇstevilo ˇclenov v vsoti, ki jih vzamemo. ♢

Slika 6. Napaka pri raˇcunanju integrala (13) z asimptotsko metodo ob prisotnosti stacionarnih toˇck v odvisnosti od ω.

3. Filonove metode

Vˇcasih moramo funkcijo f v (1) aproksimirati oziroma morda sploh ne poznamo funkcije f, poznamo le njene vrednosti v nekaj toˇckah, in jo moramo zato inter- polirati. V takih primerih bomo uporabili Filonove metode. Pogledali si bomo tri naˇcine za interpolacijo oziroma aproksimacijo:

• polinomska interpolacija,

• Remesov postopek,

• metoda najmanjˇsih kvadratov.

3.1. Ideja Filonove metode. Naj bo {τ_k}^N_k=0 mnoˇzica takihN + 1 toˇck na [a, b], da je a = τ₀ < τ₁ < τ₂ < · · · < τ_N =b. Najbolj znana izbora toˇck sta ekvidistan- tne toˇcke in ˇCebiˇseve toˇcke. Pri vseh nadaljnjih metodah za aproksimacijo bomo uporabili ˇCebiˇseve toˇcke,

{τ_k}^N_k=0 =

{︄a+b

2 +b−a 2 cos

(︃N −k

N π

)︃}︄N

k=0

, (14)

razen ˇce drugaˇce ne specificiramo. Naj bo {ϕj}ⁿ_j=0 mnoˇzica zveznih funkcij, ki sli- kajo iz [a, b] v R in zadoˇsˇcajo pogoju, da je matrika (︁

ϕ_j(x_i))︁n

i,j=0 nesingularna za poljuben izbor paroma razliˇcnih toˇckx_i ∈[a, b], i= 0,1, . . . , n. Ta pogoj se imenuje Haarov pogoj in implicira obstoj enoliˇcnega interpolacijskega polinoma nad mnoˇzico interpolacijskih toˇck{xi}ⁿ_i=0. Namesto, da bi aproksimirali celoten integrandf e^iwg, bomo aproksimirali samo f. Aproksimant funkcije f oznaˇcimo sf˜︁in naj bo oblike

f˜︁(x) =

n

∑︂

j=0

c_jϕ_j(x),

pri ˇcemer bodo koeficienti (c_j)ⁿ_j=0 doloˇceni na podlagi vrednosti funkcije f v toˇckah (τ_k)^N_k=0. Tedaj je Filonova integracijska metoda oblike

Q^F_n[f] =

∫︂ b a

f(x)e˜︁ ^iωg(x)dx

=

∫︂ b a

n

∑︂

j=0

cjϕj(x)e^iωg(x)dx

=

n

∑︂

j=0

cj

∫︂ b a

ϕj(x)e^iωg(x)dx.

(15)

Na prvi pogled napredek morda ni viden, a ˇce izberemo funkcije ϕ_j tako, da po- znamo toˇcno vrednost integralov ∫︁b

a ϕ_j(x)e^iωg(x)dx, ki jih imenujemo momenti, po- tem lahko vsoto integralov (15) enostavno izraˇcunamo. Funkcije {ϕ_j}ⁿ_j=0 izberemo glede na funkcijo g. Od funkcije g je tudi odvisno, ali sploh lahko izberemo tako mnoˇzico funkcij{ϕ_j}ⁿ_j=0, da dobimo eksaktne vrednosti momentov. ˇCe to ni mogoˇce, lahko momente vseeno izraˇcunamo s pomoˇcjo asimptotske metode, ki smo jo opisali v prejˇsnjem poglavju. Pogledali si bomo nekaj razliˇcnih metod za aproksimacijo funkcije f.

3.2. Filon-Lagrangeve metode. Pri Filon-Lagrangevi metodi uporabimo poli- nomsko interpolacijo za raˇcunanje f˜︁. Pri tem mora veljati N = n. Naˇs cilj je izraˇcunati interpolacijski polinom, ki zadoˇsˇca pogojem f(τ˜︁ ℓ) =f(τℓ), ℓ = 0, . . . , n.

Ce pri polinomski interpolaciji uporabimo standardno bazoˇ {ϕj(x)}ⁿ_j=0 ={︁

1, x, x², . . . , xⁿ}︁

,

je interpolacijski polinom f˜︁oblike f˜︁(x) =

n

∑︂

j=0

cjx^j, pri ˇcemer so koeficienti c_j pridobljeni z reˇsitvijo sistema

⎡

⎢

⎢

⎢

⎣

1 τ₀ τ₀² . . . τ₀ⁿ 1 τ₁ τ₁² . . . τ₁ⁿ ... . .. ... 1 τn τ_n² . . . τ_nⁿ

⎤

⎥

⎥

⎥

⎦

⎡

⎢

⎢

⎢

⎣ c₀ c₁ ... cn

⎤

⎥

⎥

⎥

⎦

=

⎡

⎢

⎢

⎢

⎣ f(τ₀) f(τ₁)

... f(τn)

⎤

⎥

⎥

⎥

⎦ .

Opomba 3.1. Pri velikih n imamo lahko teˇzave pri reˇsevanju sistema linearnih enaˇcb z Vandermondovo matriko, saj je le ta lahko zelo obˇcutljiva (ima veliko pogo- jenostno ˇstevilo), zato se ponavadi uporablja Lagrangeeva oblika interpolacijskega polinoma. Vseeno bomo raˇcunali z Vandermondovo matriko, saj so v tem primeru momenti enaki ∫︁b

a x^je^iωg(x)dx in se jih da za doloˇcene izbire funkcije g enostavno eksaktno izraˇcunati.

Lahko si pogledamo napako metode v primeru, ko g nima stacionarne toˇcke na intervalu [a, b]. Napaka je enaka

Absolutna napaka =

⃓

⃓

⃓

⃓

⃓

∫︂ b a

f(x)e^iωg(x)dx−

∫︂ b a

f˜︁(x)e^iωg(x)dx

⃓

⃓

⃓

⃓

⃓

=

⃓

⃓

⃓

⃓

⃓

∫︂ b a

(︂

f(x)−f(x)˜︁ )︂

e^iωg(x)dx

⃓

⃓

⃓

⃓

⃓ .

Ker g nima stacionarne toˇcke na intervalu integracije, lahko vstavimo (4) namesto I[f] in I[f˜︁] in dobimo

Absolutna napaka =

⃓

⃓

⃓

⃓

⃓

⃓ 1 iω

[︄e^iωg(b) g^′(b)

(︂

f(b)−f˜︁(b))︂

− e^iωg(b) g^′(a)

(︂

f(a)−f˜︁(a))︂

]︄

(16)

+

∞

∑︂

m=1

−1 (−iω)^m+1

[︄e^iωg(x) g^′(x)

(︂

f_m(x)−f˜︁_m(x))︂

]︄b

a

⃓

⃓

⃓

⃓

⃓

⃓ .

Ker smo pri izbiri interpolacijskih toˇck zahtevali τ₀ = a in τ_n = b, po definiciji polinomske interpolacije velja f(a) =f˜︁(a) in f(b) =f˜︁(b) in tako je prvi ˇclen v (16) enak 0. Vidimo, da je absolutna napaka oblikeO(ω⁻²). ˇCe bi interpolant spremenili tako, da bi se ujemal tudi v odvodih, bi napako v odvisnosti od ω ˇse zmanjˇsali.

Primer 3.2. Poglejmo si primer integrala oblike (17)

∫︂ b a

f(x)e^iωxdx,

kjer je f zvezna funkcija. Vidimo, da je g(x) = x, in v tem primeru si ob naboru funkcij ϕj(x) =x^j, j = 0,1, . . . , n,pomagamo z eksaktno vrednostjo integrala

Mj(x) :=

∫︂

x^je^iωxdx=− (︃i

ω )︃1+j

Γ(1 +j,−iωx), j ∈N,

kjer je

Γ(j, z) =

∫︂ ∞ z

t^j−1e^−tdt.

Te integrale lahko za j ∈Nraˇcunamo rekurzivno z uporabo integracije po delih. Za j = 0,1, . . . ,5 dobimo:

M0(x) =−ie^iωx ω ,

M₁(x) =−i(ωx+i)e^iωx

ω² ,

M₂(x) =−

(︁iω²x²−2ωx−2i)︁

e^iωx

ω³ ,

M₃(x) =−

(︁iω³x³−3ω²x²−6iωx+ 6)︁

e^iωx

ω⁴ ,

M₄(x) =−i(︁

ω⁴x⁴+ 4iω³x³−12ω²x²−24iωx+ 24)︁

e^iωx

ω⁵ ,

M₅(x) =−

(︁iω⁵x⁵−5ω⁴x⁴−20iω³x³+ 60ω²x²+ 120iωx−120)︁

e^iωx

ω⁶ .

Tako lahko izrazimo

Q^F_n[f] =

∫︂ b a

f(x)e˜︁ ^iωxdx

=

∫︂ b a

n

∑︂

j=0

c_jx^je^iωxdx

=

n

∑︂

j=0

c_j

∫︂ b a

x^je^iωxdx

=

n

∑︂

j=0

c_j(︁

M_j(b)−M_j(a))︁

. (18)

Z uporabo (18) lahko ocenimo vsak integral oblike (17). Recimo, da nas zanima vrednost integrala (6),

I =

∫︂ 1

−1

cos(x)e^iωxdx.

Primerjajmo absolutno napako pri raˇcunanju vrednosti tega integrala za razliˇcne stopnje interpolacijskih polinomov. Na sliki 7 vidimo, da se z naraˇsˇcanjem ω in stopnje interpolacijskega polinoma absolutna napaka zmanjˇsuje, a pride zaradi so- dosti funkcije cos(x) do pojava, da je interpolacijski polinom sode stopnje (2k) bolj natanˇcen kot interpolacijski polinom, ki ima eno stopnjo veˇc (2k+ 1). ♢ 3.3. Remesov postopek. Pri Remesovem postopku iˇsˇcemo polinom najboljˇse ena- komerne aproksimacije za f na [a, b], ki je doloˇcen kot reˇsitev problema

minp∈Pn

∥f−p∥_∞,[a,b],

pri ˇcemer namPⁿ oznaˇcuje podprostor polinomov stopnje manj ali enakon.

Slika 7. Absolutna napaka pri aproksimaciji integrala (6) v odvi- snosti od ω ∈ N,1 ≤ ω ≤ 100, kjer s polinomsko interpolacijo inter- poliramo funkcijo cos(x), glede na razliˇcne stopnje interpolacijskega polinoma pri izbiri ˇCebiˇsevih toˇck (14).

Izrek 3.3. Naj bo f ∈C([a, b]). ˇCe je polinomp∈Pn takˇsen, da residual r=f−p doseˇze svojo normo∥r∥_∞,[a,b], alternirajoˇce v vsaj n+ 2 toˇckah τ_ℓ ∈ [a, b], a ≤ τ₀ <

· · · < τ_n+1 ≤ b, potem je p polinom najboljˇse enakomerne aproksimacije za f na [a, b].

Dokaz. Recimo, dap ni polinom najboljˇse enakomerne aproksimacije za f na [a, b].

Potem obstaja q∈Pn, da je

⃓⃓f(τ_ℓ)−q(τ_ℓ)⃓

⃓≤ max

x∈[a,b]

⃓⃓f(x)−q(x)⃓

⃓

=∥f−q∥_∞,[a,b]<∥f −p∥_∞,[a,b]=⃓

⃓f(τ_ℓ)−p(τ_ℓ)⃓

⃓. Torej velja

⃓⃓f(τ_ℓ)−q(τ_ℓ)⃓

⃓<⃓

⃓f(τ_ℓ)−p(τ_ℓ)⃓

⃓ in zato

−sign(f(τℓ)−p(τℓ))(f(τℓ)−p(τℓ))< f(τℓ)−q(τℓ)<sign(f(τℓ)−p(τℓ))(f(τℓ)−p(τℓ)).

Podobno velja

sign(f(τ_ℓ)−p(τ_ℓ))(f(τ_ℓ+1)−p(τ_ℓ+1))< f(τ_ℓ+1)−q(τ_ℓ+1)

<−sign(f(τ_ℓ)−p(τ_ℓ))(f(τ_ℓ+1)−p(τ_ℓ+1)),

saj smo predpostavili, da residual f − p na mnoˇzici toˇck (τ_ℓ)ⁿ⁺²_ℓ=0 alternira v pred- znaku. Brez ˇskode za sploˇsnost naj bo sign(f(τ_ℓ)−p(τ_ℓ)) = 1. Sledi

f(τ_ℓ)−q(τ_ℓ)< f(τ_ℓ)−p(τ_ℓ) p(τ_ℓ)−q(τ_ℓ)<0

in

f(τ_ℓ+1)−p(τ_ℓ+1)< f(τ_ℓ+1)−q(τ_ℓ+1) p(τ_ℓ+1)−q(τ_ℓ+1)>0.

Razlika p−q ∈Pn ima na (τ_ℓ, τ_ℓ+1), ℓ = 0,1, . . . , n vsaj eno niˇclo, torej mora imeti p−qvsajn+ 1 razliˇcnih niˇcel. To je moˇzno le, ˇce jep≡q, kar vodi v protislovje. □

Iskanje polinoma najboljˇse enakomerne aproksimacije se prevede na iskanje mnoˇzice toˇck {τ₀, . . . , τ_n+1}, na kateri residual alternirajoˇce doseˇze svojo normo.

Algorithm 1 Remesov postopek

1: Vhodni podatki: funkcija f, interval [a, b], stopnjan, toleranca ϵ

2: Izberi mnoˇzico toˇck E0 ={τℓ, a≤τ0 <· · ·< τn+1 ≤b}

3: konec = 0

4: k = 0

5: while konec = 0 do

6: Poiˇsci polinom p^∗_k∈Pn, ki zadoˇsˇca pogojem

f(τℓ)−p^∗_k(τℓ) = (−1)^ℓm, ℓ = 0,1, . . . , n+ 1.

(19)

7: Poiˇsˇci ekstrem residuala r_k =f −p^∗_k, to je poiˇsˇci u∈[a, b], da bo

⃓⃓r_k(u)⃓

⃓=∥r_k∥_∞,[a,b].

8: if ⃓

⃓r_k(u)⃓

⃓− |m|< ϵ then ▷ Ce je napaka dovolj majhnaˇ

9: konec = 1 ▷ konˇcaj postopek

10: else

11: Naredimo zamenjavo toˇck v mnoˇzici Ek, tako da v nabor vstavimo toˇcko uin odstranimo tisto sosednjo toˇcko toˇckeu, tako da se ohrani alternacija residuala. Dobimo novo mnoˇzicoE_k+1 ={τ_ℓ, a≤τ₀ <· · ·< τ_n+1 ≤b}in poveˇcamo k za 1.

Opomba 3.4. Ko reˇsujemo sistem enaˇcb (19) v algoritmu, iˇsˇcemo polinom p∈Pⁿ, ki se najbolje prilega funkciji f v diskretnem pribliˇzku za neskonˇcno normo:

p∈minPn

i=0,1,...,n+1max

⃓⃓f(τ_i)−p(τ_i)⃓

⃓.

Neznanke v tem sistemu n+ 2 linearnih enaˇcb so koeficienti polinoma p^∗_k in m.

Ce zapiˇsemoˇ p^∗_k =c₀+c₁x+· · ·+c_nxⁿ, je sistem enaˇcb (19) v matriˇcni obliki enak

⎡

⎢

⎢

⎢

⎢

⎢

⎣

1 1 τ₀ . . . τ₀ⁿ

−1 1 τ₁ . . . τ₁ⁿ 1 1 τ₂ . . . τ₂ⁿ ... ... ... . .. ... (−1)ⁿ⁺¹ 1 τn+1 . . . τ_n+1ⁿ

⎤

⎥

⎥

⎥

⎥

⎥

⎦

⎡

⎢

⎢

⎢

⎢

⎢

⎣ m c₀ c₁ ... cn

⎤

⎥

⎥

⎥

⎥

⎥

⎦

=

⎡

⎢

⎢

⎢

⎢

⎢

⎣ f(τ₀) f(τ₁) f(τ₂)

... f(τn+1)

⎤

⎥

⎥

⎥

⎥

⎥

⎦ .

Remesov algoritem je iterativen postopek, ki ga ponavljamo, dokler polinom ni izraˇcunan na predpisano natanˇcnost. Da se pokazati, da postopek vedno konver- gira [1]. V naˇsih primerih bomo za doloˇcitev aproksimanta f˜︁izvedli le en korak Remesovega postopka, pri ˇcemer bomo za zaˇcetno mnoˇzico toˇck E₀ izbrali ˇCebiseve toˇcke (14) z N =n+ 1. Tako dobimo primerljivo ˇcasovno zahtevnost s polinomsko interpolacijo.

Primer 3.5. Poglejmo si absolutno napako pri raˇcunanju vrednosti integrala (6), I =

∫︂ 1

−1

cos(x)e^iωxdx,

kjer funkcijo f(x) = cos(x) aproksimiramo s polinomom stopnje n, ki ga dobimo z enim korakom Remesovega postopka za n = 1,2, . . . ,5. Pri raˇcunanju aproksima- cije integrala si zopet pomagamo z vrednostmi momentov M_j(x), j = 0, 1, . . . , n+ 1.

Na sliki 8 vidimo, da se z naraˇsˇcanjem ω in veˇcanjem stopnje f˜︁absolutna napaka zmanjˇsuje. Zaradi sodosti funkcije cos(x) pride do pojava, kjer je absolutna napaka pri raˇcunanju integrala, kjer f aproksimiramo z f˜︁sode stopnje (n = 2k), po na- tanˇcnosti primerljiva z absolutno napako, kjer namestof uporabimo f, ki ima eno˜︁

stopnjo veˇc (n = 2k+ 1). ♢

3.4. Metoda najmanjˇsih kvadratov. Pri aproksimaciji funkcije po metodi naj- manjˇsih kvadratov potrebujemo drugo normo, ki jo definiramo s pomoˇcjo skalarnega produkta na nekem intervalu [a, b]:

∥f∥₂ =√︁

⟨f, f⟩.

Obiˇcajno izberemo ⟨f, g⟩ = ∫︁b

a f(x)g(x)dx, v naˇsih primerih pa bomo uporabili diskretni (semi)skalarni produkt

⟨f, g⟩=

N

∑︂

i=0

f(τ_i)g(τ_i), (τ_i)^N_i=0 ∈[a, b],

pri ˇcemer bomo za τ0, τ1, . . . , τN izbrali ekvidistantne toˇcke. Za izbrano funkcijo f iˇsˇcemo f˜︁∈S = lin{ϕ₀, ϕ₁, . . . , ϕ_n}, da velja

⃦

⃦

⃦f −f˜︁

⃦

⃦

⃦2

= min

s∈S∥f −s∥₂.

Aproksimant f˜︁imenujemo element najboljˇse aproksimacije po metodi najmanjˇsih kvadratov.

Izrek 3.6. Funkcija f˜︁∈S je element najboljˇse aproksimacije po metodi najmanjˇsih kvadratov za f natanko tedaj, ko je f −f˜︁⊥S.

Slika 8. Absolutna napaka pri aproksimaciji vrednosti integrala (6) v odvisnosti od ω ∈ N,1 ≤ ω ≤ 100, kjer aproksimant f˜︁doloˇcimo z enim korakom Remesovega postopka.

Dokaz se nahaja v viru [1]. Iz izreka (3.6) pa sledi konstrukcija. Da bo f˜︁−f⊥S, mora biti funkcija f − f˜︁pravokotna na ϕ_i za vse i = 0, . . . , n. Ker je f˜︁ ∈ S, jo zapiˇsemo kot f˜︁ = ∑︁n

j=0c_jϕ_j, pogoj pravokotnosti pa je ekvivalenten pogoju

⟨f −f , ϕ˜︁ _i⟩= 0. Sledi:

⟨f−

n

∑︂

j=0

c_jϕ_j, ϕ_i⟩= 0

⟨f, ϕ_i⟩ −

n

∑︂

j=0

c_j⟨ϕ_j, ϕ_i⟩= 0

n

∑︂

j=0

cj⟨ϕj, ϕi⟩=⟨f, ϕi⟩.

Dobili smo sistem linearnih enaˇcb, ki se v matriˇcni obliki zapiˇse kot

⎡

⎢

⎢

⎢

⎣

⟨ϕ₀, ϕ₀⟩ ⟨ϕ₁, ϕ₀⟩ . . . ⟨ϕ_n, ϕ₀⟩

⟨ϕ₀, ϕ₁⟩ ⟨ϕ₁, ϕ₁⟩ . . . ⟨ϕ_n, ϕ₁⟩

... . .. ...

⟨ϕ₀, ϕ_n⟩ ⟨ϕ₁, ϕ_n⟩ . . . ⟨ϕ_n, ϕ_n⟩

⎤

⎥

⎥

⎥

⎦

⎡

⎢

⎢

⎢

⎣ c₀ c₁ ... c_n

⎤

⎥

⎥

⎥

⎦

=

⎡

⎢

⎢

⎢

⎣

⟨f, ϕ₁⟩

⟨f, ϕ₂⟩ ...

⟨f, ϕ_n⟩

⎤

⎥

⎥

⎥

⎦ .

Primer 3.7. Zopet raˇcunamo vrednost integrala (6), I =

∫︂ 1

−1

cos(x)e^iωxdx,

le da tokrat funkcijof(x) = cos(x) aproksimiramo po metodi najmanjˇsih kvadratov, kjer uporabimo 100 ekvidistantnih toˇck v diskretnem skalarnem produktu. Apro- ksimant iˇsˇcemo v podprostoru polinomovPk za razliˇcne stopnje k = 1,2, . . . ,6. Za bazne polinome izberemo kar potence, to je ϕ_j(x) = x^j, j = 0,1, . . . , k. Polinom najboljˇse aproksimacije po metodi najmanjˇsih kvadratov p˜︁_k ∈ Pk, p˜︁_k = ∑︁k

j=0c_jx^j je tako doloˇcen z reˇsitvijo sistema linearnih enaˇcb

⎡

⎢

⎢

⎢

⎣

⟨1,1⟩ ⟨x,1⟩ . . . ⟨x^k,1⟩

⟨1, x⟩ ⟨x, x⟩ . . . ⟨x^k, x⟩

... ... . .. ...

⟨1, x^k⟩ ⟨x, x^k⟩ . . . ⟨x^k, x^k⟩

⎤

⎥

⎥

⎥

⎦

⎡

⎢

⎢

⎢

⎣ c₀ c₁ ... ck

⎤

⎥

⎥

⎥

⎦

=

⎡

⎢

⎢

⎢

⎣

⟨cos(x),1⟩

⟨cos(x), x⟩

...

⟨cos(x), x^k⟩

⎤

⎥

⎥

⎥

⎦

, k = 1,2, . . .6.

Pribliˇzek za integral dobimo preko formule (18). Na sliki 9 so prikazane absolutne napake pri aproksimaciji vrednosti integrala glede na stopnjo aproksimacijskega po- linoma p˜︁k ∈P^k, k = 0,1, . . . ,6. Vidimo, da je precejˇsnja razlika v absolutni napaki med aproksimacijskimi polinomi sode in lihe stopnje.

Slika 9. Absolutna napaka pri aproksimaciji vrednosti integrala (6), kjer f aproksimiramo s polinomi p˜︁_k ∈ P^k po metodi najmanjˇsih kva- dratov.

♢

4. Zakljuˇcek

Spoznali smo asimptotsko metodo za integrale oblike (1), kjer funkcija g nima stacionarnih toˇck na intervalu integracije in njeno nadgradnjo v primerih, ko ima eno ali veˇc stacionarnih toˇck. Ker moramo veˇckrat v praksi funkcijo f spremeniti oziroma poznamo le nekaj njenih vrednosti, smo si pogledali tri metode za aproksi- macijo funkcije f, ki nam ob primerni funkciji g in dobro izbrani bazi za aproksi- macijo olajˇsajo delo in zmanjˇsajo integracijsko napako. Vse predstavljene metode smo implementirali in preizkusili na numeriˇcnih primerih. Asimptotske metode so se izkazale kot najboljˇse. Med Filonovimi metodami sta interpolacija in Remesov postopek vrnila primerljive rezultate pri primerljivi ˇcasovni zahtevnosti, medtem ko so bili rezultati, pridobljeni z aproksimacijsko funkcijo po metodi najmanjˇsih kvadratov, najslabˇsi.

Slovar strokovnih izrazov asymptotic methods asimptotske metode

Filon methods Filonove metode

highly oscillatory function hitro oscilirajoˇca funkcija least squares method metoda najmanjˇsih kvadratov polynomial interpolation polinomska interpolacija Remez algorithm Remesov postopek

Literatura

[1] J. Kozak,Numeriˇcna analiza, DMFA – zaloˇzniˇstvo, Ljubljana, 2008, 417 str.

[2] F.W.J. Olver, Error Bounds for Stationary Phase Approximations, SIMA 5 (1974), strani 19–29.

[3] E. Stein in T. Murphy,Harmonic Analysis: Real-Variable Methods, Orthogonality, and Oscil- latory Integrals, PMS–43, Princeton University Press, Princeton, 1993.

[4] J. Trivedi,A Survey Of Numerical Quadrature Methods For Highly Oscillatory Integrals, ma- gistrsko delo, Graduate & Postdoctoral Studies, Western University, 2019.

[5] Integral calculator, [ogled 13. 8. 2021], dostopno na https://www.integral-calculator.

com.