%% TP note 2024 
% pendule non lineaire
%y" = -sin y, y(0)=0,
% s'ecrit sous forme de système
% Y' = F(Y)
%avec Y = [y_1, y_2]^T
%F(Y) = [y_2;, -sin(y_1)]^T
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Question 1
% attention aux point virgules qui séparent les colonnes 
[t,y]=ode45(@(t,y) [y(2);-sin(y(1))], [0:0.01:64],[17*pi/18;0]);
%% tracé des  composantes
%
close all;% ferme toutes les figures.
figure;
plot(t,y(:,1),t,y(:,2));
legend('y','dy/dt');
xlabel('t');
title('schéma ode45');
%
%% tracé de la trajectoire
%dans l'espace des phases (y,y')
% 
figure;
plot(y(:,1),y(:,2),'.-');
xlabel('y');
ylabel('dy/dt');
title('plan de phase schéma ode45');
axis equal;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Question 2
% schéma Euler symplectique
% q_{n+1} = q_n + h p_{n+1} 
% p_{n+1} = p_n - h sin(q_n)
h= 1e-2;
t=(0:h:64);
N=length(t);
q=zeros(1,N);
p=zeros(1,N);
q(1)=17*pi/18;
p(1)=0;
for k=1:N-1
    p(k+1) = p(k) - h*sin(q(k));
    q(k+1) = q(k) + h*p(k+1);
end;
figure;
plot(t,q(:),t,p(:));
legend('y','dy/dt');
xlabel('t');
title('schéma implicite-explicite');
%
figure;
plot(q(:),p(:),'-');
axis equal;
xlabel('q');
ylabel('p');
title('plan de phase schéma implicite-explicite');
% les trajectoires se referment mieux qu'avec ode45
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Question 3
% conservation énergie
%test conservation de l'énergie
E= -cos(17*pi/18);
disp(' pour ode45  max |E(t) - E(0)|=');
disp(max(abs(0.5*y(:,2).^2 - cos(y(:,1))-E)));
%
disp(' pour le schema implicite-explicite   max |E(t) - E(0)|=');
disp(max(abs(0.5*p(:).^2 - cos(q(:))-E)));
figure;
plot(t,0.5*p.^2 - cos(q));
hold on
plot(t,0.5*y(:,2).^2 - cos(y(:,1)));
plot(t,E*ones(length(t),1));
xlabel('t');
ylabel('Energie');
legend('energie schema implicite-explicite','energie schema ode45','constante');
title('conservation energie');

% Le schéma implicite explicite conserve mieux l'énergie que ode45
% qui semble augmenter l'énergie de façon non physique
% cela explique que la trajectoire dans l'espace de phase se referme mieux
% avec le schéma implicite explicite.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Question 4
% calcul periode
function [val,isterm,dir] = retour(t,y)
% val est la valeur de la vitesse
% dir =-1 permet de chercher un maximum d'amplitude
% isterm = 0 permet de ne pas s'arreter au premier maximum
  val = y(2);
  isterm = 0;
  dir = -1;
end;
%option detection d'evenement
options = odeset('Events',@retour);
%resolution et calcul des instants de retour avec ode45
[t,y,te,ye,ie] = ode45(@(t,y) [y(2);-sin(y(1))], [0:0.01:64],[17*pi/18;0],options);
% instants de retour
te
disp('la periode moyenne pour ode45 est'); disp(mean(diff(te)));
% pour le schéma explicite implicite
% il suffit d'etudier les instants ou la vitesse s'annule en changeant de signe
% et en passant de positif à négatif de facon a detecter les maximums 
% d'amplitude
t_max = t(find(diff(sign(p)) < 0))
disp('la periode moyenne pour le schema implicite-explicite est');
 disp(mean(diff(t_max)));
 % La periode calculee est beaucoup plus  de la periode theorique avec 
 % le schéma implicite-explicite meme s'il est seulement d'ordre un.
 % on a vu en DM que c'est un schéma qui conserve les aires dans le plan de phase.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Question 5
%  condition initiale modifiee
% en fait la condition initiale 19*pi/20 correspondant à un angle de 171°
% ne  change rien aux resultats précédents. ODE45 de Matlab est plus précis que
% RK45 de  Python. Pourtant c'est en théorie le meme algorithme...voir les 
% complements plus bas, mais ne pas en tenir compte dans la correction
[t,y]=ode45(@(t,y) [y(2);-sin(y(1))], [0:0.01:64],[19*pi/20;0]);
figure;
plot(t,y(:,1),t,y(:,2));
legend('y','dy/dt');
xlabel('t');
title('schéma ode45, condition initiale y(0)= 19pi/20');
%
%% tracé de la trajectoire
%dans l'espace des phases (y,y')
% 
figure;
plot(y(:,1),y(:,2),'.-');
xlabel('q');
ylabel('p');
axis equal;
title('plan de phase schéma ode45, condition initiale y(0)= 19pi/20');
%% schéma Euler symplectique
% q_{n+1} = q_n + h p_{n+1} 
% p_{n+1} = p_n - h sin(q_n)
h= 1e-2;
t=(0:h:64);
N=length(t);
q=zeros(1,N);
p=zeros(1,N);
q(1)=19*pi/20;
p(1)=0;
for k=1:N-1
    p(k+1) = p(k) - h*sin(q(k));
    q(k+1) = q(k) + h*p(k+1);
end;
figure;
plot(t,q(:),t,p(:));
legend('y','dy/dt');
xlabel('t');
title('schéma implicite-explicite condition initiale y(0)= 19pi/20 ');
figure;
plot(q(:),p(:),'-');
axis equal;
xlabel('q');
ylabel('p');
title(' phase schéma implicite-explicite condition initiale y(0)= 19pi/20');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% COMPLEMENTS condition initiale 175.5° 
% c'est à partir de cet angle qu'on voit des aberrations.
[t,y]=ode45(@(t,y) [y(2);-sin(y(1))], [0:0.01:128],[175.5*pi/180;0]);
figure;
plot(t,y(:,1),t,y(:,2));
legend('y','dy/dt');
xlabel('t');
title('schéma ode45, condition initiale y(0)= 175.5°');
%
%% tracé de la trajectoire
%dans l'espace des phases (y,y')
% 
figure;
plot(y(:,1),y(:,2),'.-');
xlabel('q');
ylabel('p');
axis equal;
title('plan de phase schéma ode45, condition initiale y(0)= 175.5°');
%% schéma Euler symplectique
% q_{n+1} = q_n + h p_{n+1} 
% p_{n+1} = p_n - h sin(q_n)
h= 1e-2;
t=(0:h:128);
N=length(t);
q=zeros(1,N);
p=zeros(1,N);
q(1)=19*pi/20;
p(1)=0;
for k=1:N-1
    p(k+1) = p(k) - h*sin(q(k));
    q(k+1) = q(k) + h*p(k+1);
end;
figure;
plot(q(:),p(:),'-');
axis equal;
xlabel('q');
ylabel('p');
title(' phase schéma implicite-explicite condition initiale y(0)= 175.5°');
% le schéma ode45 donne des resultats aberrants, après deux periodes le pendule
% fait un tour complet  sur lui-même et continue de tourner sans s'arrêter.
% Le pendule gagne de l'énergie au bout de deux tours ???
%  Le schéma implicite-explicite n'a pas ce comportement aberrant. 
% Il conserve mieux l'energie.