Edit on GitLab

Interpolacija z zlepki

Interpolacija s polinomi

Poglejmo si najprej primer interpolacije s polinomom. Interpolirali bomo funkcijo $e^{2x}$ v tokah $-1, 0, 1, 2$ z newtonovim interpolacijskim polinomom. Nato napako ocenimo z oceno

\[\begin{equation}\label{eq:napaka} f(x) - p_n(x) = \frac{f^{n+1}(\xi)}{(n+1)!}\Pi_{k=1}^n(x-x_i) \end{equation}\]

in oceno primerjamo z numerično napako.

Newtonova interpolacija

Naj bodo $x_1, \ldots x_{n+1}$ vrednosti neodvisne spremenljivke in $y_1,\ldots y_{n+1}$ vrednosti neznane funkcije. Interpolacijski polinom je polinom, ki zadošča enačbam

\[\begin{equation}\label{eq:enacbe} p(x_1) = y_1, p(x_2)=y_2, \ldots p(x_{n+1})=y_{n+1} \end{equation}\]

Newtonov interpolacijski je interpolacijski polinom

\[p(x) = a_0 + a_1(x-x_1) +a_2(x-x_1)(x-x_2) + \ldots a_n(x-x_1)\ldots (x-x_{n-1}),\]

ki je zapisan v bazi

\[1, x-x_1, (x-x_1)(x-x_2), \ldots (x-x_1)\ldots (x-x_{n-1}).\]

Koeficiente $a_i$ je tako lažje izračunati, kot če bi bil polinom zapisan v standardni bazi. Poleg tega je računanje vrednosti polinoma v standardni bazi lahko numerično nestabilno, če so vrednosti $x_i$ relativno skupaj.

Koeficiente $a_i$ lahko poiščemo bodisi tako, da rešimo spodnje trikotni sistem, ki ga dobimo iz enačb \eqref{eq:enacbe} ali pa z deljenimi diferencami

using Plots
x = [-1, 0, 1, 2]
y = exp.(2x)
slika = scatter(x, y)
-1 0 1 2 0 10 20 30 40 50 y1

Ustvarimo tabelo deljenih diferenc in zapišemo Newtonov interpolacijski polinom

using NumMat
p = deljene_diference(x,y)
plot!(slika, p, -1, 2, label="interpolacijski polinom")
plot!(slika, x->exp(2x), -1, 2, label="exp(2x)")
-1 0 1 2 0 10 20 30 40 50 y1 interpolacijski polinom exp(2x)

Razlika med funkcijo in interpolacijskim polinomom lahko ocenimo s formulo \eqref{eq:napaka}, če upoštevamo, da je $f^{(4)}=16e^{2x}\le 16e^4\simeq 873$. Potem je

\[|p(x)-f(x)|\le 36.4 (x+1)x(x-1)(x-2)\le 36.4\]

medtem ko dejansko napako najlepše prikažemo grafično

plot(x->p(x)-exp(2x), -1, 2, label="napaka", title="Napaka polinomske interpolacije")
-1 0 1 2 -1 0 1 2 3 4 Napaka polinomske interpolacije napaka

Ocena, ki smo jo dobili s formulo \eqref{eq:napaka} je tako precej pesimistična, čeprav je vsaj red velikosti približno pravilen.

Tipi in večlična dodelitev omogoča izražanje podobno kot v matematičnemu jeziku

Na primeru Newtonovega polinoma lahko ilustriramo, kako izkoristimo tipe in večlično dodelitev, da lahko transparentno definiramo operacije, ki so v matematiki naravne, v programskih jeziki pa pogosto njihova uporaba ni več tako intuitivna. Za Newtonov interpolacijski polinom smo tako definirali tip NewtonovPolinom in metodo polyval, s katero lahko izračunamo vrednost Newtonovega interpolacijskega polinoma v dani točki. Lahko gremo še dlje in definiramo, da se vrednosti tipa NewtonovPolinom obnašajo kot funkcije

(p::NewtonovPolinom)(x) = polyval(p, x)

tako da lahko vrednosti polinoma dobimo tako, da kličemo polinom sam, kot bi bil funkcija

p = NewtonovPolinom([1,0,0],[1, 2, 3])
p(1.23)

Rungejev pojav

Pri nizkih stopnjah polinomov se napaka interpolacije zmanjšuje, če povečamo število interpolacijskih točk. A le do neke mere, če stopnja polinoma preveč povečamo, začne napaka spet naraščati. Interpolirajmo funkcijo $\sin(3x)+\cos(2x)$ v točkah $x_i=\frac{i-1}{2}$, za $i=1\ldots 13$.

