Interpolacija

Navigacija:
Lagrangeov polinom
Newtonov polinom
Čebišev polinom
Spline interpolacija

Korisno je, a nekada i neophodno zamijeniti funkciju $f$ sa nekom funkcijom $p$ koja je "bliska" u nekom smislu sa funkciji $f,$ a čiji analitički oblik možemo lako izračunati, za ovakvu funkciju $p$ kažemo da aproksimira funkciju $f$ i pišemo $f\approx p.$ Kako će se definisati bliskost ovih funkcija zavisi od metrike uvedene u prostoru kojem pripadaju funkcije, te stoga imamo različite tipove zadataka teorije aproksimacije. Funkcija $p,$ u načelu, zavisi od parametara $(a_0,a_1,\ldots,a_n)$ i optimalna bliskost postiže se izborom ovih parametara. Ako je $p$ linearna funkcija parametara $a_i,\,i=0,1,\ldots,n$, aproksimacija je linearna u suprotnom je nelinearna. Pri linearnoj aproksimaciji $p$ se traži u obliku generalisanog polinoma $$ p(x)=a_0\phi_0(x)+\cdots+a_n\phi_n(x), $$ gdje su $\phi_0,\ldots,\phi_n$ linearno nezavisne funkcije koje čine tzv. osnovni sistem funkcija. Funkcije $\phi_i,i=0,\ldots,n$ najčešće biramo iz neke specijalne klase funkcija (stepene, racionalne, trigonometrijske, eksponencijalne funkcije,$\ldots$). Na primjer, ako osnovni sistem čine cijeli nenegativni stepeni argumenta $x,$ $\phi_0(x)=1,\phi_1(x)=x,\ldots,\phi_n(x)=x^n,$ funkcija $p$ je algebarski polinom stepena $n$ $$ p_n(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0, $$ indeks $n$ u oznaci $p_n$ označava stepen polinoma. Ako sada izaberemo funkcije $\phi_i,\,i=0,\ldots,2n-1$ na sljedeći način: $\phi_0(x)=1,\phi_1(x)=\cos x,\phi_2(x)=\sin x,\ldots,\phi_{2n-2}(x)=\cos nx,\phi_{2n-1}(x)=\sin nx,$ $p$ je trigonometrijski polinom stepena $n$ $$ p_n(x)=a_0+a_1\cos x+b_1\sin x+\ldots+a_n\cos x+b_n\sin x. $$

Nastavak će biti posvećen je aproksimaciji funkcije polinomom, a jedan od razloga primjene polinoma u apoksimaciji funkcija je uniformna aproksimacija neprekidnih funkcija. Za datu funkciju neprekidnu na segmentu $[a,b]$, postoji polinom koji je "blizu" te funkcije koliko to mi želimo. Ovaj rezultat je precizno formulisan u sljedećoj Weierstrassovoj teoremi.

Teorema (Weierstrass) Pretpostavimo da je funkcija $f$ definisana i neprekidna na segmentu $[a,b].$ Za svako $\varepsilon>0,$ postoji polinom $p,$ sa osobinom $$\left| f(x)-p(x) \right|<\varepsilon,\:\forall\,x\in[a,b].\label{inter4}$$


Grafici funkcija $f,f+\epsilon, f-\epsilon$ i polinoma $p$


Na prethodnoj slici predstavljeni su grafici fukcije $f$ i polinoma $p.$ Vidimo da se polinom nalazi izmedju grafika funkcija $f-\epsilon$ i $f+\epsilon.$ Po Weierstrassovoj teoremi za proizvoljnu vrijednost $\epsilon$ postoji polinom $p,$ takav da je izmedju $f-\epsilon$ i $f+\epsilon.$ Drugim riječima za datu funkciju $f$ uvijek postoji polinom "blizak" ovoj funkciji, pa funkciju $f$ možemo aproksimirati sa takvim polinomom.

Zadatak interpolacije. Podijelimo segment $[a,b]$ tačkama $x_0,x_1,\ldots,x_n$ na sljedeći način $$ a=x_0<x_1<\ldots<x_n=b. $$ Tačke $x_i,\,i=0,1,\ldots,n$ su čvorovi interpolacije. Neka su i poznate vrijednosti funkcije $f$ u tačkama $x_0,x_1,\ldots, x_n,\:x_i\neq x_j,$ tj. znamo $f(x_0),f(x_1),\ldots,f(x_n).$ Zadatak interpolacije je odrediti, u našem slučaju algebarski polinom $p_n$ stepena $n,$ $$ p_n(x)=a_nx^n+a_{n-1}x^{n-1}+\ldots+a_1x+a_0, $$ takav da je $$ p(x_i)=f(x_i),\,i=0,1,\ldots,n.\quad (\star\star) $$

Polinom $p_n$ je odredjen ako mu znamo koeficijente $a_i,\,i=0,1, \ldots,n,$ kojih ima $n+1,$ a toliko ima i uslova $p(x_i)=f(x_i).$ Na ovaj način dobijamo sistem $(n+1)\times(n+1)$ \begin{align*} a_nx^n_0+a_{n-1}x^{n-1}_0+\ldots+a_1x_0+a_0=&f(x_0)\\ a_nx^n_1+a_{n-1}x^{n-1}_1+\ldots+a_1x_1+a_0=&f(x_1)\quad (\star)\\ \quad\vdots&\\ a_nx^n_n+a_{n-1}x^{n-1}_n+\ldots+a_1x_n+a_n=&f(x_n) \end{align*}

ili u matričnom obliku

\begin{gather*} \left(\begin{array}{ccccc} x^{n}_0 & x^{n-1}_{0}&\cdots&x_0&1\\ x^{n}_1 & x^{n-1}_{1}&\cdots&x_1&1\\ &&\vdots&&\\ x^{n}_n & x^{n-1}_{n}&\cdots&x_n&1\\ \end{array} \right) \left(\begin{array}{c} a_n\\a_{n-1}\\\vdots\\a_0 \end{array} \right)= \left(\begin{array}{c} f(x_0)\\f(x_{1})\\\vdots\\f(x_n) \end{array} \right) \end{gather*}

Teorema Ako su $x_0,x_1,\ldots, x_n$ različite tačke i $f(x_0),f(x_1),\ldots,f(x_n)$ proizvoljni brojevi, tada postoji jedinstven polinom $p_n$ koji zadovoljava uslove $(\star\star).$

Matrica sistema $(\star)$ je Vandermondeova, pa je determinanta matrice sistema različita od nule (vrijedi $x_i \neq x_j$), tj. sistem $(\star)$ ima jedinstveno rješenje, odnosno traženi polinom $p_n$ je jedinstven.

Osim konkretnog računanja interpolacionog polinoma kojim se aproksimira neka funkcija, bitno je i odrediti grešku te aproksimacije, ova se greška uobičajeno zove greška interpolacije.

Teorema Neka su $x_0,x_1,\ldots,x_n$ čvorovi interpolacije i neka je $x$ proizvoljna tačka iz oblasti definisanosti funkcije $f.$ Ako funkcija $f$ ima izvod $(n+1)$-og reda na segmentu $[a,b]$ koji sadrži sve čvorove interpolacije i tačku $x$, tada je greška interpolacije u tački $x$ data sa $$R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_n(x),$$ gdje su $\xi\in[a,b],$ tačka koja zavisi od izbora tačaka $x,x_0,\ldots,x_n$ i funkcije $f,$ $\omega_n=(x-x_0)(x-x_1)\cdot\ldots\cdot(x-x_n)$ i $R_n(x)=f(x)-p_n(x).$

Zadatak. Izračunati interpolacioni polinom za funkciju $f,$ ako je poznata vrijednost funkcije u čvorovima interpolacije $A_1(1,1),\,A_2(2,3),\,A_3(3,7).$

Rješenje. Rješenje. Pošto su date tri tačke $n=2,$ pa ćemo računati interpolacioni polinom u obliku $p_2(x)=a_2x^2+a_1x+a_0.$

Poslije uvrštavanja vrijednosti čvorova datih u tabeli u polinom $p_2$ dobijamo sistem \begin{align*} a_2+a_1+a_0&=1\\ 4a_2+2a_1+a_0&=3\\ 9a_2+3a_1+a_0&=7. \end{align*}

In [ ]:
from IPython import get_ipython;  # briše varijable, funkcije ...
get_ipython().magic('reset -sf')  #

import numpy as np
from matplotlib import pyplot
from numpy import linalg
from numpy import polynomial
from matplotlib import pyplot
from matplotlib import rc
pyplot.style.use('seaborn-whitegrid')
noviparametri={'figure.figsize':(15,7),'font.size': 13,'axes.grid': True,'lines.markersize':10, 'lines.linewidth':2,'figure.dpi':300}
pyplot.rcParams.update(noviparametri)



Polynomial=np.polynomial.Polynomial
solve=linalg.solve
plot=pyplot.plot

a =([1,1,1,],[1,2,4],[1,3,9]);   # funkcija Polynomial vraca polinom od koeficijenta, ali redoslijed koeficijenta u listi je   
b=([1,3,7]);                     # ao, a1,..., an, zbog toga je izmijenjen i redoslijed u listi a 

koeficijenti=[]
koeficijenti=solve(a,b);
print( 'Koeficijenti interpolacionog polinoma su a0=',koeficijenti[0],' ,a1=',koeficijenti[1],' ,a2=',koeficijenti[2])

### prvi način dobijanja polinoma ######
p1=Polynomial(koeficijenti)

### drugi način dobijanja polinoma  #######
broj=len(koeficijenti)
def p2(x):
    poly=0
    for i in range(broj):
        poly=poly+koeficijenti[i]*x**(i)
    return poly

# Tacke interpolacije
xInterpol=[1,2,3]
yInterpol=[1,3,7]

# x koordinate tacaka potrebnih za crtanje grafika
xCrtanje=np.linspace(.8,3.2,100)

# crtanje oba grafika u jednoj vrsti 
fig, ax1 = pyplot.subplots(1,2, sharey=True)

ax1[0].plot(xCrtanje, p1(xCrtanje), 'b',xInterpol,yInterpol,'ro')
ax1[0].set(xlabel='Prvi način računanja polinoma od koeficijenata', ylabel='')

ax1[1].plot(xCrtanje, p2(xCrtanje), 'g',xInterpol,yInterpol,'ro')
ax1[1].set(xlabel='Drugi način računanja polinoma od koeficijenata', ylabel='')

pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/interpolacija3.png')
pyplot.close(fig)
#pyplot.show()
#pyplot.close()



Zadatak. Interpolirati funkciju $f:[a,b]\rightarrow\textbf{R},$ čije su vrijednosti u čvorovima mreže poznate $A_0(0,-3),\,A_1(1,2),\,A_2(2,-3),\,A_3(2.5,1),\,A_4(3,4),\,A_5(4,-1).$

In [ ]:
import numpy as np
from matplotlib import pyplot
from numpy import linalg
from numpy import polynomial
from matplotlib import rc

noviparametri = {'figure.figsize': (15, 7), 'axes.grid': True,
             'lines.markersize': 10, 'lines.linewidth': 2,
             'font.size': 15, 'mathtext.fontset': 'stix',
             'font.family': 'STIXGeneral','figure.dpi':300}
pyplot.rcParams.update(noviparametri)

Polynomial=np.polynomial.Polynomial
solve=linalg.solve
plot=pyplot.plot

# Lista sa koordinatama interpolacionih tacaka
l1=[[0,-3],[1,2],[2,-3],[2.5,1],[3,4],[4,-1]];

# Duzina 
duz1=len(l1)

# Vektor slobodnih clanova, ujedno vektor y koordinata interpolacionih tacaka
vektorSlob=[l1[i][1] for i in range(duz1)]

# Vektor x koordinata interpolacionih tacaka
xInterpol=[l1[i][0] for i in range(duz1)]

# Matrica sistema
matricaSistema=[[(l1[j][0])**i for i in range(duz1)]  for j in range(duz1)] 

rjesenje=solve(matricaSistema,vektorSlob)

# Prvi način dobijanja koeficijenta polinoma
p1=Polynomial(rjesenje)

# Drugi način dobijanja koeficijenta polinoma
broj=len(rjesenje)
def p2(x):
    poly=0
    for i in range(broj):
        poly=poly+rjesenje[i]*x**(i)
    return poly

# Krajnje tacke segmenta za crtanje
xa=-.1;xb=4.5

# x koordinate tacaka za crtanje grafika
xCrtanje=np.linspace(xa,xb,250)

fig=pyplot.figure()
plot(xCrtanje,p1(xCrtanje),'b',xInterpol,vektorSlob,'ro')
pyplot.legend(('Int. polinom','Int. tacke'),loc='lower center',fontsize=20)
pyplot.xlabel('Prvi način računanja polinoma od koeficijenata', fontsize=20)
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/interpolacija4a.png')
pyplot.close(fig)

fig=pyplot.figure()
plot(xCrtanje,p2(xCrtanje),'g',xInterpol,vektorSlob,'ro')
pyplot.legend(('Int. polinom','Int. tacke'),loc='lower center',fontsize=20)
pyplot.xlabel('Drugi način računanja polinoma od koeficijenata', fontsize=20)
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/interpolacija4b.png')
pyplot.close(fig)



Zadatak. Interpolirati funkciju $f(x)=e^x\cos(0.75 x)$ na $[-0.5,11.5]$ u 6 tačaka.

In [ ]:
import numpy as np
from matplotlib import pyplot
from numpy import linalg
from numpy import polynomial
from matplotlib import rc
noviparametri = {'figure.figsize': (15, 7), 'axes.grid': True,
             'lines.markersize': 10, 'lines.linewidth': 2,
             'font.size': 15, 'mathtext.fontset': 'stix',
             'font.family': 'STIXGeneral','figure.dpi':300}
pyplot.rcParams.update(noviparametri)

Polynomial=np.polynomial.Polynomial
solve=linalg.solve
plot=pyplot.plot

def f1(x):
    return np.exp(.1*x)*np.cos(x*.75)

xa=-.5;xb=11.5; br=6;
iks=np.linspace(xa,xb,br)

b=[f1(iks[i]) for i in range(br)]

a=[[(iks[j])**i for i in range(br)]  for j in range(br)]
    
rjesenje=solve(a,b);

p1=Polynomial(rjesenje)
x=np.linspace(xa-.5,xb+.5,250)
y2=f1(x);
fig=pyplot.figure()
plot(x,p1(x),'b',iks,b,'ro',x,y2,'g--')
pyplot.legend(('Int. polinom','Int. tacke','Funkcija $f(x)$'),loc='lower center', fontsize=20)
pyplot.xlabel('Na slici su prikazani grafici funkcije i interpolacionog polinoma',fontsize=20)
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/interpolacija5.png')
pyplot.close(fig)



Napomena. U kodu za prethodni zadatak granice segmenta interpolacije, broj tačaka i samu funkciju jednostavno je zamjeniti promjenom xa, xb, br i f1.


Interpolacioni polinom u Lagrangeovom obliku

Navigacija:
Vrati se na početak
Newtonov polinom
Čebišev polinom
Spline interpolacija

Interpolacioni polinom $p_n$ u Lagrangeovom obliku računamo po formuli $$p_n(x)=f(x_0)\frac{(x-x_1)\cdots(x-x_n)}{(x_0-x_1)\cdots(x_0-x_n)}+\ldots+ f(x_i)\frac{(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)} {(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n)}+\ldots+ f(x_n)\frac{(x-x_0)\cdots(x-x_{n-1})}{(x_n-x_0)\cdots(x_n-x_{n-1})},$$ gdje su $f(x_i)$ vrijednosti funkcije (ili neke vrijednosti dobijene mjerenjem, eksperimentom,...) koju interpoliramo u tačkama mreže $x_i,\,i=0,1,\ldots,n.$ Polinomi $l_i$ su kardinalne ili Lagrangeove funkcije, polinomi su $n$ tog reda za koje vrijedi $$l_i(x_j)=\delta_{ij}=\begin{cases}1,\:j=i\\\\0,\:j\neq i \end{cases}$$

Grafici kardinalnih funkcija $l_{i-2},l_{i-1},l_i,l_{i+1},l_{i+2}$
Sa slike se vidi da je svaka predstavljena kardinalna funkcija $p_i$ u tački interpolacije $x_i$ jednaka $1,$ dok u ostalim tačkama interpolacije ima vrijednost $0$ (tačke označene crvenom bojom, obje vrijednosti).

Zadatak. Izračunati interpolacioni polinom u Lagrangeovom obliku za funkciju $f:[a,b]\rightarrow\textbf{R},$ ako je poznata vrijednost funkcije u čvorovima interpolacije $A_1(1,1),A_2(2,3),A_3(3,7).$

Rješenje. Interpolacioni u Lagrangeovom obliku računamo po formuli $$p_2(x)=f(x_0)\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0+x_2)}+f(x_1)\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1+x_2)}+f(x_2)\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2+x_1)}.$$

