close all;

%Inga egyenlet
%x'=y
%y'=-sin(x)

% Energiafuggveny: 0.5*y*y + 1-cos(x)

inga_de = @(x)[x(2); -sin(x(1))];

%Kezdeti ertekek specialisan:
%x(0)=pi/2
%y(0)=0

%Negy (egyenlore meg csak harom) 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
%SEuler === "Semi-Euler" modszer matrixa
%VEuler === 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?
%3) Nezzuk meg az energia fuggvenyeket is!
StepSize=0.1;
NumOfSteps=100;

x_null=pi/2;
y_null=0;
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
EEuler=zeros(NumOfSteps,2);

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

% X = x +h*y
% Y = y -h*sin(x)

% Ide jon a megoldas.
for i=1:(NumOfSteps-1)
EEuler(i+1,1)=EEuler(i,1)+StepSize*EEuler(i,2);
EEuler(i+1,2)=EEuler(i,2)-StepSize*sin(EEuler(i,1));
end
%Energiafuggveny
% Energiafuggveny: 0.5*y*y + 1-cos(x)
NRG_EE=0.5.*EEuler(:,2).*EEuler(:,2)+1.-cos(EEuler(:,1));
% Ide jon a megoldas.

% % %Implicit Euler method

% Ezt majd jovo oran.

%"Semi-euler"
SEuler = zeros(NumOfSteps,2);

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


% X = x + h*y
% Y = y - h*sin(X)

% Ide jon a megoldas.
for i=1:(NumOfSteps-1)
SEuler(i+1,1)=SEuler(i,1)+StepSize*SEuler(i,2);
SEuler(i+1,2)=SEuler(i,2)-StepSize*sin(SEuler(i+1,1));
end
%Energiafuggveny
NRG_SE=0.5.*SEuler(:,2).*SEuler(:,2)+1.-cos(SEuler(:,1));
% Ide jon a megoldas.

VEuler = zeros(NumOfSteps,2);

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

% X = x + h * y - (h*h/2)*sin(x)
% Y = y - h/2 * sin(x) - h/2*sin(X)

% Ide jon a megoldas.
for i=1:(NumOfSteps-1)
VEuler(i+1,1)=VEuler(i,1)+StepSize.*VEuler(i,2)-StepSize.*StepSize.*0.5.*sin(VEuler(i,1));
VEuler(i+1,2)=VEuler(i,2)-StepSize.*0.5.*sin(VEuler(i,1))-StepSize.*0.5.*sin(VEuler(i+1,1));
end
%Energiafuggveny
NRG_VE=0.5.*VEuler(:,2).*VEuler(:,2)+1.-cos(VEuler(:,1));
% Ide jon a megoldas.

% % (1) EZT StepSize=0.1 es NumOfSteps=100 al
% % (2) ES StepSize=0.001;  % Kis lepeskoz NumOfSteps=100000; % Ez sok lepes, de az ossz ido csak 100
% 
% %t-szerint
 figure
 plot(time,EEuler(:,1),'red',time,SEuler(:,1),'blue',time,VEuler(:,1),'magenta','DisplayName', 't szerint');
 axis equal
 xlabel('time')
 ylabel('x')
 legend('Explicit Euler', 'Semi Euler', 'Verlet');
% 
% %x-y szerint
 figure
 plot(EEuler(:,1),EEuler(:,2),'red',SEuler(:,1),SEuler(:,2),'blue',VEuler(:,1),VEuler(:,2),'magenta','DisplayName', 'x-y szerint')
 axis equal
 xlabel('x')
 ylabel('y')
 legend('Explicit Euler', 'Semi Euler', 'Verlet');
% 
% %Energiak
 figure
 plot(time,NRG_EE,'red',time,NRG_SE,'blue',time,NRG_VE,'magenta','DisplayName', 'Energiak')
 axis equal
 xlabel('time')
 ylabel('energia')
 legend('Explicit Euler', 'Semi Euler', 'Verlet');
% 
% %t-x-y szerint
 figure
 plot3(time,EEuler(:,1),EEuler(:,2),'red',time,SEuler(:,1),SEuler(:,2),'blue',time,VEuler(:,1),VEuler(:,2),'magenta','DisplayName', 't-x-y szerint')
 legend('Explicit Euler', 'Semi Euler', 'Verlet');
 xlabel('time')
 ylabel('x')
 zlabel('y');
 grid on

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

% % (3) EZT StepSize=0.1 es NumOfSteps=1000 el
% 
% %t-szerint
% figure
% plot(time,SEuler(:,1),'blue',time,VEuler(:,1),'magenta','DisplayName', 't szerint');
% axis equal
% 
% %x-y szerint
% figure
% plot(SEuler(:,1),SEuler(:,2),'blue',VEuler(:,1),VEuler(:,2),'magenta','DisplayName', 'x-y szerint')
% axis equal
% 
% %Energiak
% figure
% plot(time,NRG_SE,'blue',time,NRG_VE,'magenta','DisplayName', 'Energiak')
% 
% %t-x-y szerint
% figure
% plot3(time,SEuler(:,1),SEuler(:,2),'blue',time,VEuler(:,1),VEuler(:,2),'magenta','DisplayName', 't-x-y szerint')
% grid on