using Plots
using NumMat
x = 0:0.5:6
f(x) = sin(3x)+cos(2x)
scatter(x, f.(x), title="Interpolacija s polinomom stopnje 12", grid=false)
plot!(f, 0, 6; label="f(x)")
p = deljene_diference(Array(x), f.(x))
plot!(p, 0, 6; label="polinom stopnje 12")
0 1 2 3 4 5 6 -2 -1 0 1 Interpolacija s polinomom stopnje 12 y1 f(x) polinom stopnje 12

Večja napaka se v tem primeru pojavi blizu roba interpolacijskega intervala 

plot(x->p(x)-f(x), 0, 6; title="Napaka interpolacije", label="p12(x)-f(x)")
0 1 2 3 4 5 6 -0.4 -0.2 0.0 0.2 Napaka interpolacije p12(x)-f(x)
Interpolacija s polinomi visokih stopenj je lahko problematična

V prejšnjem primeru opazimo, da se napaka na robu znatno poveča. Povečanje je posledica velikih oscilacij, ki se pojavijo na robu, če interpoliramo s polinomom visoke stopnje na ekvidistančnih točkah. To je znan pojav pod imenom Rungejev pojav in se pojavi, če interpoliramo s polinomom visoke stopnje v ekvidistančnih točkah. Problemu se lahko izognemo, če namesto ekvidistančnih točk uporabimo Čebiševe točke.

Zlepki

Zakaj zlepki

Pogosto je bolje uporabiti različne definicije funkcije na različnih intervalih, kot eno funkcijo z veliko parametri. Funkcijam, ki so sestavljene iz več različnih definicij pravimo zlepki. Če interpoliramo z zlepkom, se izognemo Rungejevemu pojavu in ne dobimo velikih oscilacij na robu.

Linearen zlepek

Interpoliraj funkcijo $sin$ z linearnim zlepkom. Numerično oceni napako.

Hermitovi zlepki

Za točke $x_1, \ldots x_n$ imamo podane vrednosti funkcije $f_1, \ldots f_n$ in vrednosti odvoda $df_1, \ldots df_n$. Poiskati želimo gladek zlepek sestavljen iz polinomov $S_i(x)$ na intervalu $[x_i, x_{i+1}]$, ki interpolira dane podatke. To pomeni, da so izpolnjene naslednje enačbe

\[\begin{equation}\label{eq:hermite} \begin{array}{ll} S_i(x_i) = f_i & S_{i}(x_{i+1})=f_{i+1}\cr S_i'(x_i) = df_i& S_{i}'(x_{i+1})=df_{i+1} \end{array} \end{equation}\]

Hermitov zlepek je vedno zvezen in zvezno odvedljiv, ker se vrednosti funkcije in vrednosti odvoda predpisane v točkah, kjer se dva predpisa dotikata in se zato avtomatično ujemajo.

Newtonova interpolacija vrednosti in odvodov

Če poleg funkcijskih vrednosti podamo še odvode, lahko še vedno uporabimo deljene diference za izračun Newtonovega interpolacijskega polinoma. Pri tem vrednosti neodvisne spremenljivke $x_i$, v katerih je podan tudi odvod, v tabeli deljenih diferenc napišemo večkrat (točko ponovimo tolikokrat, kolikor je podanih višjih odvodov). Poglejmo si primer. Interpolirajmo podatke $x_1=0, x_2=1, f(x_1)=0, f(x_2)=0$ in $f'(x_1)=1, f''(x_1)=-4$. V tabelo deljenih diferenc zapišemo $0,0,0,1$ za vrednosti $x$ in $0, 1, --4, 0$ za vrednosti $f$.

\[\begin{array}{c|cccc} x& f & f[\cdot,\cdot]& f[\cdot,\cdot, \cdot] & f[\cdot,\cdot, \cdot, \cdot]\cr \hline 0 & 0& 1& \frac{-4}{2!}& 1\cr 0 & 1& 1& -1& \cr 0 & -4& 0& & \cr 1 & 0& & & \cr \end{array}\]

Hermitove zlepke sedaj lahko preprosto poiščemo tako da na vsakem intervalu uporabimo deljene diference za vrednosti in odvode v krajščih (enačbe \eqref{eq:hermite}).

zlepek = hermitov_zlepek([0, 1, 2, 4], [0, 1, 2, 0], [0, -1, 0, 1])
plot(zlepek, title="Heremitov zlepek")
0 1 2 3 4 0.0 0.5 1.0 1.5 2.0 Heremitov zlepek

Laplaceovi $\mathcal{C}^1$ zlepki