In [ ]:
import numpy as np
from matplotlib import pyplot
from matplotlib import rc

noviparametri = {'figure.figsize': (15, 7), 'axes.grid': True,
             'lines.markersize': 10, 'lines.linewidth': 2,
             'font.size': 15, 'mathtext.fontset': 'stix',
             'font.family': 'STIXGeneral','figure.dpi':300}
pyplot.rcParams.update(noviparametri)
plot=pyplot.plot

# Interpolacione tacke
data=np.array([[1,1],[2,3],[3,7]])
# x koordinate interpolacionih tacaka
xData=data[:,0]
# y koordinate interpolacionih tacaka
yData=data[:,1]

# funkcija za dobijanje interpolacionog polinoma p_2 u Lagrangeovom obliku
def p2(x):
    return yData[0]*(x-xData[1])*(x-xData[2])/((xData[0]-xData[1])*(xData[0]-xData[2])) \
           +yData[1]*(x-xData[0])*(x-xData[2])/((xData[1]-xData[0])*(xData[1]-xData[2])) \
            +yData[2]*(x-xData[0])*(x-xData[1])/((xData[2]-xData[0])*(xData[2]-xData[1])) 

fig=pyplot.figure()
xCrtanje=np.linspace(.9,3.1,100)
plot(xCrtanje,p2(xCrtanje),'b',xData,yData,'ro')
pyplot.legend(('Int. polinom','Int. tačke'),fontsize=20) 
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/interpolacija6.png')
pyplot.close(fig)



