%% Trigonometrikus interpoláció
% 1. f(x)=|x-\pi| függvény esetén határozzuk meg az interpolációs trigonometrikus polinomot 4 alappont használatával!

n=3;
x=[0:n]*2*pi/(n+1);
f=abs(x-pi);

plot(x,f,'o')

a0=sum(f)/(n+1)
a1=f*cos(x)'*2/(n+1)
a2=f*cos(2*x)'*2/(n+1)

b1=f*sin(x)'*2/(n+1)

t=@(x)a0+a1*cos(x)+b1*sin(x)+a2*cos(2*x);

hold on
fplot(t,[0,2*pi])
hold off

%
% 2. A [0,2\pi]intervalllumon 6 ekvidisztáns pontban egy függvény értékei rendre az alábbiak: 0,1,1,0,-1,-1. Határozzuk meg a trigonometrikus interpolációs polinomot papíron kézzel számolva a valós együtthatókat megadó képletek segítségével, a komplex együtthatókat megadó mátrix-vektor szorzás segítségével a gyors Fourier-transzformációt alkalmazva, ill. a Matlab fft parancsával!

clear all
n=5;
f=[0,1,1,0,-1,-1];
x=[0:n]*2*pi/(n+1);
plot(x,f,'o')

a0=sum(f)/(n+1)
a1=f*cos(x)'*2/(n+1)
a2=f*cos(2*x)'*2/(n+1)
a3=f*cos(3*x)'/(n+1)

b1=f*sin(x)'*2/(n+1)
b2=f*sin(2*x)'*2/(n+1)

t=@(x)a0+a1*cos(x)+b1*sin(x)+a2*cos(2*x)+b2*sin(2*x)+a3*cos(3*x);

hold on
fplot(t,[0,2*pi])
hold off

%%%%%%%%%%%%%%%%%%%

w=exp(-1i*2*pi/(n+1));
[xi,yi]=meshgrid(0:n,0:n);
F=w.^(xi.*yi)
c=F*f'/(n+1)

F_kicsi=F([1,2,3],[1,3,5])
D=diag(w.^[0:2])
f1=f([1,3,5])'
f2=f([2,4,6])'

F_kicsi_f1=F_kicsi*f1;
DF_kicsi_f2=D*(F_kicsi*f2);

c1=(F_kicsi_f1+DF_kicsi_f2)/(n+1)
c2=(F_kicsi_f1-DF_kicsi_f2)/(n+1)

%%%%%%%%%%%%%%%%%%%

c=fft(f')/(n+1)

%
% 3. A következő kód egy 2\pi periodikus zajos jelből ad meg 256 ekvidisztáns pontban egy mintavételi vektort. Diszkrét Fourier-transzformációt használva, az előző feladat kódjának átírásával simítsuk ki a jelet!

x=linspace(0,2*pi,257);
x=x(1:256);
f=sin(3*x).*cos(4*x)-sin(x).^2;
f=f+rand(1,256)*0.1;
f=f+rand(1,256)*0.3;
plot(x,f)
xlabel('idő')
ylabel('jelerősség')

f_hat=fft(f);
plot(1:30,f_hat(1:30))

%
%% Numerikus differenciálás
% 4. Igazoljuk, hogy az alábbi képlet negyedrendben közelíti egy kellően sima függvény x_0-beli első deriváltját! Hogyan alkalmazható ebben az esetben a Richardson-extrapoláció?
% \frac{-f_2+8f_1-8f_{-1}+f_{-2}}{12h}

format long
syms x;
f=exp(sin(x^2));
x0=1;
der=diff(f);
der_pontos=eval(subs(der,x0))
h=0.1;
kozelites=eval((-subs(f,x0+2*h)+8*subs(f,x0+h)-8*subs(f,x0-h)+subs(f,x0-2*h))/(12*h))


format short

%
% 5. Illesszünk az (x_{-1},f_{-1}),(x_0,f_0),(x_1,f_1)  pontokra interpolációs polinomot és adjuk meg annak első és második deriváltját x_0-ban! Mit kaptunk?

syms  x fa fb fc a b c h
a=b-h
c=b+h

f=fa+(fb-fa)/(b-a)*(x-a)+((fc-fb)/(c-b) - (fb-fa)/(b-a))/(c-a)*(x-a)*(x-b);
fd=diff(f)
expand(subs(fd,b))

fdd=diff(f,2)
expand(subs(fdd,b))


%
%% Numerikus integrálás
% 6. Adjuk meg az N^{4,k}_{\text{zárt}} Newton-Cotes-együtthatókat, és adjunk meg közelítő integrálértéket velük segítségével az alábbi integrálra:
% \int_{0}^\pi \sin x\, dx

syms x
l1=(x-1/4)*(x-1/2)*(x-3/4)*(x-1)/((0-1/4)*(0-1/2)*(0-3/4)*(0-1));
int(l1,0,1)

%
% 7. Vizsgáljuk meg a fenti integrál összetett Simpson-formulás közelítése során a konvergenciarendet! Mekkora legyen a távolság az osztópontok között a {10}^{-6}-os hiba eléréséhez?

clear all
format long
n=2;
h=pi/n;
x=0:h/2:pi;
f=sin(x)

Ipontos=2;
I=h/6*(f(1)+f(end)+4*sum(f(2:2:end-1))+2*sum(f(3:2:end-2)))
hiba=I-Ipontos

format short

% Matlab numerikus integráló parancsa:

quad(@sin,0,pi)

%
% 8. Írjuk át az előző programot úgy, hogy az összetett érintőformulát hajtsa végre. Ellenőrizzük le a konvergenciarendet!
%
% 9. Írjuk fel a négypontos Gauss-Csebisev-formulát a [-1,1] intervallumon! Mekkora a formula pontossági rendje? Ellenőrizzük le!

format long
x=cos(pi/8+[0:3]*pi/4);

syms t
p1=(4*t^7+4*t^6-3*t^5+6*t^4-t^3+4*t^2-t+5);
p=p1/sqrt(1-t^2);
I_pontos=int(p,-1,1)
eval(I_pontos)

q=matlabFunction(p1);
I=sum(q(x))*pi/4

%
% 10. Írjuk fel a hárompontos összetett Gauss-Legendre-formulát az e^{-x^2} függvény [0,1] intervallumon vett integráljának közelítésére (L_4=(5x^3-3x)/2)! Mekkora a formula konvergenciarendje? Ellenőrizzük le!

clear all
syms x
solve(5*x^3-3*x)

f=exp(-x^2);
I_pontos=eval(int(f,0,1))
f=matlabFunction(f);

int((x-0)*(x-sqrt(3/5))/(-sqrt(3/5)-0)/(-sqrt(3/5)-sqrt(3/5)),-1,1)

n=16
h=1/n;
I=h/2*(sum(f(h/2:h:1))*8/9+sum(f(h/2+sqrt(3/5)*h/2:h:1))*5/9+sum(f(h/2-sqrt(3/5)*h/2:h:1))*5/9)
hiba=I_pontos-I

%
%