%% TP4 
%% EX 1 
% oscillateur harmonique
%y" = -y, y(0)=0,
% s'ecrit sous forme de système
% Y' = A Y
%avec Y = [y_1, y_2]^T
%A = [0 ,1; -1, 0]
% attention aux point virgules qui séparent les colonnes 
[t,y]=ode45(@(t,y) [y(2);-y(1)], [0,2*pi],[0;1]);
%% tracé des  composantes
% on voit sin et cos
close all;% ferme toutes les figures.
subplot(2,1,1);
plot(t,y(:,1));
xlabel('t');
ylabel('y');
subplot(2,1,2);
plot(t,y(:,2));
xlabel('t');
ylabel('dy/dt');
%
%% tracé de la trajectoire
% on voit un cercle dans l'espace des phases (y,y')
% car y(t)^2 + y'(t)^2 = Cte
figure;
plot(y(:,1),y(:,2),'.-');
xlabel('q');
ylabel('p');
axis equal;
%% simulation en grand temps
%on fait  1000 révolutions. Ca prend un peu de temps.
[t,y]=ode45(@(t,y) [y(2);-y(1)], [0,2000*pi],[0;1]);
%% tracé de la trajectoire
% le schéma dissipe de l'énergie. Les trajectoires ne se referment pas
% exactement et spiralent peu à peu vers l'origine.
%c'est pour cela que le graphe est épaissi.
figure;
plot(y(:,1),y(:,2));
axis equal;
xlabel('q');
ylabel('p');
 %% précision plus élevée
 % on augmente la précision du solveur
 % en espérant que les trajectoires se referment mieux.
 % C'est mieux mais si on zoome ce n'est pas parfait
 % et c'est long à calculer!
 options = odeset('RelTol', 1e-5);
 [t,y]=ode45(@(t,y) [y(2);-y(1)], [0,2000*pi],[0;1],options);
 figure;
 plot(y(:,1),y(:,2));
 xlabel('y');
 ylabel('dy/dt');
%% schéma Euler symplectique
% q_{n+1} = q_n + h p_{n+1} 
% p_{n+1} = p_n - h q_n
h= 2e-1;
t=(0:h:2000*pi);
N=length(t);
q=zeros(1,N);
p=zeros(1,N);
q(1)=0;
p(1)=1;
for k=1:N
    p(k+1) = p(k) - h*q(k);
    q(k+1) = q(k) + h*p(k+1);
end;
figure;
plot(q(:),p(:),'-');
axis equal;
xlabel('q');
ylabel('p');
hold on;
plot(cos( t ), sin( t ),'r');
legend('schéma Euler implicite-explicite', 'cercle unité');
% les trajectoires se referment
% et le schéma ne perd plus d'énergie.
% Il est d'ordre 1 donc peu précis ce qui explique
% que la trajectoire n'est pas parfaitement circulaire.