Zadatak Interpolirati funkciju $f(x)=x\sin x$ interpolacionim poinomom u Lagrangeovom obliku na segmentu $[-1,15]$ sa 24 čvorova.

In [ ]:
import math 
import numpy as np
from matplotlib import pyplot 
from numpy import linalg
from matplotlib import rc
noviparametri = {'figure.figsize': (15, 7), 'axes.grid': True,
             'lines.markersize': 10, 'lines.linewidth': 1.75,
             'font.size': 15, 'mathtext.fontset': 'stix',
             'font.family': 'STIXGeneral','figure.dpi':300}
pyplot.rcParams.update(noviparametri)

# Funkcija za generisanje f(x)=x sin x i odgovarajucih nizova za x i y koordinate
def funkcija(x):
    return x*np.sin(x)
xInterpol=np.linspace(-1,15,24)
yInterpol=funkcija(xInterpol)
   
# Funkcija za generisanje polinoma iz podataka
def pol(x):
    poly=0;
    for i in range(len(xInterpol)):
        def baza(x):
            brojnik=1;
            nazivnik=1;
            for j in range(len(xInterpol)):
                if j!=i:
                    brojnik=brojnik*(x-xInterpol[j])
                    nazivnik=nazivnik*(xInterpol[i]-xInterpol[j])
            return brojnik/nazivnik    
        poly=poly+yInterpol[i]*baza(x)        
    return poly
    
