%Rugo egyenlet
%x'=y
%y'=-x
%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)

daspect([1,1,1]);
close all;
StepSize=0.2;  % Ez csunya nagy lepeskoz
NumOfSteps=50;
x_null=0;
y_null=1;
t_null=0;

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

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

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

   %Ide jon a megoldas
   for i=1:NumOfSteps-1
       EEuler(i+1,1)=StepSize*EEuler(i,2)+EEuler(i,1);
       EEuler(i+1,2)=-StepSize*EEuler(i,1)+EEuler(i,2);
   %time(i+1)=time(i)+StepSize;
   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);
time=zeros(NumOfSteps,1);

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

   %Ide jon a megoldas
   for i=1:NumOfSteps
       analytic(i,1)=cos(-pi/2+time(i));
       analytic(i,2)=-sin(-pi/2+time(i));
   time(i+1)=time(i)+StepSize;
   end

%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);
time=zeros(NumOfSteps,1);

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

   %Ide jon a megoldas
   for i=1:NumOfSteps-1
       IEuler(i+1, 1)=(IEuler(i,1)+StepSize*IEuler(i,2))/(1+(StepSize^2));
       IEuler(i+1, 2)=(IEuler(i,2)-StepSize*IEuler(i,1))/(1+(StepSize^2));
  % time(i+1)=time(i)+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);
time=zeros(NumOfSteps,1);

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

   %Ide jon a megoldas
   %Segitseg:
       %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)
       %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)
    for i=1:NumOfSteps-1
        %????????
       TEuler(i+1,1)=(TEuler(i,1)+StepSize*(1-theta)*TEuler(i,2)+StepSize*theta*TEuler(i,2)+(StepSize^2)*theta*(theta-1)*TEuler(i,1))/(1+(StepSize^2)*(theta^2));
       TEuler(i+1,2)=(TEuler(i,2)+StepSize*(theta-1)*TEuler(i,1)-StepSize*theta*TEuler(i,1)-(StepSize^2)*theta*(1-theta)*TEuler(i,2))/(1+(StepSize^2)*(theta^2));
  % time(i+1)=time(i)+StepSize;
   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
%Ebbol levezetheto:
% X*X + Y*Y + h*X*Y = x*x + y*y +h*x*y
% Invariansan hagyja az X^2 + X^2 +h*x*y = const ellipsziseket

SVEuler = zeros(NumOfSteps,2);
time=zeros(NumOfSteps,1);
SVEuler(1,1)=x_null;
SVEuler(1,2)=y_null;
time(1)=0;

   %Ide jon a megoldas
   for i=1:NumOfSteps-1
       SVEuler(i+1,1)=SVEuler(i,1)+StepSize*SVEuler(i,2);
       SVEuler(i+1,2)=SVEuler(i,2)-StepSize*SVEuler(i,1)-StepSize^2*SVEuler(i,2);
   time(i+1)=time(i)+StepSize;
   end
   
%Euklideszi norma
NRG_SVE=SVEuler(:,1).*SVEuler(:,1)+SVEuler(:,2).*SVEuler(:,2);

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

 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');
 
 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');
 
 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', 'x-y-t szerint')
 grid on
legend('EE', 'IE', 'Anal', 'theta0.5', 'VE');