%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         3.gyakorlat         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nemlin. din. rsz-ek alapjai %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         2016.09.28.         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     Numerikus módszerek     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Paraméterek:
h=0.1;                      % lepeskoz
N=500;                      % lepesszam

x_null=pi/2;                   % kezdeti feltételek
y_null=0;
t_null=0;

% Idõ kiszámolása.
% Minden módszer azonos diszkrét idõpontokat használ, ezért itt külön van
% szedve.
time=zeros(N,1);            % N*1-es mátrix végig 0 elemekkel
time(1)=t_null;             % kezdeti érték megadása
for a=1:(N-1)               % "idõ mátrix" feltöltése
    time(a+1)=time(a)+h;
end
% FELADAT
% Adott a következõ egyenlet: x"+sin(x) = 0
% Átírva 2 elsõfokúra:
% x'= y
% y'=-sin(x) adódik.

% Energiafüggvénye: 0.5*y^2 + 1-cos(x)

%ode45
inga_ode=@(t,xy) [xy(2); -sin(xy(1))];
[timeode45, ode45mg]=ode45(inga_ode, [t_null, h*N], [x_null, y_null]);

% Energiafüggvény
NRG_ode45=0.5*ode45mg(:,2).^2.+1.-cos(ode45mg(:,1));



% Explicit Euler módszer
% X = x + h*y
% Y = y -h*sin(x)

EEuler=zeros(N,2);

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

for a=1:(N-1)
    EEuler(a+1,1) = EEuler(a,1) + h .* EEuler(a,2);
    EEuler(a+1,2) = EEuler(a,2) - h .* sin(EEuler(a,1));
end

% Energiafüggvény
NRG_EE=EEuler(:,2).*EEuler(:,2)./2.+1.-cos(EEuler(:,1));



% Implicit Euler

% Szemi-implicit Euler
SEuler = zeros(N,2);

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


% X = x + h*y
% Y = y - h*sin(X)
for a = 1 : (N-1)
    SEuler(a+1,1)  = SEuler(a,1) + h * SEuler(a,2);
    SEuler(a+1,2)  = SEuler(a,2) - h * sin(SEuler(a,1) + h*SEuler(a,2));
end

% Energiafüggvény
NRG_SE=SEuler(:,2).*SEuler(:,2)./2.+1.-cos(SEuler(:,1));





IEuler = zeros(N,2);

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

eps = 1.0e-6; % Fixpont tétel toleranciája

% Ennek a fixpontját keressük:
%X_i+1 = x + h*(y - h*sin(X_i))=x + h* Y ;

%X_f fixpontra megoldjuk az egyenletet
%X = X_f;               % Ettõl implicit
%Y = y - h*sin(X);      % Ettõl implicit

for a = 1 : (N-1)
    
    % Fixpont tetel:
    X(1) = IEuler(a,1);
    i=1;
    done=false;
    
    while ~done 
        X(i+1) = IEuler(a,1)+h*(IEuler(a,2) - h*sin(X(i)));
        i=i+1;
        if abs(X(i)-X(i-1)) < eps
            done= true;
        end
    end
    i
    IEuler(a+1,1) = X(i);                               % Ettõl implicit
    IEuler(a+1,2) = IEuler(a,2) - h*sin(IEuler(a+1,1)); % Ettõl implicit
end

% Energiafüggvény
NRG_IE=IEuler(:,2).*IEuler(:,2)./2.+1.-cos(IEuler(:,1));



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        Kirajzoltatás        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



figure;
%t-szerint
subplot(2,2,1);
hold on;
plot(time,EEuler(:,1),'red');
plot(time,IEuler(:,1),'green');
plot(time,SEuler(:,1),'blue');
plot(timeode45, ode45mg(:,1),'magenta');
hold off;
% egybe minden:
%plot(time,EEuler(:,1),'red',time,IEuler(:,1),'green',time,SEuler(:,1),'blue', timeode45, ode45mg(:,1),'magenta');
xlabel('t');
ylabel('x');
axis equal
legend('EE', 'IE', 'SE', 'ode45');

%x-y szerint
subplot(2,2,2)
hold on;
plot(EEuler(:,1),EEuler(:,2),'red');
plot(IEuler(:,1),IEuler(:,2),'green');
plot(SEuler(:,1),SEuler(:,2),'blue');
plot(ode45mg(:,1), ode45mg(:,2), 'magenta');
hold off;
% egybe minden:
%plot(EEuler(:,1),EEuler(:,2),'red',IEuler(:,1),IEuler(:,2),'green',SEuler(:,1),SEuler(:,2),'blue', ode45mg(:,1), ode45mg(:,2), 'magenta')
xlabel('x');
ylabel('y');
axis equal
legend('EE', 'IE', 'SE', 'ode45');

%Energiak
subplot(2,2,3);
hold on;
plot(time,NRG_EE,'red');
plot(time,NRG_IE,'green');
plot(time,NRG_SE,'blue');
plot(timeode45, NRG_ode45,'magenta');
hold off;
% egybe minden:
%plot(time,NRG_EE,'red',time,NRG_IE,'green',time,NRG_SE,'blue', timeode45, NRG_ode45,'magenta')
xlabel('t');
ylabel('Energy');
axis equal
legend('EE', 'IE', 'SE', 'ode45');

%t-x-y szerint
subplot(2,2,4);
hold on;
plot3(time,EEuler(:,1),EEuler(:,2),'red');
plot3(time,IEuler(:,1),IEuler(:,2),'green');
plot3(time,SEuler(:,1),SEuler(:,2),'blue');
plot3(timeode45, ode45mg(:,1), ode45mg(:,2),'magenta');
hold off;
% egybe minden:
%plot3(time,EEuler(:,1),EEuler(:,2),'red',time,IEuler(:,1),IEuler(:,2),'green',time,SEuler(:,1),SEuler(:,2),'blue', timeode45, ode45mg(:,1), ode45mg(:,2),'magenta')
grid on
xlabel('t');
ylabel('x');
zlabel('y');
legend('EE', 'IE', 'SE', 'ode45');
