#!/usr/bin/env python
# coding: utf-8

# <br>
# <center><div style="font-size:45px" id="pocetak"> Interpolacija      </div> </center>   <br>
# 
# 
# Navigacija:<br>
# [Lagrangeov polinom](#lagrange)<br>
# [Newtonov polinom](#njutn)<br>
# [Čebišev polinom](#cebisev)<br>
# [Spline interpolacija](#spline)<br>
# 
# 
# 
# 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.
# <br><br>
# <div style=" border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6;  background-color:#FAFAFA; padding: 5px 5px 5px 15px">
# <b> Teorema </b>(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}$$     
# </div><br><br>
# <center style="font-size:15px"><img src="Slike/interpolacija1.png" width=800px>  
#     <caption> Grafici funkcija $f,f+\epsilon, f-\epsilon$ i polinoma $p$</caption> 
# </center><br><br>
# 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.  
# 

# <b>Zadatak interpolacije.</b> 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*}<br>
# <div style="border-left: 2px solid black; border-right: 1px solid #E6E6E6; border-top: 1px solid #E6E6E6; border-bottom: 1px solid #E6E6E6; background-color:#FAFAFA; padding: 5px 5px 5px 15px" >
# <b> Teorema </b>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).$  
# </div><br>
# 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. <br><br>
# <div style="border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6; border-bottom: 1px solid #E6E6E6; background-color:#FAFAFA; padding: 5px 5px 5px 15px" >
# <b> Teorema </b>
# 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).$
# </div>
# <br>
# <div style="border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6;background-color:#FAFAFA; padding: 5px 5px 5px 15px" >
# <b>Zadatak.</b> 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).$
# </div>
# 
# 
# 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()


# <center style="font-size:15px"><img src="Slike/interpolacija3.png" width=1200px>  
#     <caption></caption> 
# </center><br><br>

# <div style="border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6; background-color:#FAFAFA; padding: 5px 5px 5px 15px" >
# <b> Zadatak.</b> 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).$
# </div>

# 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)


# <center style="font-size:15px"><img src="Slike/interpolacija4a.png" width=1200px>  
#     <center style="font-size:15px"><img src="Slike/interpolacija4b.png" width=1200px>  
#     <caption></caption> 
# </center><br><br>

# <div style="border-left: 2px solid black; border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6;background-color:#FAFAFA; padding: 5px 5px 5px 15px" >
# <b> Zadatak.</b> Interpolirati funkciju $f(x)=e^x\cos(0.75 x)$ na $[-0.5,11.5]$ u 6 tačaka.
# </div>

# 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)


# <center style="font-size:15px"><img src="Slike/interpolacija5.png" width=1200px>  
#    <caption></caption> 
# </center><br><br>

# <div style="border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6; background-color:#FAFAFA; padding: 5px 5px 5px 15px" >
# <b>Napomena.</b> U kodu za prethodni zadatak granice segmenta interpolacije, broj tačaka i samu funkciju jednostavno je zamjeniti promjenom xa, xb, br i f1.
# </div>

# <br>
# <div style="font-size:35px" id="lagrange"> Interpolacioni polinom u Lagrangeovom obliku</div>
# <br>
# 
# Navigacija:<br>
# [Vrati se na početak](#pocetak)<br>
# [Newtonov polinom](#njutn)<br>
# [Čebišev polinom](#cebisev)<br>
# [Spline interpolacija](#spline)
# <br><br>
# 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}$$
# 
# 
# 
# 
# <center style="font-size:20px"><img src="Slike/interpolacija2.png" width=800px>  
#     <caption > Grafici kardinalnih funkcija $l_{i-2},l_{i-1},l_i,l_{i+1},l_{i+2}$</caption> 
# </center>
# 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 <font color="red"> crvenom </font> bojom, obje vrijednosti).   <br><br>
# 
# 
# <div style=" border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6;background-color:#FAFAFA; padding: 5px 5px 5px 15px" >
# <b>Zadatak.</b> 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).$
# </div><br>
# 
# 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)


# <center style="font-size:15px"><img src="Slike/interpolacija6.png" width=1200px>  
#    <caption></caption> 
# </center><br><br>

# <div  style=" border-left: 2px solid black; border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6; background-color:#FAFAFA; padding: 5px 5px 5px 15px" >
# <b>Zadatak </b> Interpolirati funkciju $f(x)=x\sin x$ interpolacionim poinomom u Lagrangeovom obliku na segmentu $[-1,15]$ sa 24 čvorova.  
# </div>

# 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)


# <center style="font-size:15px"><img src="Slike/interpolacija7.png" width=1200px>  
#    <caption></caption> 
# </center><br><br>

# <br><div style="font-size:35px" id="njutn"> Interpolacioni polinom u Newtonovom obliku</div><br>
# 

# 
# 
# Navigacija:<br>
# [Vrati se na početak](#pocetak)<br>
# [Langrangeov polinom](#lagrange)<br>
# [Čebišev polinom](#cebisev)<br>
# [Spline interpolacija](#spline)
# <br><br>
# 
# 

# 
# 
# 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$)<br><br>
# \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*} 
# 
# 
# 
# 
# 
# <div  style="border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;   border-bottom: 1px solid #E6E6E6; background-color:#FAFAFA;  padding: 5px 5px 5px 15px" >
# <b>Zadatak </b> Interpolirati funkciju $f(x)=e^{-0.1x}\sin 5x$ na $[-0.1,2.5]$ koristeći podijeljene razlike u 11 čvorova.  
# </div>

# 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)


# <center style="font-size:15px"><img src="Slike/interpolacija8.png" width=1200px>  
#    <caption></caption> 
# </center><br><br>