# Nizovi sa x i y kooordinatama vrijednostima trazenog polinoma
xCrtanje=np.linspace(min(xInterpol)-.05,max(xInterpol)+.05,100)
yCrtanje=pol(xCrtanje)

fig=pyplot.figure()
pyplot.plot(xCrtanje,yCrtanje,'b')
pyplot.plot(xInterpol,yInterpol,'ro')
pyplot.plot(xCrtanje,funkcija(xCrtanje),'r-.')
pyplot.legend(('Interpolacioni polinom','Interpolacione tačke','Funkcija $f(x)=x\,\sin\, x$'))
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/interpolacija7.png')
pyplot.close(fig)




Interpolacioni polinom u Newtonovom obliku

Interpolacioni polinom u Newtonovom obliku dobićemo primjenom algoritma podijeljene razlike. Interpolacioni polinom računamo u obliku $$p_n(x)=f[x_0]+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+\ldots +f[x_0,\ldots,x_n](x-x_0)(x-x_1)\cdot\ldots\cdot(x-x_{n-1}),$$ gdje je $$f[x_0]=f(x_0), \quad f[x_i,\ldots,x_j]=\frac{f[x_{i+1},\ldots,x_j]-f[x_i-x_{j-1}]}{x_j-x_i}.$$

Predstavljena je tabela podijeljenih razlika za 4 čvora ($x_0,x_1,x_2.x_3$)

\begin{align*} x &\quad f[\:\:\:] & f[\:\:\:,\:\:\:]\: &\quad f[\:\:,\:\:,\:\:] & f[\:\:,\:\:,\:\:,\:\:]\quad\quad \\\hline x_0 &\quad f[x_0] & & & \\ %\hline &\quad & f[x_0,x_1] & & \\ %\hline x_1 &\quad f[x_1] & &\quad f[x_0,x_1,x_2] & \\%\hline & & f[x_1,x_2] & & f[x_0,x_1,x_2,x_3] \\%\hline x_2 &\quad f[x_2] & &\quad f[x_1,x_2,x_3] & \\%\hline & & f[x_2,x_3] & & \\%\hline x_3 &\quad f[x_3] & & & %\vspace{.5cm}%\hline \end{align*}

Zadatak Interpolirati funkciju $f(x)=e^{-0.1x}\sin 5x$ na $[-0.1,2.5]$ koristeći podijeljene razlike u 11 čvorova.

In [ ]:
import numpy as np
from matplotlib import pyplot
#plot=pyplot.plot
#show=pyplot.show()
from numpy import polynomial, linalg
solve=linalg.solve
Polynomial=polynomial.Polynomial
from matplotlib import rc

noviparametri = {'figure.figsize': (15, 7), 'axes.grid': True,
             'lines.markersize': 10, 'lines.linewidth': 2,
             'font.size': 15, 'mathtext.fontset': 'stix',
             'font.family': 'STIXGeneral','figure.dpi':300}
pyplot.rcParams.update(noviparametri)

##### interpolacione tacke ####################################
#xInterpol=[1,2,3,4,5,6,7,7.5,2.5,.5,1]
#yInterpol=[-1,2,3,-3,2,4,-1,0,1,1,-4]

# Broj cvorova i krajnje tacke segmenta interpolacije
brojCvorova=11
xa=-.1
xb=2.5

# Nizovi sa x koordinatama za interpolaciju i crtanje
xInterpol=np.linspace(xa,xb,brojCvorova)
xCrtanje=np.linspace(xa-.05,xb+.05,250)

# Funkcija za f(x)=e^{0.1} sin 5x  
def f1(x):
    return np.sin(5*x)*np.exp(-.1*x)
# Niz sa y koordinatama interpolacionih tacaka
yInterpol=f1(xInterpol)

# Inicijalizacija matrice razlika 
razlike=np.zeros((brojCvorova,brojCvorova))
# Prva kolona sa elementima f(x-i)
for i in range(brojCvorova):
    razlike[i,0]=f1(xInterpol[i])

##### racunanje koeficijenata f[x_0,...,x_i] ################################################    
for j in range(1,brojCvorova):
    for i in range(j,brojCvorova):
        razlike[i,j]=(razlike[i,j-1]-razlike[i-1,j-1])/(xInterpol[i]-xInterpol[i-j])

