close all;

%Rugo egyenlet
%x'=y
%y'=-x

rugo_de = @(x)[x(2) -x(1)];

%Kezdeti ertekek specialisan:
%x(0)=0
%y(0)=1

%KEP megoldasa:
%x(t) = sin(t)
%y(t) = cos(t)

%VAGY AHOGY A JEGYZETBEN VAN:
%x(t) = cos(t)
%y(t) = -sin(t)

%Negy Numerikus modszert vizsgalunk ezen a KEP-an.
%Azt tapasztaljuk, hogy minosegileg maskepp viselkednek a parameterek
%fuggvenyeben.
%Kiszamoljuk a rendszerek energiait es abrazoljuk kulonfele koordinatak
%szerint az egyes megoldasokat.

%Jeloles magyarazat:

%StepSize === h
%NumOfSteps === N (ertsd: 1->N)
%time === [ t_0 t_1 ...t_n... t_N-1 ]
%EEuler === Explicit Euler modszer matrixa
%analytic === Analitikus megoldas matrixa
%IEuler === Implicit Euler modszer matrixa
%TEuler === Theta = 0.5 modszer matrixa
%SVEuler === Stormer-Verlet modszer matrixa
%NRG_*E === A rendszer energiajanak vektora az ido szerint
%NRG_AN === A rendszer energiajanak vektora az ido szerint

%Matrix jeloles:
% Matrix/Vektor elso elemenek indexe a Matlabban 1!!!
%*Euler(a,1) === x === x_n
%*Euler(a,2) === y === y_n
%*Euler(a+1,1) === X === x_n+1
%*Euler(a+1,2) === Y === y_n+1

%FELADAT:
%1) Implementaljuk az egyes modszereket a jeloleseknek megfeleloen
%2) Probaljunk ki mas lepeskoz es lepesszam kombinaciokat! Mit latunk?
StepSize=0.2;  % Ez csunya nagy lepeskoz
NumOfSteps=50; % Ez keves lepes 

x_null=0;
y_null=1;
t_null=0;

% Ido kiszamolasa.
%Minden modszer azonos diszkret ido-pontokat hasznal, ezert itt kulon van
%szedve.
time=zeros(NumOfSteps,1);
time(1)=t_null;

for a=1:(NumOfSteps-1)
    time(a+1)=time(a)+StepSize;
end

%Explicit euler method
% (X-x)/h = y
% (Y-y)/h =-x
EEuler=zeros(NumOfSteps,2);

EEuler(1,1)=x_null;
EEuler(1,2)=y_null;

for a=1:(NumOfSteps-1)
    EEuler(a+1,:) = EEuler(a,:) + StepSize .* rugo_de(EEuler(a,:));
end

%Euklideszi norma
NRG_EE=EEuler(:,1).*EEuler(:,1)+EEuler(:,2).*EEuler(:,2);


%Analytic solution
%x(t) = cos(t)
%y(t) = -sin(t)

analytic=zeros(NumOfSteps,2);

analytic(1,1)=x_null;
analytic(1,2)=y_null;

for a=1:(NumOfSteps-1)
    analytic(a+1,1) = sin(time(a+1));
    analytic(a+1,2) = cos(time(a+1));
end

%Gyakorlatilag:
%time=linspace(0,StepSize*NumOfSteps,NumOfSteps) %Valahogy így
%analytic=[sin(time) cos(time)];

%Euklideszi norma
NRG_AN=analytic(:,1).*analytic(:,1)+analytic(:,2).*analytic(:,2);

%Implicit Euler method
% (X-x)/h = Y
% (Y-y)/h =-X
IEuler = zeros(NumOfSteps,2);

IEuler(1,1)=x_null;
IEuler(1,2)=y_null;

for a = 1 : (NumOfSteps-1)
    IEuler(a+1,1)  = (IEuler(a,1)  + (StepSize *IEuler(a,2))) / (1+StepSize*StepSize);
    IEuler(a+1,2)  = (IEuler(a,2)- (StepSize * IEuler(a,1))) / (1+StepSize*StepSize);
end

%Euklideszi norma
NRG_IE=IEuler(:,1).*IEuler(:,1)+IEuler(:,2).*IEuler(:,2);