# <br>
# <div style="font-size:35px" id="cebisev">Čebiševljevi polinomi</div>
# <br>
# 
# Navigacija:<br>
# [Vrati se na početak](#pocetak)<br>
# [Lagrangeov polinom](#lagrange)<br>
# [Newtonov polinom](#njutn)<br>
# [Spline interpolacija](#spline)<br><br>
# 
# 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](#slikacebisev2)).  
# 
# <center style="font-size:15px" id="slikacebisev1"><img src="Slike/runge1.png" width=1200px>  
#     <caption> Interpolacija Rungeove funkcije korištenjem ekvidistantnih čvorova</caption> 
# </center><br>
# 
# Ovaj nedostatak interpolacije polinomom sa ravnomjerno rasporedjenim čvorovima može se izbjeći korištenjem nula Čebiševljevih polinoma kao čvorova. <br><br>
# 
# <div  style="border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6;background-color:#FAFAFA;padding: 5px 5px 5px 15px" >
# <b>Zadatak.</b> Interpolirati funkciju $f(x)=\cos \pi x$ sa 9 čvorova na segmentu $[-1,1].$<br>
# (a) Koristiti ravnomjerno rasporedjene čvorove.<br>
# (b) Zatim koristiti za čvorove nule Čebiševljevog  polinoma $T_9.$ 
# </div><br>
# 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.
# <center style="font-size:15px"><img src="Slike/cebisev1.png" width=1600px>  
#     <caption> Grafici Čebiševih polinoma</caption> 
# </center><br>
# 
# 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()


# <center style="font-size:15px" id="slikacebisev3">
#       <img src="Slike/cebisev2.png" width=1200px> 
#       <caption> Grafik funkcije i interpolacionog polinoma na ekvidistantnim čvorovima   </caption>
#       <img src="Slike/cebisev2a.png" width=1200px>  
#       <caption> Grafik funkcije i interpolacionog polinoma na čvorovima  koji su nule Čebiševljevog polinoma  </caption>
#       <img src="Slike/cebisev3.png" width=1200px> 
#       <caption> Grafik greške $f-p_8$ sa ekvidistantnim čvorovima   </caption>
#       <img src="Slike/cebisev3a.png" width=1200px> 
#        <caption> Grafik greške $f-p_8$ sa čvorovima  koji su nule Čebiševljevog polinoma   </caption>
#       <img src="Slike/cebisev4.png" width=1200px>  
#       <caption> Grafik polinoma $\omega_9$ sa ekvidistantnim čvorovima  </caption>
#       <img src="Slike/cebisev4a.png" width=1200px>
#       <caption> Grafik polinoma $\omega_9$ sa čvorovima koji su nule Čebiševljevog polinoma </caption>
# </center><br>

# 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](#slikacebisev1)).
# 
# 
# 
# <center style="font-size:15px" id="slikacebisev2">
#       <img src="Slike/runge2.png" width=1200px>  
#       <img src="Slike/runge3.png" width=1200px>  
#     <img src="Slike/runge4.png" width=1200px>  
#     <caption> Interpolacija Rungeove funkcije gdje su čvorovima interpolacije nule Čebiševljevog polinoma</caption> 
# </center><br>

#  <center  style="font-size:35px" id="spline"> Spline interpolacija </center> <br>
#  
#  
# Navigacija:<br>
# [Vrati se na početak](#pocetak)<br>
# [Lagrangeov polinom](#lagrange)<br>
# [Newtonov polinom](#njutn)<br>
# [Čebišev polinom](#cebisev)<br><br>
# 
#  
#  
#  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. 
# 
# 
# 
# <p style="font-size:25px">Linearni spline</p><br>
# 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. 
# 
# 
# <center style="font-size:20px" id="slikacebisev1"><img src="Slike/spline1.png" width=1200px>  
#     <caption> Interpolacija linearnom funkcijom </caption> 
# </center><br>
# 
# 
# 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).$$
# 

# <div  style="border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6;background-color:#FAFAFA;padding: 5px 5px 5px 15px" >
# <b>Zadatak.</b> Interpolirati funkciju $f(x)=e^{-x}\sin 3x$ linearnim splineom sa 14 čvorova na segmentu $[0,5].$<br>
# </div><br>

# 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()


# <center style="font-size:20px" id="slikaspline1"><img src="Slike/linspline1.png" width=1200px>  
#     <caption> Interpolacija linearnim splineom </caption> 
# </center><br>

#  <p style="font-size:25px">Kubni spline</p><br>
# 
# 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.
# 
# 
# <center style="font-size:20px" id="slikacebisev1"><img src="Slike/spline2.png" width=1200px>  
#     <caption> Interpolacija kubnim splineom </caption> 
# </center><br><br>
# 
# 
# 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*}  

# <div  style="border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6;background-color:#FAFAFA;padding: 5px 5px 5px 15px" >
# <b>Zadatak.</b> Interpolirati funkciju $f(x)=e^{-x}\sin 3x$ prirodnim kubnim splineom sa 14 čvorova na segmentu $[0,5].$<br> 
# </div><br>

# 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()


# <center style="font-size:20px" id="slikaspline1"><img src="Slike/kubnispline1.png" width=1200px>  
#     <caption> Interpolacija kubnim splineom </caption> 
# </center><br>

# <div style="font-size:35px"> Literatura </div><br>
# 
# [Bah77]  N.S. Bahvalov, <i>Numerical Methods: Analysis, Algebra, Ordinary Differential Equations</i>, MIR, USSR, 1977. ISBN: 9780714712079.<br><br>
# [Sci15] R. Scitovski, <i>Numerička matematika,</i> Sveučilišta Josipa Jurja Strossmayera u Osijeku – Odjel za matematiku, R.Hrvatska, 2015. ISBN: 978-953-6931-79-8.<b><br>