##########  polinom od koeficijenata #######################################
def pol1(iks,lista,xlista):
    duzina=len(xlista)-1
    suma=lista[duzina,duzina]
    for i in range(duzina-1,-1,-1):
        suma=suma*(iks-xlista[i])+lista[i,i]
    return suma

fig=pyplot.figure()
pyplot.plot(xInterpol,yInterpol,'ro')
pyplot.plot(xCrtanje,pol1(xCrtanje,razlike,xInterpol),'g')
pyplot.plot(xCrtanje,f1(xCrtanje),'b--')
pyplot.legend(('interpolacione tačke','interpolacioni polinom','funkcija $f(x)=e^{0.1x}\sin 5x$'))
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/interpolacija8.png')
pyplot.close(fig)




Čebiševljevi polinomi

Navigacija:
Vrati se na početak
Lagrangeov polinom
Newtonov polinom
Spline interpolacija

Na sljedećoj slici predstavljeni su grafici Rungeove fukcije $f(x)=\frac{1}{1+x^2}$ i interpolacionih polinoma $p_8$ i $p_{10}.$ Korišteni su ekvidistantni čvorovi, sa grafika se vidi da je greška veća pri interpolaciji polinomom većeg stepena (pogledati grafik interpolacije Rungeove funkcije uz korištenje nula Čebiševljevih polinoma kao čvorova interpolacije, kliknuti na link).

Interpolacija Rungeove funkcije korištenjem ekvidistantnih čvorova

Ovaj nedostatak interpolacije polinomom sa ravnomjerno rasporedjenim čvorovima može se izbjeći korištenjem nula Čebiševljevih polinoma kao čvorova.

Zadatak. Interpolirati funkciju $f(x)=\cos \pi x$ sa 9 čvorova na segmentu $[-1,1].$
(a) Koristiti ravnomjerno rasporedjene čvorove.
(b) Zatim koristiti za čvorove nule Čebiševljevog polinoma $T_9.$

Rješenje. Interpolacioni polinom računamo u obliku $p_8(x)=(x-x_0)(x-x_1)\cdots(x-x_8)$ u oba slučaja.

(a) U prvom slučaju čvorove biramo $x_i=-1+\frac{2}{8}i,\,i=0,\ldots,8.$

Greška interpolacije je $R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x),$ gdje je $\xi\in(a,b)$ i $\omega_{n+1}(x)=(x-x_0)(x-x_1)\cdots(x-x_n).$

In [ ]:
import numpy as np
from matplotlib import pyplot
from numpy import linalg
from numpy import polynomial
from matplotlib import rc

noviparametri = {'figure.figsize': (15, 7), 'axes.grid': True,
             'lines.markersize': 10, 'lines.linewidth': 2,
             'font.size': 15, 'mathtext.fontset': 'stix',
             'font.family': 'STIXGeneral','figure.dpi':300}
pyplot.rcParams.update(noviparametri)

Polynomial=np.polynomial.Polynomial
solve=linalg.solve

brojCvorova=9
xData=np.linspace(-1,1,brojCvorova)
yData=np.zeros(brojCvorova)

##########  interpolacija ####################
def f1(x):
    return np.cos(x*np.pi)

#iks=np.linspace(-1,1,brojCvorova)


matricaSistema=np.zeros((brojCvorova,brojCvorova))
vektorSlob=np.zeros(brojCvorova)

for i in range(brojCvorova):
    vektorSlob[i]=f1(xData[i])
    for j in range(brojCvorova):
        matricaSistema[j,i]=(xData[j])**i
# Drugi način formiranja matrice sistema a i vektora slobodnih clanoova b
#b=[f1(iks[i]) for i in range(n1+1)]
#a=[[(iks[j])**i for i in range(n1+1)]  for j in range(n1+1)]        

rjesenje=solve(matricaSistema,vektorSlob);

p1=Polynomial(rjesenje)
xCrtanje=np.linspace(-1.05,1.05,250)
yCrtanje=f1(xCrtanje);

fig=pyplot.figure()
pyplot.plot(xCrtanje,p1(xCrtanje),'b',xData,vektorSlob,'ro',linewidth=1)
pyplot.plot(xCrtanje,yCrtanje,'g--',linewidth=2)
pyplot.xlim(-1.05,1.05)
pyplot.ylim(-1.05,1.05)
pyplot.legend(('Interpolacioni polinom','Cvorovi interpolacije','Funkcija $f(x)$'),loc='lower center')
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/cebisev2.png')
pyplot.close(fig)
#pyplot.show()

##### greska ######
def f(x3,n,nule):
    if n==1:
        return (x3-nule[0])*(x3-nule[1])
    else:
        return (x3-nule[n])*f(x3,n-1,nule)

nule1=np.zeros(brojCvorova)

fig=pyplot.figure()
pyplot.plot(xCrtanje,(f1(xCrtanje)-p1(xCrtanje)),'g',xData, nule1,'ro',linewidth=2) 
pyplot.legend(('Grafik greske','Cvorovi interpolacije'),loc="upper center")
pyplot.xlim(-1.015,1.015)
pyplot.ylim(-.0002,.0004)
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/cebisev3.png')
pyplot.close(fig)
#pyplot.show()

fig=pyplot.figure()
xCrtanje1=np.linspace(-1.015,1.015,250)
pyplot.xlim(-1.01,1.01)
pyplot.ylim(-.02,.02)
pyplot.plot(xCrtanje1,f(xCrtanje1,brojCvorova-1,xData),'b',xData,yData,'ro')
pyplot.legend(('Polinom $\omega_9$','Cvorovi interpolacije'),loc='upper center')
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/cebisev4.png')
pyplot.close(fig)
#pyplot.show()

(b) Čebiševljev polinom $T_n$ definišemo sa $$T_n(x)=\cos(n \arccos x),\quad n=0,1,\ldots$$ pri čemu je $x\in[-1,1].$ Iz identiteta $$\cos(n+1)\theta+\cos(n-1)\theta=2\cos x\cos n\theta,$$ za $\theta(x)=\arccos x$ dobijamo $$T_{n+1}(x)=2xT_n(x)-T_{n-1}(x),\:n=1,2,\ldots$$ Za $n=0$ i $n=1$ uzimamo da je $T_0(x)=1,$ i $T_1(x)=x,$ respektivno, pa je sada \begin{align*} T_2(x)&=2x^2-1\\ T_3(x)&=4x^3-3x\\ T_4(x)&=8x^4-8x^2+1\\ T_5(x)&=16x^5-20x^3+5x\\ &\vdots \end{align*} Grafici nekoliko Čebiševljevih polinoma ($T_2,\,T_3,\,T_4,\,T_5$) dati su na sljedećoj slici.

