///
/// modele du pendule (amorti ou non)
/// avec force exterieure de type sinusoidale
///
k=evstr(x_dialog('valeur de k (frottement)','0'));
omega=evstr(x_dialog('valeur de omega (frequence forcée)','1.2'));
omega0=evstr(x_dialog('valeur de omega0 (frequence propre)','2'))
A=evstr(x_dialog('amplitude de la force exterieure','1'));
x0=evstr(x_dialog('angle initial du pendule (vitesse nulle)','0.5'));
//
Tmax=50;dt0=0.001;
deff('ydot=f(t,y)','ydot=[y(2);A*cos(omega*t)-k*y(2)-omega0^2*y(1)]');
//
//
//////  avec ode (en noir)
//
X0=[x0;0];
t=0:dt0:Tmax
x=ode(X0,0,t,f);
figure(0)
clf()
plot2d(t,x(1,:));
figure(1)
clf()
plot2d(x(1,:),x(2,:));

dt=evstr(x_dialog('choix du pas de temps avec Euler','0.001'));
//
//////  avec Euler (en bleu)
//
xe=[x0;0];Xetot=xe;
for i=2:length(t)
  xe=xe+dt*f(t(i-1),xe);
  Xetot=[Xetot,xe];
end
figure(0)
plot2d(t,Xetot(1,:),2);
legends(['avec ode','avec Euler'],[1,2])
//
figure(1)
plot2d(Xetot(1,:),Xetot(2,:),2);
legends(['avec ode','avec Euler'],[1,2])


