%3. gyakorlat: rugó egyenlet numerikus megoldásai:

% Paraméterek:

x_null=0;                   % kezdeti feltételek
y_null=1;

% Idõ kiszámolása.
% Minden módszer azonos diszkrét idõpontokat használ
t_null=0; %kezdeti idõpont
t_max=100; %eddig nézem a folyamatot
h=0.1; %lépésköz
time=t_null:h:t_max; %idõpillanatok ,ahol megoldom az egyenletet
N=length(time); %lépésszám

% FELADAT
% Adott a következõ egyenlet: x"+x = 0
% Átírva 2 elsõfokúra:
% x'= y
% y'=-x adódik.

% A DER analitikus megoldása:
% dx(t) = cos(t)
% dy(t) = -sin(t)


% Analitikus megoldás:
analytic=[sin(time)', cos(time)'];

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



% Explicit Euler módszer
% (X-x)/h = y
% (Y-y)/h =-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 .* (-EEuler(a,1));
end

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



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

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

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

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



% Szemi-implicit Euler:
   % (X-x)/h =y
   % (Y-y)/h =-X
SEuler = zeros(N,2);

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

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

% Euklideszi norma négyzet:
NRG_SE=SEuler(:,1).*SEuler(:,1)+SEuler(:,2).*SEuler(:,2);



% ode45:
rugo_ode=@(t,xy) [xy(2); -xy(1)]; %csillapítás nélkül

%b=0.5; %csillapítás %kicsi: rugózik, nagy: csak leül a rendszer
%rugo_ode=@(t, xy) [xy(2); -b*xy(2)-xy(1)]; %csillapított

[timeode45, ode45mg]=ode45(rugo_ode, [t_null, t_max], [x_null, y_null]);
%options=odeset('RelTol', 1e-9, 'AbsTol', 1e-12 ); %az ode45 pontosításására (alap: 1e-3, 1e-6) %inkább az AbsTol számít, ha valamihez kell közel lenni (egyens pont), és a RelTol ha egymáshoz (periódus)?
%[timeode45, ode45mg]=ode45(rugo_ode, [t_null, t_max], [x_null, y_null], options);

% Euklideszi norma négyzet:
NRG_ode45=ode45mg(:,1).*ode45mg(:,1)+ode45mg(:,2).*ode45mg(:,2);



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

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

%x-y szerint
subplot(2,2,2)
plot(analytic(:,1), analytic(:,2), 'black');
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');
xlabel('x');
ylabel('y');
axis equal
legend('Anal', 'EE', 'IE', 'SE', 'ode45');

%Energiak
subplot(2,2,3);
plot(time, NRG_AN, 'black');
hold on;
plot(time, NRG_EE, 'red');
plot(time, NRG_IE, 'green');
plot(time, NRG_SE, 'blue');
plot(timeode45, NRG_ode45, 'magenta');
xlabel('t');
ylabel('Energia');
axis equal
legend('Anal', 'EE', 'IE', 'SE', 'ode45');

%t-x-y szerint
subplot(2,2,4);
plot3(time, analytic(:,1), analytic(:,2), 'black');
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');
grid on
xlabel('t');
ylabel('x');
zlabel('y');
legend('Anal', 'EE', 'IE', 'SE', 'ode45');