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.
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*}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.
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*}
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()
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)
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)
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}$$
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)}.$$
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)
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 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*}
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)
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).
Ovaj nedostatak interpolacije polinomom sa ravnomjerno rasporedjenim čvorovima može se izbjeći korištenjem nula Čebiševljevih polinoma kao čvorova.
(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).$
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.
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.$$
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()
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).
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
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
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).$$
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()
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.
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*}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()
[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.