Grafici Čebiševih polinoma

Interpolirajmo ponovo funkciju $f(x)=\cos x,$ ali neka su sada čvorovi interpolacije nule Čebiševljevog polinoma, tj. $$x_i=\cos \frac{(2k+1)\pi}{2(n+1)},\,k=0,\ldots,n,\:n=8.$$

In [ ]:
import numpy as np
from matplotlib import pyplot
from numpy import linalg
from numpy import polynomial
from matplotlib import rc

newparams = {'figure.figsize': (15, 7), 'axes.grid': True,
             'lines.markersize': 10, 'lines.linewidth': 2,
             'font.size': 15, 'mathtext.fontset': 'stix',
             'font.family': 'STIXGeneral','figure.dpi':300}
pyplot.rcParams.update(newparams)

Polynomial=np.polynomial.Polynomial
solve=linalg.solve
plot=pyplot.plot

brojCvorova=9
xData=np.zeros(brojCvorova)
for i in range(brojCvorova):
    xData[i]=np.cos((2*i+1)*np.pi/(2*brojCvorova) )

def proizvod(x,n,nule):                    
    if n==1:
        return (x-nule[0])*(x-nule[1]) 
    else:
        return (x-nule[n])*proizvod(x,n-1,nule)

###############################################
##########  interpolacija #####################
###############################################

def f1(x):
    return np.cos(x*np.pi)

### Formiranja matrice sistema a i vektora slobodnih clanova b ####
matricaSistema=np.zeros((brojCvorova,brojCvorova))
vektorSlob=np.zeros(brojCvorova)
for i in range(brojCvorova):
    vektorSlob[i]=f1(xData[i])
    for j in range(brojCvorova):
        matricaSistema[j,i]=(xData[j])**i

##### Interpolacioni polinom  ##########        
rjesenje=solve(matricaSistema,vektorSlob);
p1=Polynomial(rjesenje)

##### Crtanje grafika ##################
xCrtanje=np.linspace(-1.05,1.05,250)
yCrtanje=f1(xCrtanje);

fig=pyplot.figure()
plot(xCrtanje,p1(xCrtanje),'b',xData,vektorSlob,'ro',linewidth=1)
plot(xCrtanje,yCrtanje,'g--',linewidth=2)
pyplot.xlim(-1.05,1.05)
pyplot.ylim(-1.05,1.05)
pyplot.legend(('Interpolacioni polinom','Čvorovi interpolacije','Funkcija $f(x)$'),loc='lower center')
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/cebisev2a.png')
pyplot.close(fig)
#pyplot.show()

##### grafik greske #####
nule1=np.zeros(brojCvorova)

fig=pyplot.figure()
pyplot.plot(xCrtanje,(f1(xCrtanje)-p1(xCrtanje)),'g',xData, nule1,'ro',linewidth=2) 
pyplot.legend(('Grafik greške','Čvorovi interpolacije'),loc="upper center")
pyplot.xlim(-1.001,1.001)
pyplot.ylim(-.0001,.0001)
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/cebisev3a.png')
pyplot.close(fig)
#pyplot.show()
    
###########  polinom \omega_9 ############################
fig=pyplot.figure()
xCrtanje1=np.linspace(-1.,1.,250)
pyplot.plot(xCrtanje1,proizvod(xCrtanje1,brojCvorova-1,xData),'b',xData,nule1,'ro')
pyplot.legend(('Polinom $\omega_n$','Nule polinoma'),loc='lower center')
pyplot.xlim(-1.001,1.001)
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/cebisev4a.png')
pyplot.close(fig)
#pyplot.show()

Grafik funkcije i interpolacionog polinoma na ekvidistantnim čvorovima Grafik funkcije i interpolacionog polinoma na čvorovima koji su nule Čebiševljevog polinoma Grafik greške $f-p_8$ sa ekvidistantnim čvorovima Grafik greške $f-p_8$ sa čvorovima koji su nule Čebiševljevog polinoma Grafik polinoma $\omega_9$ sa ekvidistantnim čvorovima Grafik polinoma $\omega_9$ sa čvorovima koji su nule Čebiševljevog polinoma

Na sljedećoj slici su grafici interpolacije Rungeove funkcije uz korištenje nula Čebiševljevih polinoma kao čvorova interpolacije. Sa grafika se vidi da se poćavanjem broja čvorova smanjuje greška interpolacije, a što nije bio slučaj kod interpolacije sa ekvidistantnim tačkama (kliknuti na link).

Interpolacija Rungeove funkcije gdje su čvorovima interpolacije nule Čebiševljevog polinoma

Spline interpolacija

Navigacija:
Vrati se na početak
Lagrangeov polinom
Newtonov polinom
Čebišev polinom

Kada želimo postići manju grešku interpolacije funkcije polinomom jednostavno povećamo stepen polinoma, medjutim ovakav pristup ne daje uvijek željene rezultate, primjer interpolacija Rungeove funkcije. Drugi način rješavanja problem interpolacije je korištenje neekvidistantnih čvorova interpolacije, primjer nule Čebiševljevih polinoma kao čvorovi interpolacije. Sljedeći pristup je je korištenje polinoma nižeg reda, što nas vodi do spline interpolacije. Segment $[a,b]$ na kome želimo aproksimirati neku funkciju $f$ podijelimo čvorovima $$a=x_0<x_1<\ldots<x_n=b.$$ Na svakom podsegmentu $[x_{i-1},x_i],\,i=1,\ldots,n$ dovoljno malom, funkciju $f$ možemo aproksimirati dovoljno dobro, čak i polinomima stepena $1$ ili $0.$ Traženi spline $s$ treba da zadovoljava uslov $$s(x_i)=f(x_i),\,i=0,\ldots,n.$$ Splineove se u opštem slučaju definišemo na sljedeći način. Za funkciju $\varphi$ defnisanu na segmentu $[a,b]$ kažemo da je spline $k$-tog reda u donosu na čvorove $x_i,\:i=0,\ldots, n$ ako ima sljedeće osobine

  1. $\varphi$ je polinom $k$-tog stepena na svakom od segmenta $[x_i,x_{i+1}],\,i=0,\ldots,n-1;$
  2. $\varphi\in C^{k-1}[a,b],$

