% ressort02.m
% fd  6 avril 2005, 08 janvier 2006, 09, 10 mars 2006, 25 fevrier 2025  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% title ('Integration d un oscillateur harmonique amorti') ;

tmax=4*pi 
n=floor(10*tmax) + 1
dt=.1
mm=1. , cc=.05 , kk=1.
  
t=zeros(n+1,1);
x=zeros(n+1,1);
y=zeros(n+1,1);

x0 = 1 ;
y0 = 0 ;
xx = x0 ;
yy = y0 ;
ga = -mm \ (cc * yy + kk*xx)  ;
ss = mm + .5*dt*cc + .25*dt*dt*kk ; 
for i=1:n;
t(i+1)=t(i)+dt;
end ;


x(1) = x0 ; y(1) = y0 ; 
for i=1:n;
xt = xx + dt*yy + .25*dt*dt*ga ; 
yt = yy + .5*dt * ga ; 
ga = ss \ ( - kk*xt - cc*yt ) ;
xn = xt + .25*dt*dt*ga ; 
yn = yt + .5*dt*ga ; 
x(i+1) = xn ; 
y(i+1:i+1,1) = yn ; 
xx = x(i+1) ; 
yy = y(i+1) ; 
end ;


temps_final = t(n+1:n+1)

hold off
plot(x(1,1) , y(1,1), 'ob;depart;') ; hold on 
plot(x(:,1) , y(:,1), '*k;trajectoire;') ; hold on 
plot(x(n+1,1) , y(n+1,1), 'or;arrivee;') ; hold on 
title ('oscillateur harmonique amorti ; plan de phase') ;
axis ( [ -1.1   , 1.1  , -1.1 , 1.1 ]) ;
pause (1.5) 
print('02b-plan-de-phase.jpg', '-djpg', '-S891,630') ;

hold off 
xmax=tmax+.1
plot(t(:,1) , x(:,1), '*b;position;') ; hold on 
plot(t(:,1) , y(:,1), 'or;impulsion;') ; 
axis ( [ -.1   , xmax  , -1.1 , 1.1 ]) ;
title ('oscillateur harmonique amorti ; evolution temporelle') ;
print('02a-instationnaire.jpg', '-djpg', '-S891,630') ; 