Za točke $x_1, \ldots x_n$ so podane vrednosti funkcije $f_1, \ldots f_n$. Predpostavimo, da je $n$ liho število. Poiskati želimo zlepek sestavljen iz kubičnih polinomov $S_i(x)$ na intervalu $[x_{2i-1}, x_{2i+1}]$, ki interpolira dane podatke. Poleg tega zahtevamo, da se v točkah, v katerih se predpisi stikajo ujemajo odvodi, da dobimo zvezno odvedljiv zlepek.

Število parametrov naj bo enako številu zahtev

Pri interpolaciji s polinomi lahko zadostimo toliko zahtevanim pogojem, kolikor ima polinom koeficientov oz kolikor je dimenzija prostora poplinomov. Za polinome 3. stopnje je dimenzija 4, kar pomeni, da lahko in moramo predpisati 4 zahteve, bodisi za vrednosti polinoma ali vrednosti odvodov v 4 različnih točkah. Na vsakem intevalu $[x_{2i-1}, x_{2i+1}]$ so 3 vrednosti predpisane (v točkah $x_{2i-1}$, $x_{2i}$, $x_{2i+1}$), to pomeni, da lahko dodamo še eno zahtevo, s katero dosežemo zvezno odvedljivost zlepka.

Omenjene zahteve lahko zapišemo z naslednjimi enačbami

\[\begin{array}{lll} S_i(x_{2i-1}) = f_{2i-1} & S_i(x_{2i})=f_{2i} & S_i(x_{2i+1})&=f_{2i+1}\cr S_{i-1}'(x_{2i-1}) = S_i'(x_{2i-1}) & S_{i}'(x_{2i+1}) = S_{i+1}'(x_{2i+1}) & \end{array}\]

To je sicer 5 zahtev za vsak polinom $S_i$, a če preštejemo vse enačbe in koeficiente, vidimo, da manjka le še ena enačba za polinom na robu $S_1$ ali $S_{(n-1)/2}$. Predpišemo lahko vrednost prvega ali drugega odvoda v eni od robnih točk $x_1$ ali $x_{n}$[1].

Koeficiente polinomov $S_i$ lahko tako računamo zaporedoma tako, da izpolnijo naslednje 4 enačbe

\[\begin{array}{ll} S_i(x_{2i-1}) = f_{2i-1} & S_i(x_{2i})=f_{2i}\cr S_i(x_{2i+1}) =f_{2i+1} & S_i'(x_{2i-1})=S_{i-1}'(x_{2i-1})\cr \end{array}\]

Za prvi polinom $S_1$ pa zadnjo enačbo nadomestimo z enačbo

\[S_1''(x_1)=0\]
Običajno se uporabljajo naravni zlepki

Lagrangeovega $C^1$ zlepka, kot smo ga opisali v tem razdelku običajno ne uporabljamo, saj je lahko za iste podatke sestavimo naravni kubični zlepek, ki ima boljlše lastnosti. Naravni zlepek je dvakrat zvezno odvedljiv in med vsemi funkcijami, ki interpolirajo dane točke ima najmanjšo povprečno ukrivljenost.

Kaj smo se naučili
  • deljene diference in Newtonov polinom lahko uporabimo tudi, če so poleg vrednosti podani tudi odvodi
  • zaradi Rungejevega pojava interpolacija s polinomi visokih stopenj na ekvidistančnih točkah ni najboljša izbira
  • zlepki so enostavni za uporabo, učinkoviti (malo operacij za izračun) in imajo v določenih primerih boljše lastnosti kot polinomi visokih stopenj
  • poznamo različne vrste zlepkov glede na zveznost in podatke, ki jih potrebujemo za določitev (linearni zlepek, Hermitov zlepek in Lagrangeov zlepek)

Koda

Julia omogoča, da definirmo metode za transparentno risanje in računanje z zlepki PlotRecipies.jl

NumMat.deljene_diference — Method
p::NewtonovPolinom{T} = deljene_diference(x::Array{T}, f::Array{U})

Izračuna koeficiente Newtonovega interpolacijskega polinoma, ki interpolira podatke y(x[k])=f[k]. Če se katere vrednosti x večkrat ponovijo, potem metoda predvideva, da so v f poleg vrednosti, podani tudi odvodi.

Primer

Polinom $x^2-1$ interpolira podatke x=[0,1,2] in y=[-1, 0, 3] lahko v Newtonovi obliki zapišemo kot $1 + x + x(x-1)$

julia> p = deljene_diference([0, 1, 2], [-1, 0, 3])
NewtonovPolinom{Float64,Int64}([-1.0, 1.0, 1.0], [0, 1])

