%% 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-')