%% *Numerikus integrálás*
% *1. 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
l0=(x-1/4)*(x-1/2)*(x-3/4)*(x-1)/((0-1/4)*(0-1/2)*(0-3/4)*(0-1));
N40=int(l0,0,1)
N44=N40
l1=(x-0)*(x-1/2)*(x-3/4)*(x-1)/((1/4-0)*(1/4-1/2)*(1/4-3/4)*(1/4-1));
N41=int(l1,0,1)
N43=N41
N42=(1-2*(N40+N41))

I=eval((pi-0)*(N40*sin(0)+N41*sin(pi/4)+N42*sin(pi/2)+N43*sin(3*pi/4)+N44*sin(pi)))

%% 
% *2. 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? Alkalmazzuk a quad parancsot az integrál 
% kiszámítására!*
%%
clear all
format long
n=4;
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)
%% 
% *3. Írjuk át az elõzõ programot úgy, hogy az összetett érintõformulát 
% hajtsa végre! Ellenõrizzük le a konvergenciarendet!*
% 
% *4. Í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=eval(int(p,-1,1))

q=matlabFunction(p1);
I=sum(q(x))*pi/4
%% 
% *5. Írjuk fel a négypontos ö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=(35x^4-30x^2+3)/8$*)! 
% Mekkora a formula konvergenciarendje? Ellenõrizzük le!*
%%
clear all
syms x
zh=eval(solve(35*x^4-30*x^2+3))
w=prod(x-zh);
for i=1:4
suly(i)=int(w/(x-zh(i))/subs(diff(w),zh(i)),-1,1);
end
suly=eval(suly)

f=exp(-x^2);
I_pontos=evalf(int(f,0,1)) % Itt Matlabbal kell kiszámolni a pontos integrált is.
f=matlabFunction(f);

n=2
h=1/n;
I=0;
for i=1:4
I=I+sum(f(h/2*(1+zh(i)):h:1))*suly(i);
end
I=I*h/2;
hiba=I_pontos-I

%% Kezdetiértékfeladatok megoldása
% *6. Oldjuk meg az *$y^\prime=\text{arctg}y, y(0)=1$* kezdetiértékfeladatot 
% az explicit Euler-módszer segítségével a [0,1] intervallumon! Adjunk becslést 
% arra, hogy mekkora legyen a lépéstávolság ahhoz, hogy a hiba már kisebb legyen, 
% mint *$10^{-4}$*!*
%%
f=@(x,y)atan(y);
x0=0, y0=1, xmax=1, N=20

h=(xmax-x0)/N;

x=zeros(1,N+1); y=x;
x(1)=x0;
y(1)=y0;

for i=1:N
    x(i+1)=x(i)+h;
    y(i+1)=y(i)+h*f(x(i),y(i));
end

plot(x,y,'bo-')
%% 
% *7. Oldjuk meg az *$y^\prime=x^2(1+y), y(0)=3$* kezdetiértékfeladatot 
% a RK4-módszer segítségével! Vizsgáljuk meg a módszer konvergenciarendjét! Írjuk 
% át a módszert a Heun-módszerre, és így is vizsgáljuk meg a rendet!*
%%
ypontos=@(x)4*exp(x.^3/3)-1;
f=@(x,y)x.^2.*(1+y);
x0=0, y0=3,
xmax=1, N=80, xx=linspace(x0,xmax,500);

h=(xmax-x0)/N;

x=zeros(1,N+1); y=x;
x(1)=x0;
y(1)=y0;

for i=1:N
    x(i+1)=x(i)+h;
    k1=f(x(i),y(i)); k2=f(x(i)+h/2,y(i)+k1*h/2); k3=f(x(i)+h/2,y(i)+k2*h/2); k4=f(x(i)+h,y(i)+k3*h);
    y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);
end

plot(x,y,'bo-',xx,ypontos(xx),'r-')

hiba=y(end)-ypontos(xmax)

%% 
% *8. Oldjuk meg az elõzõ feladatot a Matlab általános KEF-megoldó parancsával: 
% az ODE45-tel!*
%%
f=@(x,y)x.^2.*(1+y);
x0=0, y0=3, xmax=1

ode45(f,[x0,xmax],y0)
pause

% vagy

[x,y]=ode45(f,[x0,xmax],y0);
plot(x,y)
%% 
% *9. Oldjuk meg az alábbi csillapított rezgõmozgást leíró kezdetiértékfeladatot 
% a [0,15] intervallumon a RK4-módszer segítségével!*
% 
% $$my^{\prime\prime}+cy^\prime+ky=0$$
% 
% *Legyen m=20kg, k=20N/m, c=5, 40, 200 Ns/m, kezdõsebesség 0m/s, kezdõ elmozdulás 
% a nyugalmi helyzetbõl 1m!*
%%
m=20, k=20, c=5
f=@(x,y)[y(2);(-c*y(2)-k*y(1))/m];
x0=0, y0=[1;0], xmax=15, N=80

h=(xmax-x0)/N;

x=zeros(1,N+1); 
y=zeros(2,N+1);
x(1)=x0;
y(:,1)=y0;

for i=1:N
    x(i+1)=x(i)+h;
    k1=f(x(i),y(:,i)); k2=f(x(i)+h/2,y(:,i)+k1*h/2); k3=f(x(i)+h/2,y(:,i)+k2*h/2); k4=f(x(i)+h,y(:,i)+k3*h);
    y(:,i+1)=y(:,i)+h/6*(k1+2*k2+2*k3+k4);
end

plot(x,y,'bo-')
%% 
% *10. Oldjuk meg az alábbi reakciókinetikai kezdetiértékfeladatot a [0,50] 
% idõintervallumon! A kezdeti érékek rendre: 1,1 és 0. Próbáljuk ki a megoldást 
% a Matlab ode45, ode23s függvényeivel ill. implicit Euler-módszerrel!*
% 
% $$c_1^\prime=-0.013c_1-1000c_1c_3\\ c_2^\prime=-2500c_2c_3\\ c_3^\prime=-0.013c_1-1000c_1c_3-2500c_2c_3$$
%%
f=@(x,y)[-0.013*y(1)-1000*y(1)*y(3);-2500*y(2)*y(3);-0.013*y(1)-1000*y(1)*y(3)-2500*y(2)*y(3)];

% ode45(f,[0,50],[1,1,0])

% ode23s(f,[0,50],[1,1,0])

N=50;
h=50/N;

x=zeros(1,N+1); 
x(1)=0;
y=zeros(3,N+1);
y(:,1)=[1;1;0];

J=@(y)[-0.013-1000*y(3)-1/h,0,-1000*y(1);0,-2500*y(3)-1/h,-2500*y(2);-0.013-1000*y(3),-2500*y(3),-1000*y(1)-2500*y(2)-1/h];

for i=1:N
    x(i+1)=x(i)+h;
    z=y(:,i);
    for k=1:10
    z=z-J(z)\(y(:,i)+h*f(x(i+1),z)-z);
    end
    y(:,i+1)=z;
end

plot(x,y,'bo-')