%Theta
%y(n+1) = y(n)+ h*( theta *f(n+1)+  (1-theta)* f(n))
%Ha theta = 0 akkor explicit euler
%Ha theta = 1 akkor implicit euler
%Minden thetara elsorendu kiveve theta=1/2 amikor masodrendu
%theta=0.5 trapez modszer
%Jelen esetben egy atlagolas:
% (X-x)/h = (y+Y)*0.5
% (Y-y)/h =-(x+X)*0.5

TEuler = zeros(NumOfSteps,2);

TEuler(1,1)=x_null;
TEuler(1,2)=y_null;
theta=0.5;

for a = 1 : (NumOfSteps-1)
    %x(n+1) = ( x(n) + h *(1-theta)* y(n) + h*theta*y(n)-h*h*theta*(theta-1)*x(n) )/(1+h*h*theta*theta)
    %y(n+1) = ( y(n) + h *(theta-1)* x(n) - h*theta*x(n)-h*h*theta*(1-theta)*y(n) )/(1+h*h*theta*theta)
    TEuler(a+1,1)  = ( TEuler(a,1) + StepSize *(1-theta)* TEuler(a,2) + StepSize*theta*TEuler(a,2)+StepSize*StepSize*theta*(theta-1)*TEuler(a,1) )/(1+StepSize*StepSize*theta*theta);
    TEuler(a+1,2)  = ( TEuler(a,2) + StepSize *(theta-1)* TEuler(a,1) - StepSize*theta*TEuler(a,1)-StepSize*StepSize*theta*(1-theta)*TEuler(a,2) )/(1+StepSize*StepSize*theta*theta);
end

%Euklideszi norma
NRG_TE=TEuler(:,1).*TEuler(:,1)+TEuler(:,2).*TEuler(:,2);

%Stormer-Verlet
% (X-x)/h = y
% (Y-y)/h =-X
%=>
%X = x+ h*y
%Y = y -h*x -h*h*y
%ERDEKESSEG: Ebbol levezetheto:
% X*X + Y*Y + h*X*Y = x*x + y*y +h*x*y
% Es invariansan hagyja az X^2 + X^2 +h*x*y = const ellipsziseket

SVEuler = zeros(NumOfSteps,2);

SVEuler(1,1)=x_null;
SVEuler(1,2)=y_null;


for a = 1 : (NumOfSteps-1)
    SVEuler(a+1,1)  = ( SVEuler(a,1) + StepSize * SVEuler(a,2));
    SVEuler(a+1,2)  = ( SVEuler(a,2) - StepSize * SVEuler(a,1) - StepSize*StepSize*SVEuler(a,2));
end

%Euklideszi norma
NRG_SVE=SVEuler(:,1).*SVEuler(:,1)+SVEuler(:,2).*SVEuler(:,2);

%t-szerint
figure
plot(time,EEuler(:,1),'red',time,analytic(:,1),'blue',time,IEuler(:,1),'green',time,TEuler(:,1),'cyan', time, SVEuler(:,1),'black','DisplayName', 't szerint');
axis equal
legend('EE', 'Anal', 'IE', 'theta0.5', 'VE');

%x-y szerint
figure
plot(EEuler(:,1),EEuler(:,2),'red',analytic(:,1),analytic(:,2),'blue',IEuler(:,1),IEuler(:,2),'green',TEuler(:,1),TEuler(:,2),'cyan',SVEuler(:,1),SVEuler(:,2),'black','DisplayName', 'x-y szerint')
axis equal
legend('EE', 'Anal', 'IE', 'theta0.5', 'VE');

%Energiak
figure
plot(time,NRG_EE,'red',time,NRG_AN,'blue',time,NRG_IE,'green',time,NRG_TE,'cyan',time,NRG_SVE,'black','DisplayName', 'Energiak')
axis equal
legend('EE',  'Anal','IE', 'theta0.5', 'VE');

%t-x-y szerint
figure
plot3(time,EEuler(:,1),EEuler(:,2),'red',time,IEuler(:,1),IEuler(:,2),'green',time,analytic(:,1),analytic(:,2),'blue',time,TEuler(:,1),TEuler(:,2),'cyan',time,SVEuler(:,1),SVEuler(:,2),'black','DisplayName', 't-x-y szerint')
grid on
legend('EE', 'IE', 'Anal', 'theta0.5', 'VE');