
% Sean Costello et Sandra Brella, 2014
%-------------------------------------------------------------------------
% SD demo 1: Un exemple de simulation du pendule de la figure 2.1
% QUESTION: Simuler le pendule de Fig 2.1, en partant de la position
% a) horizontale
% b) verticale (vers le haut). Simuler le systeme pendant au moins 20 s.
%-------------------------------------------------------------------------

% Nettoyage de l'espace de travail
clear all; close all; clc;

%-------------------------------------------------------------------------
% Definition des parametres physiques 
%-------------------------------------------------------------------------
g = 9.81;   % Constante de Gravitation de la Terre [m/s^2
l = 1;      % Longueur Pendule [m]

%-------------------------------------------------------------------------
% L'equation dynamique doit etre une fonction sous la forme dx/dt = f(t,x)
% Etats: x = [theta, dtheta], angles [rad]
%-------------------------------------------------------------------------
mydyn = @(t,x) [x(2); -g*sin(x(1))/l]; 

%-------------------------------------------------------------------------
% Definition des conditions initiales 
%-------------------------------------------------------------------------
x0 = [pi/2;0]; %horizontal
%x0 = [pi;0]; %vertical (vers le haut)

%-------------------------------------------------------------------------
% Integration Numerique de l'equation dynamique du Pendule
%-------------------------------------------------------------------------
odeoptions = odeset('RelTol',1e-6,'AbsTol',1e-6,'Stats','on'); 
tspan = [0,20]; %horizon temporel de la simulation en [s]
[tsim, xsim] = ode45(mydyn,tspan,x0,odeoptions); 

%-------------------------------------------------------------------------
% Affichage des resultats
%-------------------------------------------------------------------------
figure('Name','theta0 = pi/2'),
plot(tsim, xsim)
legend('theta (rad)','dtheta (rad/s)')
xlabel('t (sec)')

%-------------------------------------------------------------------------
% Explication des resultats
%-------------------------------------------------------------------------
% Pour a), mouvement harmonique non amorti.
%
% Pour le cas b), le resultat n'est pas conforme a la solution analytique, 
% laquelle voudrait que le pendule reste a l'equilibre (instable)
% en position verticale. En effet, la miniscule erreur numerique pendant
% l'integration (et une telle erreur est inevitable) suffit pour perturber
% l'equilibre instable du pendule.