gdje je sa $C^{k}[a,b]$ označen skup svih funkcija definisanih na segment $[a,b]$ koje na tom skupu imaju neprekidan izvod $k-1$-tog reda.

Linearni spline


Najjednostavniji i najranije korišteni spline je linearni interpolacijski spline, tj. spline prvog reda odnosno stepena. Grafik ovog splinea je izlomljena linija (slika ispod). Ovaj spline je neprekidna dio po dio linearna funkcija. Pretpostavimo da je $f∶[a, b]\mapsto \mathbb{R}$ neprekidna funkcija i da poznajemo njene vrijednosti $f(x_i),\,i = 0,\ldots, n$ u $n + 1$ jednom čvoru.

Interpolacija linearnom funkcijom

Funkciju $f$ interpoliraćemo neprekidnom po dijelovima linearnom funkcijom $\varphi:[a,b]\mapsto \mathbb{R}.$ Lako se vidi da ovakva funkcija postoji i da je jedinstvena. Za $i=1,\ldots,n-1$ definišimo

$$\varphi\bigg|_{[x_i,x_{i+1}]}=\varphi_i,$$

gdje je $\varphi_i$ linearna funkcija koja prolati kroz tačke $(x_i,f(x_i)$ i $(x_{i+1},f(x_{i+1})$ te vrijedi

$$\varphi_i(x)=f(x_i)+\frac{f(x_{i+1})-f(x_i)}{x_{i+1}+x_i},\:x\in[x_i,x_{i+1}]$$

i $$\varphi(x_0)=\varphi_0(x_0)=f(x_0),\quad \varphi(x_n)=\varphi_{n-1}(x_n)=f(x_n).$$

Zadatak. Interpolirati funkciju $f(x)=e^{-x}\sin 3x$ linearnim splineom sa 14 čvorova na segmentu $[0,5].$

In [ ]:
import numpy as np
from matplotlib import pyplot
from numpy import linalg
from numpy import polynomial
from matplotlib import rc
pyplot.style.use('seaborn-whitegrid')
newparams={'figure.figsize':(22,12),'font.size':(20),
           'lines.markersize':10,'axes.grid':True,
           'lines.linewidth':2,'figure.dpi':300}
pyplot.rcParams.update(newparams)

a=0; b=5; n=14;n1=1000;

iks=np.linspace(a,b,n)
iks2=np.linspace(a,b,n1)

def f1(x):
    return np.exp(-x)*np.sin(3*x)

ipsilon=f1(iks);
    
def f2(x):
    i=int(x*(n-1)/(b-a))
    if i<n-1:
        return (ipsilon[i+1]-ipsilon[i])*(x-iks[i])/(iks[i+1]-iks[i])+ipsilon[i]
    else:
        return f1(b)

spline=np.zeros(n1)
for j in range(n1):
    spline[j]=f2(iks2[j])

y2=f1(iks2)
fig=pyplot.figure()
pyplot.plot(iks,ipsilon,'bo',iks2,y2,'g-.',linewidth=2);
pyplot.plot(iks2,spline,'r-',linewidth=2);
pyplot.legend(('Tačke interpolacije','Funkcija','Linearni spline'))
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/linspline1.png')
pyplot.close(fig)
#pyplot.show()

Interpolacija linearnim splineom

Kubni spline


Linearni interpolacijski spline i kvadratni interpolacijski spline (spline drugog stepena) nisu pogodni za primjenu zato što nisu glatke, odnosno dovoljno glatke funkcije u čvorovima interpolacije. Zato se u praksi koriste i splineovi većeg reda. Posebno je značajan kubni interpolacioni spline (spline trećeg stepena), kratko kubni spline. Ovo su splineovi najmanjeg mogućeg stepena koji imaju neprekidan drugi izvod. Kroz primjene pokazalo se da splineovi većeg stepena rijetko daju bitnija poboljšanja.

Cilj nam je da neprekidnu funkciju $f:[a,b]\mapsto\mathbb{R},$ čije vrijednosti $f(x_i),\:i=0,\ldots n$ znamo u $n+1$ čvoru $$a=x_0<x_1<\ldots<x_n=b,$$ interpolirati funkcijom $C:[a,b]\mapsto\mathbb{R},$ $$ C(x)=C_i(x),\quad x\in\left[x_{i},x_{i+1}\right],\quad i=0,\ldots,n-1, $$

gdje su $C_i,\:i=0,\ldots,n-1,$ kubni polinomi. Na svakom podsegmentu $[x_{i},x_{i+1}],\:i=0,\ldots,n-1$ kubni polinomi su oblika $C_i(x)=c_{i,3}x^3+c_{i,2}x^2+c_{i,1}x+c_{i,0},$ dakle treba na svakom od $n$ podsegmenata odrediti $4$ koeficijenta, ukupno $4n$ koeficijenata na $[a,b],$ tj. treba nam $4n$ uslova postaviti za sistem od $4n$ jednačine. Ovo ćemo uraditi na sljedeći način, svaki kubni polinom $C_i$ u krajnjim tačkama podsegmenta $[x_{i},x_{i+1}],\:i=0,\ldots,n-1$ je jednak vrijednosti funkcije $f,$ tj. zadovoljava interpolacione uslove, zatim prvi i drugi izvodi funkcije $C$ su neprekidni u unutrašnjim čvorovima, tj. u $x_1,\ldots,x_{n-1}.$ Prethodne uslove možemo i ovako zapisati

(U1) $C_i(x_{i})=f(x_{i}),\:i=0,\ldots,n-1;$

(U2) $C_i(x_{i+1})=f(x_{i+1}),\:i=0,\ldots,n-1;$

(U3) $C'_i(x_{i+1})=C'_{i+1}(x_{i+1}),\:i=0,\ldots,n-2;$

(U4) $C''_i(x_{i+1})=C''_{i+1}(x_{i+1}),\:i=0,\ldots,n-2.$

(U1) i (U2) daju $2n$ jednačina, (U3) i (U4) po $n-1,$ sada ukupno imamo $4n-2$ jednačine, dakle nedostaju dvije jednačine. Dodatne dvije jednačine dobijamo ispunjavanjem jednog od sljedeća tri uslova

(U5a) $C''_0(x_0)=C''_{n-1}(x_n)=0;$

(U5b) $f$ i $C$ su periodične funkcije na $[a,b];$

(U5c) $C'_0(a)=f'(a),\quad C'_{n-1}(b)=f'(b).$

U nastavku koristićemo (U5a) i ovakav spline zovemo prirodni kubni spline.

Interpolacija kubnim splineom


Označimo vrijednosti funkcije $f$ u čvorovima interopolacije sa $f_i,$ tj. $$f_i=f(x_i),\:i=0,\ldots,n,$$ zatim dužinu podsegmenta sa $$h_{i+1}=x_{i+1}-x_i,\:i=0,\ldots,n-1,$$ te sa $M_i$ momente traženog splinea $C$ $$M_0:=C''_0(x_0),\:M_n:=C''_{n-1}(x_n),\:M_i:=C''_i(x_i),\:i=1,\ldots,n-1.$$

Na svakom podsegmentu $[x_i,x_{i+1}]$ kubne funkcije $C_i$ računamo po formuli

\begin{multline*} C_i(x)=M_i\frac{(x_{i+1}-x)^3}{6h_{i+1}}+M_{i+1}\frac{(x-x_i)^3}{6h_{i+1}} +\left[\frac{f_{i+1}-f_i}{h_{i+1}}-\frac{h_{i+1}}{6}(M_{i+1}-M_i) \right](x-x_i) +f_i-M_i\frac{h^2_{i+1}}{6} \end{multline*}

gdje su $M_i,\:i=1,\ldots,n-1$ rješenje sistema

\begin{multline*} \frac{h_i}{6}M_{i-1}+\frac{h_i+h_{i+1}}{3}M_i+\frac{h_{i+1}}{6}M_{i+1}= \frac{f_{i+1}-f_i}{h_{i+1}}-\frac{f_i-f_{i-1}}{h_i},\:i=1,\ldots,n-1;\:M_0=M_n=0. \end{multline*}

Zadatak. Interpolirati funkciju $f(x)=e^{-x}\sin 3x$ prirodnim kubnim splineom sa 14 čvorova na segmentu $[0,5].$

In [ ]:
import numpy as np
from matplotlib import pyplot
from numpy import linalg
from numpy import polynomial
from matplotlib import rc
pyplot.style.use('seaborn-whitegrid')
newparams={'figure.figsize':(22,12),'font.size':(20),
           'lines.markersize':10,'axes.grid':True,
           'lines.linewidth':2,'figure.dpi':300}
pyplot.rcParams.update(newparams)

a=0;b=5;n=14; n1=1000;

def f1(x):
    return np.exp(-x)*np.sin(3*x)

iks=np.linspace(a,b,n+1)

y=np.zeros(n+1)
y=f1(iks)   

iks1=np.linspace(a,b,n1)
y1=f1(iks1)

ni=np.zeros(n-1)
mi=np.zeros(n-1)
f=np.zeros(n-1)
matricasistema=np.zeros((n-1,n-1))

ni[0:(n-1)]=(iks[2:n+1]-iks[1:(n)])/(iks[2:n+1]-iks[0:(n-1)])
mi[0:(n-1)]=(iks[1:(n)]-iks[0:(n-1)])/(iks[2:n+1]-iks[0:(n-1)])
f[0:(n-1)]=(6/(iks[2:n+1]-iks[0:(n-1)]))\
      *((y[2:n+1]-y[1:n])/(iks[2:n+1]-iks[1:n])-(y[1:n]-y[0:n-1])/(iks[1:n]-iks[0:n-1]))

matricasistema[0][0]=2; matricasistema[0][1]=ni[1]
matricasistema[n-2][n-3]=mi[n-2];matricasistema[n-2][n-2]=2;

for i in range(1,n-2):
    matricasistema[i][i-1]=mi[i+1]
    matricasistema[i][i]=2
    matricasistema[i][i+1]=ni[i+1]

momenti=linalg.solve(matricasistema,f)  
momenti=np.append(momenti,[0])
momenti=np.append([0],momenti)

def f2(x):
    i=int(x*(n)/(b-a))

    if i<n:
        return y[i]+((y[i+1]-y[i])/(iks[i+1]-iks[i])-(iks[i+1]-iks[i])*(2*momenti[i]+momenti[i+1])/6)*(x-iks[i])\
             +.5*momenti[i]*(x-iks[i])**2+(momenti[i+1]-momenti[i])*(x-iks[i])**3/(6*(iks[i+1]-iks[i]))
    else:
        return f1(b)
    
spline=np.zeros(n1)
for j in range(n1):
    spline[j]=f2(iks1[j])

fig=pyplot.figure()
pyplot.plot(iks,y,'bo',markersize=8)
pyplot.plot(iks1,f1(iks1),'g-.',linewidth=2)
pyplot.plot(iks1,spline,'r')
pyplot.legend(('Interpolacione tačke','Funkcija','Kubni spline'))
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/kubnispline1.png')
pyplot.close(fig)
#pyplot.show()

Interpolacija kubnim splineom

Literatura

[Bah77] N.S. Bahvalov, Numerical Methods: Analysis, Algebra, Ordinary Differential Equations, MIR, USSR, 1977. ISBN: 9780714712079.

[Sci15] R. Scitovski, Numerička matematika, Sveučilišta Josipa Jurja Strossmayera u Osijeku – Odjel za matematiku, R.Hrvatska, 2015. ISBN: 978-953-6931-79-8.