UNIVERZA V LJUBLJANI
FAKULTETA ZA MATEMATIKO IN FIZIKO Finanˇ cna matematika – 1. stopnja
Anja Trobec
Newtonovi interpolacijski polinomi v veˇ c spremenljivkah
Delo diplomskega seminarja
Mentor: prof. dr. Marjetka Knez
Ljubljana, 2021
Kazalo
1. Uvod 5
1.1. Posploˇsitev polinomov ene spremenljivke 6
2. Newtonov interpolacijski polinom v eni spremenljivki 6
2.1. Deljene diference 7
3. Interpolacija na mreˇzi 9
3.1. Interpolacijske toˇcke 9
3.2. Multiindeksna notacija 10
3.3. Strukturirana indeksna mnoˇzica in deljene diference 11 3.4. Primer interpolacije funkcije dveh spremenljivk 12 3.5. Algoritem za raˇcunanje deljenih diferenc v dveh spremenljivkah 16 4. Obstoj in enoliˇcnost interpolacijskega polinoma na mreˇzi 16 5. Formula deljenih diferenc v Rd za mreˇzne interpolacijske toˇcke 18 6. Interpolacijski polinomi na poljubno izbranih interpolacijskih toˇckah v
ravnini 20
7. Interpolacijski prostor minimalne stopnje 24
8. Zakljuˇcek 25
Slike 27
Tabele 27
Literatura 28
Newtonovi interpolacijski polinomi v veˇc spremenljivkah Povzetek
Algoritme in tehnike reˇsevanja problema interpolacije v eni spremenljivki lahko razˇsirimo na reˇsevanje v veˇc spremenljivkah z razliˇcnimi posploˇsitvami in nadgra- dnjami le-teh. Obliko Newtonove baze in algoritem deljenih diferenc za iskanje pripadajoˇcih koeficientov lahko neposredno posploˇsimo za interpolacijo na mreˇznih toˇckah. Ta posploˇsitev je tenzorska ali pa v obliki omejitve skupne stopnje. Pri tej posploˇsitvi uporabljamo multiindeksno notacijo za sklicevanje na interpolacijske toˇcke. Mnoˇzicam interpolacijskih toˇck, za katere se uporablja tenzorski prostor, ime- nujemo polne mnoˇzice. Sem sodijo t. i. ˇskatlaste mnoˇzice toˇck in trikotne mnoˇzice toˇck. Poljubno izbrane interpolacijske toˇcke ne zagotavljajo enoliˇcne interpolacije, vendar v nekaterih primerih s posploˇsitvijo baze, ki jo imenujemo Newton-Sauerjeva baza, pridemo do preprostega trikotnega linearnega sistema, ki vrne ustrezne koe- ficiente interpolacijskega polinoma. Za ustrezno ˇstevilo paroma razliˇcnih interpo- lacijskih toˇck nam algoritem, ki temelji na Gaussovih eliminacijah, vrne Newton- Sauerjevo bazo ali pa polinom, ki ima vrednost 0 na vseh interpolacijskih toˇckah, kar nam pove, da enoliˇcna interpolacija ni mogoˇca. Algoritem lahko uporabimo tudi, da za dan nabor interpolacijskih podatkov konstruiramo polinomski podprostor mi- nimalne stopnje, kjer je enoliˇcna interpolacija vedno mogoˇca.
Multivariate Newton interpolation polynomial Abstract
Techniques and algorithms for finding univariate Newton interpolating polynomi- als can be extended to multivariate data points by different generalizations. The Newton basis format and divided-difference algorithm for coefficients can be genera- lized in a straightforward way when interpolating at nodes on a grid. Two different approaches, the tensor product case or the triangular case, are usually considered.
The multi-index notation is used to refer to nodes. Node configurations where the tensor product case applies are called lower sets, and they include n-dimensional rectangles and triangles. Arbitrary distinct nodes do not ensure unique interpola- ting polynomials but, when possible, a different basis generalization, which we call Newton–Sauer basis, results in a nice triangular linear system for the coefficients.
For the right number of distinct nodes, an algorithm based on Gaussian elimination produces either this Newton–Sauer basis or a polynomial that is zero on all nodes, showing that unique interpolation is impossible. The algorithm may also be used on any distinct nodes to produce a polynomial subspace of minimal degree where unique interpolation is possible.
Mathematics subject classifications (2020): 41A10 , 41A05, 41A63, 65D05 Kljuˇcne besede: interpolacija, interpolacijske toˇcke, Newtonovi polinomi, deljene diference, algoritem
Keywords: interpolation, nodes, Newton polynomial, divided difference, algorithm
1. Uvod
Namen interpolacije je najti ˇcim bolj preprosto funkcijo, ki aproksimira neko dano funkcijo. Najveˇckrat za interpolacijsko funkcijo vzamemo kar polinom, ki se v iz- branih toˇckah domene ujema z vrednostmi interpolirane funkcije. Izbranim toˇckam pravimo interpolacijske toˇcke.
Problem interpolacije v eni spremenljivki lahko s pomoˇcjo osnovnega znanja nu- meriˇcnih metod in linearne algebre razˇsirimo na interpolacijo funkcij veˇc spremen- ljivk. Cilj razˇsiritve je poiskati konstruktivne algoritme, s katerimi sestavimo ustre- zen interpolacijski polinom.
Vpraˇsanje na mestu je, ali interpolacijski polinom obstaja in ˇse veˇc, ali je enoliˇcno doloˇcen. V primeru interpolacije funkcije ene spremenljivke vemo, da nam izbira n + 1 paroma razliˇcnih interpolacijskih toˇck zagotavlja enoliˇcen obstoj interpola- cijskega polinoma stopnje n. V primeru veˇc spremenljivk pa sta tako obstoj kot enoliˇcnost odvisna od ˇstevila in lokacije izbranih interpolacijskih toˇck.
Recimo, da iˇsˇcemo interpolacijski polinom v dveh spremenljivkah oblike (1) p(x, y) =a+bx+cy+dx2+exy+f y2.
Imamo ˇsest neznanih koeficientov, zato potrebujemo ˇsest enaˇcb (po moˇznosti linearnih), ki bodo te koeficiente enoliˇcno doloˇcale. Enaˇcbe lahko dobimo iz interpolacijskih pogojev na ˇsestih razliˇcnih toˇckah v ravnini. Te lahko v ravnini izberemo na poljubno mnogo naˇcinov, vendar se bomo v nadaljevanju osredotoˇcili na tri najbolj pogoste primere izbir, ki so prikazani na sliki 1.
-2 -1 0 1 2
-2 -1 0 1 2
0 1 2
0 1 2
-2 -1 0 1 2
-2 -1 0 1 2
Slika 1. Primeri izbire interpolacijskih toˇck: (a), (b) in (c).
Ko izbrane toˇcke vstavimo v polinom (1) in mu doloˇcimo vrednosti v teh toˇckah, dobimo sistem linearnih enaˇcb za neznane koeficiente.
V primeru (a) leˇzijo izbrane interpolacijske toˇcke na kroˇznici s srediˇsˇcem v izhodiˇsˇcu in radijem 2. Od tod vidimo, da iz vsakega interpolacijskega polinoma oblike (1), ˇ
ce le ta obstaja, s priˇstevanjem ˇclena c(x2+y2−4), za poljubenc∈R, spet dobimo interpolacijski polinom. To nam pove, da ima sistem bodisi neskonˇcno mnogo reˇsitev
ali pa nima nobene reˇsitve. V primeru (b) z uporabo enostavne Newtonove baze predstavljene v nadaljevanju, ki vkljuˇcuje polinome sestavljene iz faktorjev (x−xi) in (y−yi) dobimo preprost trikotni sistem linearnih enaˇcb, ki ima vselej enoliˇcno reˇsitev. V primeru (c), kjer interpolacijske toˇcke nimajo kakˇsne posebne strukture, pa bi se lahko zatekli h Gaussovim eliminacijam na matriki linearnega sistema, ki bi nam odgovorile na vpraˇsanje obstoja in enoliˇcnosti reˇsitve. V primeru takega nabora toˇck bomo v nadaljevanju predstavili algoritem, ki temelji na Gaussovih eliminacijah. Ta bo poiskal
• posploˇsitev klasiˇcne oblike Newtonovega polinoma kot v primeru (b) in enoliˇcno doloˇcil interpolacijski polinom na izbrani mnoˇzici toˇck, ali pa
• polinom, ki ima na vseh vozliˇsˇcih vrednost 0, kot v primeru (a).
1.1 Posploˇsitev polinomov ene spremenljivke.
V vektorskem prostoru zveznih funkcij ene spremenljivke izberemo podprostor polinomov stopnje ≤n in ga oznaˇcimo s Pn,
Pn = Lin{1, x, x2, ..., xn}=
⎧
⎨
⎩
n
∑︂
i=0
aixi;ai ∈R
⎫
⎬
⎭ .
Polinome v eni spremenljivki lahko posploˇsimo na veˇc spremenljivk na dva naˇcina:
(1) Tenzorski naˇcin - skupna stopnja (n, m)
Naj bodo polinomi spremenljivkex iz prostoraPn in polinomi spremenljivke yiz prostoraPm. Ko naredimo tenzorski produkt med prostoromaPn inPm, dobimo prostor polinomov skupne stopnje≤n+m,
P2n,m =Pn⊗Pm = Lin{xiyj;i= 0, . . . , n, j = 0, . . . , m}.
Eno od moˇznih baz tega podprostora torej sestavljajo produkti monomov.
(2) Polinomi skupne stopnje najveˇc n
Drugi naˇcin je omejitev skupne stopnje polinomov. Ce vzamemo po-ˇ linome skupne stopnje najveˇc 2, bo njihova baza sestavljena iz mono- mov 1, x, y, x2, xy, y2. Tovrstne polinome predstavimo v trikotni shemi. V sploˇsnem definiramo vektorski podprostor
P2n = Lin{xiyj; 0≤i+j ≤n}=
⎧
⎨
⎩
n
∑︂
i=0 n−i
∑︂
j=0
aixiyj;ai ∈R
⎫
⎬
⎭ .
2. Newtonov interpolacijski polinom v eni spremenljivki
Za razˇsiritev interpolacije na veˇc spremenljivk je potrebno najprej dobro razumeti problem v eni spremenljivki.
Vzemimo poljuben nabor paroma razliˇcnih toˇck x0, x1, . . . , xn in pripadajoˇce
vrednosti fi = f(xi), i = 0,1, . . . , n, izbrane funkcije f v teh toˇckah. Realne vrednostixi skupaj z vrednostmifi tvorijo interpolacijske toˇcke (xi, fi), skozi katere bo potekal graf interpolacijskega polinoma. Za vsak k ∈ {0,1, . . . , n} definiramo polinom oblike
qk(x) =
k−1
∏︂
i=0
(x−xi).
Mnoˇzica vseh takih polinomov tvori bazo prostora polinomov Pn, pri ˇcemer je qk(xj) = 0 za vsej < k.
Linearna kombinacija baznih polinomov qk predstavlja Newtonov interpolacij- ski polinom
(2) p(x) =
n
∑︂
i=0
aiqi(x),
katerega koeficiente doloˇcimo preko trikotnega sistema linearnih enaˇcbp(xj) = fj, j = 0,1, . . . , n. Za iskanje le-teh je bolj kot trikotna substitucija uˇcinkovit algoritem deljenih diferenc. (Natanˇcneje predstavljeno v [7]).
2.1 Deljene diference.
Definicija 2.1. Deljena diferenca [x0, x1, . . . , xk]f je vodilni koeficient interpola- cijskega polinoma stopnje k (koeficient pri potenci xk), ki se s funkcijo f ujema v paroma razliˇcnih toˇckah x0, x1, . . . , xk.
Izrek 2.2. Naj bodo x0, x1, . . . , xk paroma razliˇcne toˇcke in f dana funkcija.
(1) Za deljene diference velja naslednja rekurzivna formula:
(3) [x0, . . . , xk]f = [x1, . . . , xk]f −[x0, . . . , xk−1]f
xk−x0 .
(2) Newtonov interpolacijski polinom, ki se v toˇckah x0, x1, . . . , xn ujema s funk- cijo f, je oblike
pn(x) = [x0]f + [x0, x1]f(x−x0) + [x0, x1, x2]f(x−x0)(x−x1)+
+. . .+ [x0, x1, . . . , xn]f(x−x0)· · ·(x−xn−1).
Dokaz.
(1) Dokaˇzimo najprej rekurzivno formulo. Naj bo q0 interpolacijski polinom, ki se s funkcijo f ujema v toˇckah x0, x1, . . . , xk−1 in polinom q1 interpolacijski polinom, ki se s funkcijof ujema vx1, x2, . . . , xk. Hitro vidimo, da ima inter- polacijski polinomp, ki se ujema s funkcijof v vseh naˇstetih interpolacijskih toˇckah, obliko
p(x) = x−xk
x0−xkq0(x) + x−x0 xk−x0q1(x).
Ko primerjamo vodilne koeficiente, dobimo ravno formulo (3).
(2) Ta del dokaˇzemo z uporabo indukcije. Ker je po definiciji [x0]f =f(x0), je oˇcitno, da je polinom p0(x) = [x0]f interpolacijski za f v toˇcki x0.
Naj bo
pi(x) = [x0]f+ [x0, x1]f(x−x0) +. . .+ [x0, . . . , xi]f(x−x0)· · ·(x−xi−1) polinom, ki interpolira funkcijo f v toˇckahx0, . . . , xi. Potem je
pi+1(x) = pi(x) + [x0, x1, x2, . . . , xi, xi+1]f(x−x0)· · ·(x−xi)
po definiciji deljene diference enak interpolacijskemu polinomu, ki interpolira funkcijo f v toˇckah x0, x1, . . . , xi+1.
□ Deljene diference obiˇcajno zapisujemo kot [xi, xi+1, xi+2, . . . , xj]f , vendar bomo zaradi preglednosti v nadaljevanju uvedli krajˇso notacijo [i, j]f.
Ko v rekurzivno formulo vpeljemo novo notacijo, dobimo
(4) [i−1, j]f = [i, j]f−[i−1, j−1]f xj −xi−1
.
S podano zaˇcetno vrednostjo [j, j]f = fj raˇcunamo, dokler ne doseˇzemo [0, j]f, ki predstavlja ravno koeficient aj v formuli (2). Raˇcunanje deljenih diferenc najlaˇzje ponazorimo v trikotni shemi, pri ˇcemer izraˇcunane koeficiente preberemo iz zgornje diagonale, kot je prikazano v spodnji tabeli. (Povzeto po [3].)
[·]f [·,·]f [·,·,·]f . . . [·, . . . ,·]f x0 [0,0]f
[0,1]f
x1 [1,1]f [0,2]f
[1,2]f . ..
x2 [2,2]f
... ... ... [0, n]f
xn−1 [n−1, n−1]f ...
[n−1, n]f xn [n, n]f
3. Interpolacija na mreˇzi 3.1 Interpolacijske toˇcke.
V tem poglavju se osredotoˇcimo na interpolacijo na mreˇzi in s tem na mreˇzne interpolacijske toˇcke. Kombinacija dveh toˇck v R ustvari mreˇzno toˇcko v R2. V sploˇsnem bodo mreˇzne toˇcke iz prostora Rd, kjer ˇstevilo d oznaˇcuje ˇstevilo spremenljivk interpolacijskega polinoma.
Ze na zaˇˇ cetku je pomembno poudariti, da so vse interpolacijske toˇcke hkrati tudi mreˇzne toˇcke, vendar pa ne velja obratno, saj interpolacijska funkcija ne upoˇsteva nujno vseh mreˇznih toˇck.
Pri prehodu na veˇc spremenljivk je potrebno najprej razloˇziti notacijo. Mreˇzne toˇcke ˇzelimo urediti na naˇcin, da jih preindiksiramo. Novo indeksacijo doloˇcimo glede na ˇstevilo pojavitev istih vrednosti i-te koordinate za i = 1,2, . . . , d, od vseh interpolacijskih toˇck. Natanˇcneje, za vsako dimenzijo doloˇcimo nov vektor, v katerem zapiˇsemo vse razliˇcne vrednosti projekcij toˇck xi na izbrano koordinatno os. V vektor jih zapiˇsemo po padajoˇcem vrstnem redu ˇstevila pojavitev.
Za laˇzjo predstavo si poglejmo postopek na primeru dveh spremenljivk. V ravnini imamo podane interpolacijske toˇcke oblike xi = (xi, yi). Te najprej projeciramo na prvo koordinatno os - torejx-os in projekcije zdruˇzimo v vektorx= (x(0), x(1), . . .).
Z najmanjˇsim indeksom x(0) oznaˇcimo tisto vrednost xi, na katero se projecira najveˇc mreˇznih toˇck in jo postavimo na prvo mesto naˇse nove ureditve. Enako nadaljujemo in urejamo po padajoˇcem ˇstevilu projekcij. Ko uredimo toˇcke nax osi, se na enak naˇcin lotimo urejanja na y osi.
Preizkusimo tovrstno ureditev na primeru iz slike 2. Najprej skuˇsamo urediti toˇcke na x osi. Pogledamo, v katero vrednost na x-osi se projecira najveˇc mreˇznih toˇck in vidimo, da se to zgodi za vrednost x = 0. Oznaˇcimo jo z najniˇzjo oznako x(0) in postavimo na prvo mesto v vektorju nove ureditve. Na enak naˇcin postopamo v naslednjih korakih. Naslednja vrednost je x = 1, oznaˇcimo jo z x(1). Nato x=−1 oznaˇcimo zx(2) in nazadnje z eno samo projecirano toˇcko vrednostix= 0.5 pripiˇsemo novo oznako x(3).
Na ta naˇcin smo dobili vektor
x= (x(0), x(1), x(2), x(3)) = (0,1,−1,0.5).
Podobno storimo na drugi koordinatni osi in dobimo vektor y= (y(0), y(1), y(2), y(3)) = (1,−1,0,−0.5).
-1.5 -1 -0.5 0 0.5 1 1.5 -1.5
-1 -0.5 0 0.5 1 1.5
Slika 2. Interpolacija na mreˇzi.
3.2 Multiindeksna notacija.
Naslednja pomembna stvar, ki jo uvedemo, je multiindeksna notacija. Prej izpeljane vektorje, ki so sestavljeni iz projekcij interpolacijskih toˇck na koordinatne osi, bomo opremili z multiindeksi.
Definiramo multiindeks λ ∈ Nd0 kot vektor oblike λ= (λ1, λ2, . . . , λd), pri ˇcemer je λi ∈ {0,1,2,3, ...}. ˇCe multiindeks λ ustreza mreˇzni toˇcki x(λ), lahko definiramo sklic na toˇcko s pomoˇcjo multiindeksa kot x(λ) = (x1(λ1), x2(λ2), . . . , xd(λd)). Red multiindeksa izraˇcunamo kot |λ|=λ1+λ2 +. . .+λd.
Na primeru iz slike 2 bi multiindeks λ = (2,1) ustrezal mreˇzni toˇcki x((2,1)) = (x(2), y(1)) = (−1,−1) in λ = (0,2) mreˇzni toˇcki x((0,2)) = (x(0), y(2)) = (0,0).
Ce bi se sklicevali na toˇˇ cko iz prostora R3, bi za λ(1,3,0) doloˇcili x(λ) = x((1,3,0)) = (x(1), y(3), z(0)). Pripomnimo ˇse, da je red multiindeksa λ = (2,0) enak 2, multiindeksa λ= (2,1) 3, multiindeksaλ= (1,3,0) pa 4.
Notacija λ < µ, pomeni, da stroga neenakost velja v vsaki koordinati, torej λi < µi. Konˇcno za vsak multiindeks λ, ki se ujema s kakˇsno toˇcko x(λ), definiramo polinom
(5) qλ(x) =
d
∏︂
m=1 λm−1
∏︂
i=0
(xm−xm(i)),
kjer sox= (x1, . . . , xd) neodvisne spremenljivke in vsakxm(i),m= 1, . . . , d,preteˇce vse vrednosti med 0≤i < λm. V primeru iz slike 2 za multiindeksλ= (2,1) dobimo
q(2,1)(x, y) = (x−x(0))(x−x(1))(y−y(0)) =x(x−1)(y−1).
Definiramo ˇse indeksno mnoˇzico I ⊂ Nd, v katero shranimo vse multiindekse in jih tako loˇcimo od obiˇcajnih mreˇznih toˇck. Nato poiˇsˇcemo Newtonov interpolacijski polinom oblike p(x) = ∑︁
λ∈Iaλqλ(x), ki funkcijo f interpolira na izbranih toˇckah {(x(λ), fλ) :λ∈I}, I ⊂Nd.
3.3 Strukturirana indeksna mnoˇzica in deljene diference.
Strukturirana indeksna mnoˇzica I nam omogoˇca posploˇsitev formule deljenih diferenc na
(6) [α−em,β]f = [α,β]f−[α−em,β−em]f xm(βm)−xm(αm−1) ,
za α≤β z αm >0, pri ˇcemer je em = (0,0, . . . ,1, . . . ,0) enotski vektor z enico na m-ti koordinati.
Definicija 3.1. Konˇcni nabor multiindeksovJ ⊆Nd0 imenujemo polna mnoˇzica, ˇce iz
β ∈J ∧ 0≤α≤β =⇒α∈J.
Velja tudi, da bo multiindeks αpred multiindeksom β, ˇce bo veljalo
(7) (|α| ≤ |β|) ali (|α|=|β| inαk> βk medtem, ko za j < k veljaαj =βj).
Pomembna posebna primera polnih mnoˇzic multiindeksov sta (1) ˇskatlasta mnoˇzica Jβ ={α:0≤α≤β}, za fiksen β, (2) trikotna mnoˇzica Jn ={α:|α| ≤n}, za fiksno stopnjon.
Trikotno mnoˇzico toˇck prepoznamo na primeru (b) slike 1. ˇCe velja |α| ≤2, torej imamo fiksno stopnjo n = 2, to pomeni, da je naˇsa polna mnoˇzica sestavljena iz toˇck J2 = {(0,0),(0,1),(0,2),(1,0),(1,1),(2,0)}. Te toˇcke predstavljajo ravno trikotno shemo na mreˇzi (glej sliko 3).
0 1 2
0 1 2
Slika 3. Trikotna shema toˇck.
Ce pa namesto omejitve skupne stopnje doloˇˇ cimo omejitev posamezne stopnje z nekim fiksnim multiindeksomβ, dobimo polno mreˇzo toˇck. Zaβ= (2,2) je mnoˇzica Jβ = {(0,0),(0,1),(0,2),(1,0),(1,1),(1,2),(2,0),(2,1),(2,2)} in je prikazana na sliki 4.
0 1 2
0 1 2
Slika 4. Polna shema toˇck.
Ce povzamemo, interpolacijski prostor, ki ga razpenjajo polinomiˇ qλ na trikotni mnoˇzici toˇck Jn za fiksno stopnjo n, je prostor polinomov Pdn, medtem ko klasiˇcen tenzorski prostor Lin{xα11xα22. . . xαdd :0≤α≤β}dobimo s kombinacijoqλrazpetih na ˇskatlasti mnoˇzici, ki jo omejujeβ.
V nadaljevanju nas bo zanimalo, kako poiˇsˇcemo koeficiente interpolacijskega polinoma s pomoˇcjo formule deljenih diferenc.
Vzamemo mnoˇzico toˇck {(x(λ), fλ) : λ ∈ J}, pri ˇcemer je J polna mnoˇzica.
V formuli (6) kot dano zaˇcetno vrednost uporabimo [λ,λ]f = fλ, ki nas preko rekurzivne zveze pripelje do koeficienta [0,λ]f =aλ Newtonovega interpolacijskega polinoma
(8) p(x) = ∑︂
λ∈J
aλqλ(x).
Vpraˇsanje je, v kakˇsnem zaporedju moramo odˇstevati enotske vektorje in v kakˇsnem vrstnem redu izbiramo elemente iz mnoˇzice J. Uvesti je treba ustrezno pravilo.
Za primer vzemimo β = (1,2,1). Da iz deljene diference [β,β]f pridemo do [0,β]f, moramo odˇsteti multiindekse e1,e2,e2 in e3 in sicer v smeri od leve proti desni. V rekurzivni formuli tako potrebujemo ˇse deljene diference [α,β]f za sledeˇce α:
(1,2,1),(0,2,1),(0,1,1),(0,0,1) in (0,0,0).
Lahko bi seveda uporabili tudi kakˇsen drug vrstni red, vendar bi pri tem dobili drugaˇcne vrednosti v vmesnih korakih algoritma deljenih diferenc.
3.4 Primer interpolacije funkcije dveh spremenljivk.
Uporabili bomo interpolacijske toˇcke iz slike 2, kjer za urejena vektorja vza- memo
x= (x(0), x(1), x(2), x(3)) = (0,1,−1,0.5) in
y= (y(0), y(1), y(2), y(3)) = (1,−1,0,−0.5).
V izbranem primeru je stopnja n = 3 in ˇstevilo spremenljivk d = 2. Po ureditvi multiindeksov dobimo ravno trikotno shemo toˇck. Iˇsˇcemo Newtonov interpolacijski polinom oblike
p(x, y) = ∑︁
0≤λ1+λ2≤3a(λ1,λ2)q(λ1,λ2)(x, y).
Zelene funkcijske vrednosti so predstavljene v matrikiˇ A:
A=
⎡
⎢
⎢
⎣ 2.75
3 3
9 10 16
5 8 2 4.25
⎤
⎥
⎥
⎦ .
Cez matriko bomo ˇsli vˇ k korakih, pri ˇcemer k teˇce od 1 do n. Vrednosti bomo raˇcunali po diagonalah navzdol in na desni strani matrike s puˇsˇcicami prikazali posodobitev matrike na naslednjem koraku. Algoritem je prikazan v Tabeli 1.
Za primer lahko pogledamo multiindeks β = (2,1), ki nam pove, da bomo odˇstevali enotske vektorje e1,e1 in e2. Odˇstevanje vektorja e1 bomo nakazali s puˇsˇcico levo, saj gre za spremembo v x koordinati, medtem ko bo odˇstevanje e2 prikazano s puˇsˇcico navzdol, saj gre za spremembno v y koordinati.
Poleg smeri puˇsˇcice bomo vsaki puˇsˇcici pribeleˇzili ˇstevec t, ki nam bo pove- dal, za koliko korakov je puˇsˇcica na tem mestu ˇze obrnjena v prikazano smer. V resnici bo to pomenilo, kako narazen sta vrednosti, ki ju odˇstejemo v imenovalcu (glej formulo (6)).
Tabela 1. Algoritem deljenih diferenc
λ2 y(λ2) (a) funkcijske vrednosti funkcijska shema λ2 1. korak
3 -0.5 2.75 [(0,3),(0,3)]f 3 ↓1
2 0 3 3 [(0,2),(0,2)]f [(1,2),(1,2)]f 2 ↓1 ←1
1 -1 9 10 16 [(0,1),(0,1)]f [(1,1),(1,1)]f [(2,1),(2,1)]f 1 ↓1 ←1 ←1
0 1 5 8 2 4.25 [(0,0),(0,0)]f [(1,0),(1,0)]f [(2,0),(2,0)]f [(3,0),(3,0)]f 0 ←1 ←1 ←1 x(λ1) : 0 1 -1 0.5
λ1: 0 1 2 3 λ1: 0 1 2 3
λ2 y(λ2) (b) po prvem koraku funkcijska shema λ2 2. korak
3 -0.5 0.5 [(0,2),(0,3)]f 3 ↓2
2 0 -6 0 [(0,1),(0,2)]f [(0,2),(1,2)]f 2 ↓2 ↓1
1 -1 -2 1 -3 [(0,0),(0,1)]f [(0,1),(1,1)]f [(1,1),(2,1)]f 1 ↓1 ←2
0 1 5 3 3 1.5 [(0,0),(0,0)]f [(0,0),(1,0)]f [(1,0),(2,0)]f [(2,0),(3,0)]f 0 ←2 ←2
x(λ1) : 0 1 -1 0.5
λ1: 0 1 2 3 λ1: 0 1 2 3
λ2 y(λ2) (c) po drugem koraku funkcijska shema λ2 3. korak
3 -0.5 13 [(0,1),(0,3)]f 3 ↓3
2 0 4 -1 [(0,0),(0,2)]f [(0,1),(1,2)]f 2 ↓2
1 -1 -2 1 4 [(0,0),(0,1)]f [(0,0),(1,1)]f [(0,1),(2,1)]f 1 ↓1
0 1 5 3 0 3 [(0,0),(0,0)]f [(0,0),(1,0)]f [(0,0),(2,0)]f [(1,0),(3,0)]f 0 ←3
x(λ1) : 0 1 -1 0.5
λ1: 0 1 2 3 λ1: 0 1 2 3
λ2 y(λ2) (d) po tretjem koraku funkcijska shema
3 -0.5 -6 [(0,0),(0,3)]f
2 0 4 2 [(0,0),(0,2)]f [(0,0),(1,2)]f
1 -1 -2 1 -2 [(0,0),(0,1)]f [(0,0),(1,1)]f [(0,0),(2,1)]f
0 1 5 3 0 6 [(0,0),(0,0)]f [(0,0),(1,0)]f [(0,0),(2,0)]f [(0,0),(3,0)]f x(λ1) : 0 1 -1 0.5
λ1: 0 1 2 3
Ce si pogledamo primer na drugem koraku za multiindeks (2,ˇ 1), ko puˇsˇcica kaˇze
←2, to prikazuje operacijo
A[(2,1)]← A[(2,1)]−A[(1,1)]
x(2)−x(2−2) = −3−1
−1−0 = 4, kar je ravno vrednost [(0,1),(2,1)]f.
Ko posodabljamo vrednosti matrike A, moramo paziti na vrstni red raˇcunanja.
Sledimo pravilu, da mora biti operacija, h kateri kaˇze puˇsˇcica, izraˇcunana pred tisto na njeni levi.
Vpraˇsanje je le ˇse, kako iz konˇcne matrike A preberemo koeficiente interpola- cijskega polinoma, tj. deljene diference [0,β]f. Koeficiente preberemo iz konˇcne tabele po diagonalah navzgor, kot je prikazano v zadnji funkcijski shemi v Ta- beli 1. Po prvem koraku dobimo koeficiente linearnega polinoma. Po drugem koraku lahko iz druge diagonale preberemo koeficiente kvadratnega polinoma in po zadnjem, tretjem koraku, preberemo ˇse koeficiente kubiˇcnega polinoma.
Tako dobimo rekurzivno zvezo med interpolacijskimi polinomi, katerih koeficiente preberemo iz konˇcne oblike tabele glede na relacijo urejenosti med multiindeksi.
Glede na obratno leksikografsko multiindeksno ureditev iz (7) je v spodnji ta- beli doloˇcen vrstni red, po katerem izbiramo koeficiente iz zadnjega koraka tabele 1:
10 6 9 3 5 8 1 2 4 7 Dobljeni polinomi so
p1(x, y) = 5 + 3x−2(y−1),
p2(x, y) = p1(x, y) + 0x(x−1) + 1x(y−1) + 4(y−1)(y+ 1),
p3(x, y) = p2(x, y) + 6x(x−1)(x+ 1)−2x(x−1)(y−1) + 2x(y−1)(y+ 1)
−6(y−1)(y+ 1)y, pri ˇcemer
• graf p1(x, y) = 7 + 3x−2y predstavlja ravnino skozi toˇcke (0,1,5),(1,1,8) in (0,−1,9),
• graf p2(x, y) = 3 + 2x−2y+xy+ 4y2 predstavlja hiperboliˇcni paraboloid skozi toˇcke (0,1,5),(1,1,8),(0,−1,9),(−1,1,2), (1,−1,10) in (0,0,3),
• zadnji kubiˇcni polinom pa gre skozi vseh 10 interpolacijskih toˇck in se poe- nostavi v p3(x, y) = 3−8x+ 4y+ 2x2+ 3xy+ 2y2+ 6x3−2x2y+ 2xy2−6y3. Polinomi p1, p2 in p3 so prikazani na slikah 5, 6, 7. Poleg tega lahko iz tabele preberemo tudi druge interpolacijske polinome, na primer poli- nome iz tenzorskega prostora. Za multiindeks β = (2,1) je reˇsitev pro- blema Newtonov interpolacijski polinom na ˇsestih interpolacijskih toˇckah (0,1,5),(1,1,8),(−1,1,2),(0,−1,9),(1,−1,10),(−1,−1,16):
p(x, y) = 5 + 3x+ 0x(x−1)−2(y−1) +x(y−1)−2x(x−1)(y−1)
= 7−2y+ 3xy+ 2x2−2x2y.
Slika 5. Polinom p1.
Slika 6. Polinom p2.
Slika 7. Polinom p3.
Dobili smo enoliˇcno doloˇcen interpolacijski polinom p ∈ P22,1. Podobno je za multiindeks β= (1,2) reˇsitev problema Newtonov interpolacijski polinom na ˇsestih interpolacijskih toˇckah (0,1,5),(1,1,8),(0,−1,9),(1,−1,10),(0,2,3),(1,2,3):
p(x, y) = 5 + 3x−2(y−1) + 1x(y−1) + 4y(y−1)(y+ 1) + 2x(y−1)(y+ 1)
= 3−2y+xy+ 2xy2+ 4y2.
Iz sheme deljenih diferenc za zapis interpolacijskega polinoma p3 ∈ P23 na 10 razliˇcnih toˇckah lahko torej dobimo ˇse interpolacijske polinome iz podprostorov P22,1
in P21,2 na ustreznih mreˇznih toˇckah.
Formula deljenih diferenc, ki smo jo pri izraˇcunu uporabili, bo natanˇcneje opisana in dokazana v poglavju 5.
3.5 Algoritem za raˇcunanje deljenih diferenc v dveh spremenljivkah.
Algoritem, ki nas pripelje do zgornjega rezultata, lahko posploˇsimo na nasle- dnji naˇcin:
Algoritem: Sploˇsen algoritem deljenih diferenc
Vhodni podatki: stopnja n, matrika funkcijskih vrednosti A, trikotna mnoˇzica Jn
for k from 1 to n for j from n to k
for allλ∈J with |λ|=j if k−λ1 >0
t=k−λ1, m= 2 else
t=k, m= 1
A[λ] = A[λ]−A[λ−em] xm(λm)−xm(λm−t) end for
end for end for
Izhod: koeficienti shranjeni v matriki A
4. Obstoj in enoliˇcnost interpolacijskega polinoma na mreˇzi V nadaljevanju bomo spoznali izrek, ki nam bo zagotovil obstoj enoliˇcnega interpolacijskega polinoma ne le na trikotni ali ˇskatlasti mnoˇzici toˇck, paˇc pa na poljubni unisolventni mnoˇzici mreˇznih toˇck. Unisolventna mnoˇzica je tista mnoˇzica, na kateri obstaja enoliˇcno doloˇcen interpolacijski polinom. ˇCe mnoˇzica ne bo unisolventna, bomo konstruirali enoliˇcen interpolacijski polinom v drugaˇce izbranem podprostoru.
Poljubno mnoˇzico mreˇznih interpolacijskih toˇck lahko opiˇsemo z {x(λ) : λ ∈ I}, pri ˇcemer je I ∈ Nd0. Za vsak multiindeks λ ∈ I formula (5) definira polinom qλ,
ki je sestavljen iz faktorjev oblike (xi − xi(j)) za vse j < λi. To velja tudi, ˇce mreˇzne toˇcke x(µ) z µ≤ λ morda niso vsebovane v mnoˇzici interpolacijskih toˇck.
Enostavno je videti, da velja
(9) qλ(x(µ)) = 0 ⇐⇒ ∃m ∈ {1, . . . , d} za katerega je µm < λm,
kar pripelje do trikotnega sistema linearnih enaˇcb, kot bomo videli v naslednjem dokazu.
Izrek 4.1 (Sploˇsen obstoj in enoliˇcnost). Naj bo I ⊆Nd0 podmnoˇzica multiindeksov, ki ustrezajo mreˇznim toˇckam {x(λ) : λ ∈ I}. Potem so polinomi iz mnoˇzice {qλ : λ∈I}linearno neodvisni in lahko poiˇsˇcemo enoliˇcno doloˇcen interpolacijski polinom sestavljen iz linearnih kombinacij elementov mnoˇzice Lin{qλ:λ∈I}.
Dokaz. Naj bo p=∑︁
λ∈Iaλqλ ≡0. Enaˇcbe p(x(µ)) =∑︂
λ∈I
aλqλ(x(µ)) = 0 za vse µ∈I
lahko zapiˇsemo v matriˇcni obliki z uporabo pravila linearne urejenosti multiindeksov iz formule (8). Definiramo kvadratno matriko G = [qλ(x(µ))], pri ˇcemer indeks µ teˇce po vrsticah in indeksλ po stolpcih. Dobimo sistem linearnih enaˇcbG[aλ] =0, pri ˇcemer sta [aλ] in 0stolpiˇcna vektorja indeksirana z λ. Dokaˇzimo, da je matrika G spodnje trikotna. Vzemimo, da indeks λ preseˇze indeks µ, |µ| ≤ |λ| in µ̸= λ.
Potem lahko najdemo µm < λm in po toˇcki (9) velja qλ(x(µ)) = 0. Pri tem so diagonalni elementi qµ(x(µ))̸= 0, kar tudi sledi iz formule (9). Torej vidimo, da je matrikaGspodnje trikotna in obrnljiva, iz ˇcesar sledi, da je [aλ] =0in{qλ:λ∈I}
je linearno neodvisna mnoˇzica. Interpolacijski polinom ustreza enaˇcbam p(x(µ)) =
∑︁
λ∈Iaλqλ(x(µ)) =fµ za vse µ∈I, torej je [aλ] =G−1[fµ]. □ Iz dokaza vidimo, da ima matrikaGbloˇcno strukturo, kar bolj poenostavi raˇcunanje, kot ˇce ima le obrnljivo trikotno. Za |µ| ≤ |λ| velja qλ(x(µ)) ̸= 0 ⇐⇒ µ = λ.
Torej je vsak diagonalni blok z multiindeksi istega reda diagonalna matrika z neniˇcelnimi diagonalnimi elementi.
Newtonova baza, ki jo bomo sestavili, bo vsebovala polinome qλ s spodaj naˇstetimi lastnostmi. Iˇsˇcemo tak polinom qλ, ki
(1) ima vrednost 0 na vseh interpolacijskih toˇckah, ki ustrezajo multiindeksom niˇzje stopnje,
(2) ima neniˇcelno vrednost v toˇckix(λ) in niˇcelno na vseh ostalih interpolacijskih toˇckah, ki ustrezajo multiindeksom iste stopnje,
(3) je natanko stopnje |λ|.
Posledica 4.2. Za vsako mnoˇzico toˇck {(x(λ), fλ) : α≤ λ ≤β} obstaja enoliˇcen interpolacijski polinom v podprostoru
(10) Lin{
d
∏︂
m=1 λm−1
∏︂
i=αm
(xm−xm(i)) :α≤λ≤β}= Lin{xµ11xµ22. . . xµnn :0≤µ≤β−α}.
Tenzorski prostor polinomov Pdβ−α, ki ga definiramo v (10), je direktna posledica izreka 4.1 z zaˇcetkom v α in je kljuˇcen za dokaz iskanja koeficientov s pomoˇcjo formule deljenih diferenc.
Poglejmo si primer na ˇsestih interpolacijskih toˇckah v ravnini, ki so prikazane na sliki 1 (a). Izbor toˇck na kroˇznici nam ne omogoˇca enoliˇcne interpolacije s kvadratiˇcnim polinomom. Lahko pa na izbrane toˇcke gledamo kot na mreˇzne toˇcke. ˇStiri toˇcke na robu kvadrata indeksiramo z x(0) = −1, x(1) = 1, y(0) = −√
3, y(1) = √ 3 preostali dve toˇcki na x-osi pa z y(2) = 0, x(2) = −2, x(3) = 2. Za dobljene mreˇzne toˇcke s kombinacijo multiindeksov dobimo Newtonove bazne polinome q00, q10, q01, q11 in zraven ˇse q22(x, y) = (x−x(0))(x−x(1))(y−y(0))(y−y(1)) in q32(x, y) = (x−x(2))q22(x, y).
Za vrednosti na teh toˇckah torej po izreku 4.1 obstaja enoliˇcen interpolacij- ski polinom v podprostoru Lin{1, x, y, xy, q22, q32}.
Dobili smo enoliˇcno interpolacijo, vendar smo pri tem uporabili polinome stopnje ≤5. V nadaljevanju bomo spoznali, kako pridemo do enoliˇcne interpolacije z uporabo najniˇzje moˇzne stopnje.
5. Formula deljenih diferenc v Rd za mreˇzne interpolacijske toˇcke V tem poglavju bomo pojasnili, kako algoritem deljenih diferenc posploˇsimo za raˇcunanje na funkcijah veˇc spremenljivk. Pri tem za interpolacijske toˇcke vselej vzamemo mreˇzne toˇcke.
Mnoˇzico interpolacijskih toˇck opiˇsemo z {(x(λ), fλ) : λ ∈ J}, kjer je J polna mnoˇzica multiindeksov. Iz posledice 4.2 sledi naslednja definicija.
Definicija 5.1. Za α,β ∈ J z α ≤ β definiramo [α,β]f kot koeficient pri ˇclenu najviˇsje stopnje interpolacijskega polinoma za interpolacijske toˇcke in pripadajoˇce vrednosti {(x(λ), fλ) :α≤λ≤β}.
Torej [α,β]f je enoliˇcno doloˇcen koeficient interpolacijskega polinoma pri monomu xµ11xµ22. . . xµdd, kjer je µ = β− α in je enak enoliˇcno doloˇcenemu koeficientu pri ˇ
clenu najviˇsje stopnje v Newtonovi obliki zapisa polinoma iz formule (8).
Izrek 5.2. Polinom p(x) =∑︁
λ∈J[0,λ]f qλ(x) zadoˇsˇca enaˇcbi p(x(λ)) =fλ za vse λ∈J.
Dokaz. Po izreku 4.1 obstaja enoliˇcno doloˇcen interpolacijski polinom oblike p
ˆ(x) =∑︁
λ∈Jaλqλ(x). Naj bo µ∈J. Pokazati moramo, da jeaµ= [0,µ]f. Za ta fiksen µ je
fµ =pˆ(x(µ)) = ∑︂
λ≤µ
aλqλ(x(µ)) + ∑︂
λ∈J,λ≰µ
aλqλ(x(µ)).
Po formuli (9) je druga vsota enaka 0. Prva vsota pa predstavlja ravno interpolacijski polinom na ˇskatlasti mnoˇzici {(x(λ), fλ) :0≤λ≤µ} inaλ= [0,µ]f po definiciji.
□ V naslednjem izreku je podana rekurzivna formula, ki jo uporabljamo v algoritmu za raˇcunanje deljenih diferenc.
Izrek 5.3. Za α,β∈J z α≤β in αm >0 velja formula [α−em,β]f = [α,β]f −[α−em,β−em]f
xm(βm)−xm(αm−1) , kjer je em = (0, . . . ,0,1,0, . . . ,0) enotski vektor.
Slika 8. Multiindeksi v dokazu izreka formule deljenih diferenc.
Dokaz. Zamislimo si enoliˇcno polinomsko interpolacijo na razliˇcnih podmnoˇzicah multiindeksov, kot je predstavljeno na sliki 8. Naj bo p+(x) interpolacijski polinom za mnoˇzico toˇck {(x(λ), fλ) : α ≤ λ ≤ β} in p-(x) interpolacijski polinom za mnoˇzico toˇck {(x(λ), fλ) :α−em ≤λ≤β−em}. Definiramo
(11) p(x) =p-(x) +
(︃ xm−xm(αm−1) xm(βm)−xm(αm−1)
)︃
(p+(x)−p-(x)).
Ker polinomap+(x) inp-(x) pripadata mnoˇzici Lin{xµ11xµ22. . . xµdd :0≤µ≤β−α}
in p(x) pripada Lin{xµ11xµ22. . . xµdd : 0 ≤ µ ≤ β −α+em}, lahko pokaˇzemo, da polinomp(x) interpolira funkcijo na uniji podmnoˇzic{(x(λ), fλ) :α−em ≤λ≤β}.
Prekrivajoˇci multiindeksi, α ≤ λ ≤ β−em, ustrezajo p+(x(λ))−p-(x(λ)) = 0, torej je
p(x(λ)) =p-(x(λ)) =fλ.
Neprekrivajoˇci indeksi, kjer je λm = αm−1, ustrezajo xm(λm)−xm(αm−1) = 0, torej je
p(x(λ)) =p-(x(λ)) =fλ.
Drugi neprekrivajoˇci indeksi z λm =βm, pa ustrezajo
p(x(λ)) =p-(x(λ)) + (1)(p+(x(λ))−p-(x(λ))) =fλ.
Zakljuˇcimo, da je p(x) interpolacijski polinom na uniji multiindeksnih podmnoˇzic in po definiciji je koeficient ob ˇclenu najviˇsje stopnje enak [α−em,β]f. Najviˇsja stopnja, ki jo dobimo pri xµ11xµ22. . . xµdd zaµ=β−α+em, se pojavi v formuli (11), kjer linearni faktor mnoˇzimo s ˇclenom p+(x)−p-(x). Koeficient pri ˇclenu najviˇsje stopnje polinomap+(x) je [α,β]f, prip-(x) pa [α−em,β−em]f. Ko to vstavimo v formulo in izrazimo koeficient poleg ˇclena najviˇsje stopnje na desni, dobimo ˇzeljeno
formulo deljenih diferenc (11). □
6. Interpolacijski polinomi na poljubno izbranih interpolacijskih toˇckah v ravnini
Vpraˇsanje, ki si ga postavimo v tem razdelku, je, ali lahko najdemo Newtonovo bazo stopnje ≤n za poljubno mnoˇzico interpolacijskih toˇck (︁n+d
d
)︁ v Rd.
S pomoˇcjo linearne algebre lahko najdemo pristop, ki nam pomaga konstrui- rati ustrezno bazo ali pa poiˇsˇce ustrezno algebraiˇcno hiperploskev, ki vsebuje izbrane interpolacijske toˇcke, na katerih vrne vrednost 0 in s tem pokaˇze, da enoliˇcna interpolacija ni mogoˇca.
Poglejmo si primer interpolacije v standardni monomski bazi. Pri dveh spre- menljivkah in stopnji polinoma n = 2 izberemo 6 poljubnih paroma razliˇcnih interpolacijskih toˇck v ravnini. Opravka imamo s polinomom oblike
(12) p(x, y) =c00+c10x+c01y+c20x2+c11xy+c02y2.
Ko v polinom (12) vstavimo izbrane interpolacijske toˇcke (xi, yi) in njihove funk- cijske vrednosti, dobimo sistem linearnih enaˇcb, ki ga lahko zapiˇsemo v matriˇcni obliki. Novo matriko imenujemo kolokacijska matrika in jo oznaˇcimo z M. ˇCe je M singularna, potem enoliˇcne reˇsitve ni mogoˇce najti. ˇCe pa je M nesingularna, torej obrnljiva, lahko s pomoˇcjo Gaussovih eliminacij naredimo LU razcep. Matriko M zapiˇsemo kot produkt zgornje trikotne matrike U in spodnje trikotne matrikeL z enicami po diagonali. Z LU razcepom tako problem oblike Mc = f prevedemo na reˇsevanje dveh trikotnih linearnih sistemov, Ld = f in Uc = d, preko katerih najdemo ustrezne koeficiente za poljubno podano mnoˇzico funkcijskih vrednosti na izbranih interpolacijskih toˇckah.
V primeru, da je matrika M nesingularna, lahko namesto LU razcepa upora- bimo podoben algoritem na matriki V =MT, ki nam bo prav tako vrnil koeficiente nove baze. Novo dobljena baza, ki jo je predstavil matematik Sauer [4], ima zelo podobne lastnosti kot Newtonova baza. Omogoˇca nam, da z reˇsitvijo enega trikotnega sistema pridemo do ustreznih koeficientov interpolacijskega polinoma za dane vrednosti na izbranih interpolacijskih toˇckah. V primeru, da je matrika
V singularna, bomo dobili polinom, ki bo vrnil vrednost 0 na vsaki interpolacijski toˇcki. Z razˇsiritvijo matrike V bomo priˇsli do enoliˇcne reˇsitve. Veˇc o tem v naslednjem poglavju, ko bomo spoznali interpolacijski podprostor minimalne stopnje.
Za toˇcke xj, j = 1,2, . . . ,(︁n+d
d
)︁, kjer je xj ∈ Rd, definiramo kvadratno matriko V = [(xj)λ]. Indeks j teˇce po stolpcih matrike, multiindeks λ pa po vrsticah. Ve- lja|λ| ≤n, pri ˇcemer nam multiindeksλpove stopnjo monoma vsake spremenljivke.
Interpolacijske toˇcke so zdaj poljubne, torej niso veˇc nujno mreˇzne toˇcke, zato uporabimo novo indeksno notacijo. Interpolacijske toˇcke zdaj povezujemo z oznako multiindeksov multi(j), ki se zaˇcne z multi(1) = (0,0, . . . ,0) in nadaljuje glede na linearno urejenost multiindeksov iz formule (7). ˇSe veˇc, zdaj lahko definiramo red interpolacijske toˇcke kot|multi(j)|.
Recimo da imamo podane toˇcke oblike xj = (aj, bj) za n = 2 in d = 2.
Bloki matrike V so urejeni glede na multiindeksno ureditev. Linearni sistem, iz katerega izraˇcunamo koeficiente polinoma (12), je cTV =fT. Pri tem sta
cT = [c00, c10, c01, c20, c11, c02] in fT = [f1, f2, f3, f4, f5, f6]
vrstiˇcna vektorja, pri ˇcemer je fj vrednost funkcije na interpolacijski toˇcki xj, ma- trika V pa je enaka
V =
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
1 1 1 1 1 1
a1 a2 a3 a4 a5 a6 b1 b2 b3 b4 b5 b6 a21 a22 a23 a24 a25 a26 a1b1 a2b2 a3b3 a4b4 a5b5 a6b6
b21 b22 b23 b24 b25 b26
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦ .
Na matriki V bomo izvajali elementarne eliminacije, s katerimi bomo elimini- rali vrednosti pod in nekatere nad diagonalo in s tem priˇsli do ˇzelene zgornje trikotne matrike W. Eliminacije, ki jih izvajamo na matriki V, oznaˇcimo z matrikami Li za i = 1, . . . , k in dobimo LkLk−1. . . L1V = W. Pri tem produkt matrik LkLk−1. . . L1 oznaˇcimo z R, torej imamo opravka s sistemom RV = W. Jasno je, da je matrika R obrnljiva in po mnoˇzenju dobimo V =R−1W. Pri tem je R−1 spodnje trikotna matrika inW zgornje trikotna matrika.
Polinom p iz (12) lahko zdaj zapiˇsemo v novi bazi kot
p(x, y) =[︁
c00 c10 c01 c20 c11 c02]︁
R−1R
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣ 1 x y x2 xy y2
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
=cTR−1R
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣ 1 x y x2 xy y2
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦ .
Ce zdaj vstavimoˇ LU razcep v sistem cTV = fT, dobimo cTR−1W = fT, kjer produkt cTR−1 oznaˇcimo z a. Ko reˇsimo sistem aW = fT, dobimo koeficiente a v novi bazi R[︁
1 x y x2 xy y2]︁
, iz matrike R pa lahko preberemo polinome sploˇsne Newtonove baze.
Pokaˇzimo zgoraj povedano na primeru. Vzemimo kar toˇcke iz slike 1 (c). Interpo- lacijske toˇcke zloˇzimo v matriko V in jo razˇsirimo z dodatkom identiˇcne matrike.
1 x y x2 xy y2
1 1 1 1 1 1 1 1 0 0 0 0 0
x 0 1 2 2 −1 −2 0 1 0 0 0 0
y 0 −1 1 2 2 1 0 0 1 0 0 0
x2 0 1 4 4 1 4 0 0 0 1 0 0
xy 0 −1 2 4 −2 −2 0 0 0 0 1 0
y2 0 1 1 4 4 1 0 0 0 0 0 1
Vsak stolpec na levi strani te sheme se ujema z interpolacijsko toˇcko, katere koordinate so prikazane v vrsticah 2 in 3. Vsak stolpec na desni predstavlja en monom iz standardne baze, seveda enako velja za vrstice na levi polovici sheme. Na tem mestu ˇse enkrat poudarimo, da matriko razdelimo na bloke glede na stopnjo multiindeksov, kar izhaja iz definicije matrike V. Po mnoˇzenju z elementarnimi operacijami, kot je opisano zgoraj, dobimo shemo
1 x y x2 xy y2
1 1 1 1 1 1 1 1 0 0 0 0 0
x 0 1 0 −23 −53 −43 0 13 −23 0 0 0 y 0 0 1 43 13 −13 0 13 13 0 0 0 x2 0 0 0 1 0 0 0 −413 −926 391 134 1978 xy 0 0 0 0 1 0 0 133 135 −439 −313 391 y2 0 0 0 0 0 1 0 −1752 −1152 529 131 521
Koeficienti na levi strani te sheme doloˇcajo matriko W, koeficienti na desni strani pa matriko R. Matrika R definira polinome rλ za vsak vrstiˇcni multiindeks λ. Na primer, r10(x, y) = 13x− 23y. Vrednosti teh polinomov v interpolacijskih toˇckah preberemo iz matrike W po vrsticah. Vsak polinom rλ
(1) ima vrednost 0 na vseh interpolacijskih toˇckah niˇzjega reda,
(2) ima vrednost 1 za predstavljajoˇco interpolacijsko toˇcko in vrednost 0 na vseh ostalih interpolacijskih toˇckah istega reda in
(3) ima stopnjo |λ|.
Polinome oblike rλ imenujemo Newton-Sauerjevi polinomi.
Nadaljujemo s primerom. Recimo, da najdemo interpolacijski polinom za
vrednosti fT = [5,6,7,8,9,10] na ˇsestih vozliˇsˇcih. Iz reˇsitve trikotnega sistema aW =fT, pridemo do koeficientov interpolacijskega polinoma v Newton-Sauerjevi bazi. Z uporabo obratnih substitucij dobimo a = [5,1,2,1,5,7] in iskan interpola- cijski polinom je enak
p(x, y) = 5 +r10(x, y) + 2r01(x, y) +r20(x, y) + 5r11+ 7r02
= 1561 (780−69x+ 15y+ 133x2−48xy+ 79y2).
Spet se vrnimo na vpraˇsanje, kaj ˇce za interpolacijske toˇcke vzamemo toˇcke na kroˇznici, kot je prikazano na sliki 1(a). Izbrane interpolacijske toˇcke so torej
(−1,−√
3),(1,−√
3),(−1,√
3),(1,√
3),(−2,0),(2,0).
V matriko V jih zloˇzimo kar v tem vrstnem redu. Seveda bi jih lahko uredili tudi drugaˇce in bi tako dobili drugaˇcne Newton-Sauerjeve polinome. Z uporabo elemen- tarnih eliminacij na matriki [V|I] dobimo konˇcno matriko
1 x y x2 xy y2
1 1 1 1 1 1 1 1 0 0 0 0 0
x 0 1 0 −12 1 32 12 12 0 0 0 0 y 0 0 1 12 1 12 12 0
√ 3
6 0 0 0
x2 0 0 0 1 0 1 −13 0 0 13 0 0 xy 0 0 0 0 1 1 16 14
√ 3 12
1 12
√ 3 12 0 y2 0 0 0 0 0 0 −4 0 0 1 0 1
Iz zadnje vrstice preberemo polinom p(x, y) = −4 + x2 + y2. Ta polinom vrne vrednost 0 v vseh interpolacijskih toˇckah, kar pomeni, da enoliˇcna reˇsitev ne obstaja. Reˇsitev tega problema bo opisana v nadaljevanju v poglavju 7.
Izrek 6.1. Naj bo {xj} mnoˇzica(︁n+d
n
)︁ interpolacijskih toˇck, ki jih zapiˇsemo v obliki matrike V = [(xj)λ]. Pri tem indeks λ teˇce po vrsticah, indeks j pa po stolpcih.
Matriko V dopolnemo z identiˇcno matriko osnovne monomske baze do [V|I]. Z uporabo elementarnih operacij na vrsticah in menjavo stolpcev matriko preoblikujemo v [W|R], kjer je W zgornje trikotna matrika z identiˇcnimi bloki po diagonali, R pa spodnje trikotna matrika. S pomoˇcjo tega algoritma
(a) dobimo ustrezno ureditev interpolacijskih toˇck xj in Newton-Sauerjeve poli- nome rλ(x) stopnje |λ|, pri ˇcemer so koeficienti v vsaki vrstici λ matrike R podani tako, da je
rλ(xj) =
{︄0; ˇce |multi(j)| ≤ |λ| in multi(j)̸=λ 1; ˇce multi(j) =λ
ali pa
(b) ne uspemo najti reˇsitve, saj so nekatere vrstice ν na levi strani matrike W niˇcelne, kar pomeni, da na desni strani za polinomerν(x)stopnje≤ddobimo vrednost rν(xj) = 0 na vseh interpolacijskih toˇckah.
7. Interpolacijski prostor minimalne stopnje
V tem razdelku bomo poiskali interpolacijski prostor minimalne stopnje, na katerem bo mogoˇce izvesti enoliˇcno interpolacijo z Newtonovim polinomom na poljubno izbranih interpolacijskih toˇckah. Pokazali bomo, kako ravnamo v primeru, ko naletimo na teˇzavo iz izreka 6.1 (b).
Ideja je naslednja, kadarkoli se nam v matriki na levi strani sheme pojavi vrstica, ki vsebuje samo vrednosti 0, to vrstico izbriˇsemo iz matrike in dodamo novo vrstico in nov stolpec, ki predstavljata naslednji monom standardne baze.
Ne smemo pozabiti izbrisati tudi stolpca na desni strani matrike, ki se ujema z izbrisano vrstico, saj ga ne potrebujemo veˇc.
Nadaljujmo s primerom, ko interpolacijske toˇcke izberemo na kroˇznici (slika 1(a)). Naslednji monom po vrsti je x3. Ko dodamo nov stolpec na desni strani matrike, moramo na levi strani dodati novo vrstico za izbrani monom. Nova dodana vrstica v naˇsem primeru bi bila
1 x y x2 xy x3
−1 1 −1 1 −8 8 0 0 0 0 0 1 .
Bloˇcne identiˇcne matrike po diagonali matrike W se skrˇcijo vsakiˇc, ko izbriˇsemo vrstico. V zgornjem primeru, se spremeni samo zadnja vrstica. Prvih pet vrstic ostane nespremenjenih in s pomoˇcjo operacij na vrsticah jo preoblikujemo v
1 x y x2 xy x3 0 0 0 0 0 1 −16 121 0 16 0 121 .
Vidimo, da za vsako izbiro interpolacijskih toˇck v krogu, obstaja enoliˇcna interpola- cija na mnoˇzici Lin{1, x, y, x2, xy, x3}, ki je ekvivalenta mnoˇzici Newton-Sauerjevih polinomov, ki jih najdemo na desni strani matrikeR z novo zadnjo vrstico. Newton- Sauerjeve polinome lahko preberemo iz sheme
1 x y x2 xy x3
1 1 1 1 1 1 1 0 0 0 0 0
0 1 0 −12 1 32 12 12 0 0 0 0 0 0 1 12 1 12 12 0
√3
6 0 0 0
0 0 0 1 0 1 −13 0 0 13 0 0 0 0 0 0 1 1 16 14
√3 12
1 12
√3 12 0 0 0 0 0 0 1 −16 121 0 16 0 121
in dobimo
r00(x, y) = 1, r10(x, y) = 12 +12x, r01(x, y) = 12 +
√ 3 6 y, r20(x, y) =−13 +13x2, r11(x, y) = 16 +14x+
√ 2
12y+121 x2+
√ 3 12xy, r30(x, y) =−16 +121 x+16x2+ 121x3.
Nadgrajeni algoritem deluje za poljubno ˇstevilo interpolacijskih toˇck in poiˇsˇce po- linomski podprostor minimalne stopnje, v katerem je interpolacijski problem vselej korekten. Enakovreden algoritem, ki nam vrne enak rezultat je predstavljen v [5].
8. Zakljuˇcek
Vsi predstavljeni algoritmi nam ˇze z osnovnim znanjem linearne algebre in numeriˇcnih metod odprejo vrata v veˇcrazseˇzno interpolacijo, ki je veliko bolj razgibana. Veˇcrazseˇzna posploˇsitev klasiˇcnega Newtonovega polinoma in algoritma deljenih diferenc je zelo uporabno orodje, kadar imamo ˇskatlasto ali trikotno mnoˇzico interpolacijskih toˇck. Prav tako za poljubno izbrane interpolacijske toˇcke z uporabo algoritma, ki temelji na elementarnih eliminacijah na matriki, dobimo zadovoljiv rezultat interpolacije. Za veˇcrazseˇzne bazne polinome stopnje n imajo lahko Newton-Sauerjevi polinomi vse koeficiente pri monomih iste stopnje neniˇcelne.
Ce algoritem uporabimo na trikotni mnoˇˇ zici interpolacijskih toˇck, nas ta pripelje do klasiˇcnega Newtonovega polinoma, pri ˇcemer Newton-Sauerjev normalizira diagonalne elemente na vrednost 1. Ko interpolacijske toˇcke ne dopuˇsˇcajo reˇsitve za polinome stopnje ≤ n, nam Newton-Sauerjev algoritem pojasni, kako reˇsimo problem, saj dobimo algebraiˇcno hiperpovrˇsino, ki vsebuje vse izbrane interpola- cijske toˇcke in problem reˇsimo z razˇsiritvijo algoritma na iskanje interpolacijskega prostora minimalne stopnje.
Algoritem deljenih diferenc in algoritem na dveh matrikah delujeta povsem razliˇcno. Ko se lotimo reˇsevanja interpolacijskega problema preko algoritma delje- nih diferenc na prvem koraku najdemo ustrezno Newtonovo bazo in na naslednjem koraku iˇsˇcemo ustrezne koeficiente za ta interpolacijski polinom. V primeru, da imamo trikotno mnoˇzico interpolacijskih toˇck vdspremenljivkah omejeno s stopnjo n, je prvi korak zelo preprost, saj bazne polinome le izpiˇsemo. Na drugem koraku pa uporabimo algoritem deljenih diferenc, ki nam vrne ustrezne koeficiente. Za algoritem deljenih diferenc porabimo mn operacij, pri ˇcemer je m = (︁n+d
d
)︁ ˇstevilo razliˇcnih vnosov na n korakih. ˇCe to primerjamo z algoritmom na dveh matrikah, pri ˇcemer je zaˇcetna matrika velikosti m×m, vidimo, da na prvem koraku za LU
razcep porabimo ˇze O(m3) operacij. Na drugem koraku za reˇsevanje trikotnega sistema potrebujemo O(m2) operacij. Zdi se, da je algoritem deljenih diferenc veliko bolj uˇcinkovit.
Algoritmi predstavljeni v diplomski nalogi so baza za nadaljno raziskovanje, razˇsirjanje in izboljˇsave veˇcrazseˇzne interpolacije s pomoˇcjo Newtonovih interpo- lacijskih polinomov na strukturirani mnoˇzici interpolacijskih toˇck ali na podobnih bazah za katerikoli mnoˇzico toˇck.
Slike
1 Primeri izbire interpolacijskih toˇck: (a), (b) in (c). 5
2 Interpolacija na mreˇzi. 10
3 Trikotna shema toˇck. 11
4 Polna shema toˇck. 11
5 Polinom p1. 15
6 Polinom p2. 15
7 Polinom p3. 15
8 Multiindeksi v dokazu izreka formule deljenih diferenc. 19
Tabele
1 Algoritem deljenih diferenc 13