Če imamo več istih vrednosti abscise x, moramo v f podati vrednosti funkcije in odvodov. Na primer polinom $p(x) = x^4 = x + 3x(x-1) + 3x(x-1)^2 + x(x-1)^3$ ima v $x=1$ vrednosti $p(1)=1, p'(1)=4, p''(1)=12$

julia> p = deljene_diference([0,1,1,1,2], [0,1,4,12,16])
NewtonovPolinom{Float64,Int64}([0.0, 1.0, 3.0, 3.0, 1.0], [0, 1, 1, 1])

julia> x = (1,2,3,4,5); p.(x).-x.^4
(0.0, 0.0, 0.0, 0.0, 0.0)
source
NumMat.hermitov_zlepek — Method
z::Zlepek = hermitov_zlepek(x, f, df)

Izračuna koeficiente kubičnega C^1 zlepka, ki interpolira podatke $p_i(x_i) = f_i, p_i(x_{i+1}) = f_{i+1}$ in $p_i'(x_i) = df_i, p_i'(x_{i+1})=df_{i+1}$.

Primer

julia> z = hermitov_zlepek([1, 2, 3], [0, 1, 0], [0, 0, 0]);

julia> z.(1:3) - [0, 1, 0]
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

 julia> polyder(z.funkcije[1], 1)
 0
source
NumMat.lagrangev_zlepek — Method
lagrangev_zlepek(x,f)

Izračuna koeficiente kubičnega $C^1$ zlepka, ki interpolira podatke $p_i(x_{2i-1})=f_{2i-1}, p_i(x_{2i})=f_{2i}$ in $p_i(x_{2i+1})=f_{2i+1}$, poleg tega pa je zvezno odvedljiv. Dodatna lasntost zlepka je $p_1''(x_1) = 0$.

Primer

julia> z = lagrangev_zlepek([1, 2, 3, 4, 5], [0, 1, 2, 1, 0]);

julia> z.(1:5) - [0, 1, 2, 1, 0]
5-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
source
NumMat.polyder — Method
dy = polyder(p::NewtonovPolinom, x)

Izračuna vrednost odvoda newtonovega polinoma p v točki x s hornerjevim algoritmom, ki je dopolnjen tako, da računa tudi odvod.

Primer

Vzemimo primer polinoma $1+x^2=1-x+x(x+1)$, katerega odvod je $2x$

julia> p = NewtonovPolinom([1,-1,1], [0,-1]);

julia> odvodp(x) = polyder(p, x); odvodp.((1, 2, 3))
(2, 4, 6)
source
NumMat.polyval — Method
polyval(p::NewtonovPolinom, x)

Izračuna vrednot newtonovega polinoma p v x s Hornerjevo methodo. Objekte tipe NewtonovPolinom lahko kličemo kot funkcijo, ki pokliče to funkcijo.

Primer

julia> p = NewtonovPolinom([0, 0, 1], [0, 1]);

julia> polyval(p, 2)
2

julia> p(2)
2
source
NumMat.NewtonovPolinom — Type
NewtonovPolinom{T}(a::Vector{T}, x::Vector{T})

Vrne newtonov interpolacijski polinom oblike a[1]+a[2](x-x[1])+a[3](x-x[1])(x-x[2])+... s koeficienti a in vozlišči x.

Primer

Poglejmo si polinom $1+x+x(x-1)$, ki je definiran s koeficienti [1,1,1] in z vozlišči $x_0=0$ in $x_1=1$

julia> p = NewtonovPolinom([1, 1, 1], [0, 1])
NewtonovPolinom{Int64}([1, 1, 1], [0, 1])

julia> p.([1,2,3])
3-element Array{Int64,1}:
  2
  5
 10
source
NumMat.Zlepek — Type
Zlepek(funkcije::Vector{F}, vozlisca::Vector{T})

Vrednost tipa Zlepek predstavlja funkcijo podano kot zlepek različnih predpisov na različnih intervalih

\[f(x) = \begin{cases} f_1(x); x\in [x_1, x_2]\cr f_2(x); x\in (x_2, x_3]\cr \vdots \end{cases}\]

Krajišča intervalov so podana v seznamu vozlisca, medtem ko so predpisi zbrani v seznamu funkcije.

Primer

julia> z = Zlepek([x->x, x->2-x], [0, 1, 2]);

julia> z(0.5), z(1.5)
(0.5, 0.5)
source
  • 1Najpogosteje zahtevamo, da je drugi odvod v robnih točkah enak 0. Če predpišemo prvi odvod, določimo tangente v robnih točkah, kar je pogosto preveč vpliva na končno obliko zlepka.