% ressort03.m
% fd  6 avril 2005, 08 janvier 2006, 09, 10  mars 2006
   % fd  27 fevrier, 01 mars 2025 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% title ('precision du Newmark dans le cas non amorti') ;
mm=1 ; cc=0 ; kk=1 ; 

tmax=4*pi 
m=6
pas=zeros(m,1);
err=zeros(m,1);
%
n=32

  hold off 
  for j=1:m;
dt=tmax/n
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,1)

%% if (j<=3) ; hold on ; else ; hold off ; end ; 
if (n==32) ; plot(t(:,1) , x(:,1), '*k;n=32;') ; hold on ; end ;
if (n==64) ; plot(t(:,1) , x(:,1), 'ob;n=64;') ; hold on ; end ;
if (n==128) ; plot(t(:,1) , x(:,1), 'dr;n=128;') ; hold on ; end ;
if (n>=129) ; plot(t(:,1) , x(:,1), '*k;;') ; hold off ; end ;
axis ( [ -.2   , tmax+.2 , -1.1 , 1.39]) ;
title ('schema de Newmark pour un oscillateur harmonique') ;
if (j==3) ;
print('03-schema-Newmark.jpg', '-djpg','-S891,630');
print('03-schema-Newmark.pdf', '-dpdf','-S891,630');
end 
pause (.3) 
%
%                       erreur
pas(j:j,1)=n ; 
err(j:j,1)=log ( abs( x(n+1:n+1,1) - 1 ) + abs( y(n+1:n+1,1) ) ) / log(2)  ;
     n=2*n;
end

pause (2) 
hold off 
x=zeros(m,1); x =log(pas(:,1)) / log(2) ;
y=zeros(m,1); y= err(:,1) ;
     xmax= 22 ; 
plot( x , y ,'*;logarithme des ereurs;') ;
axis ([ 4.9 ,  m+4+.1  ,  -13.6 , -1.9 ]) ;
title ('precision du schema de Newmark dans le cas non amorti') ;
print('03-precision-Newmark.jpg', '-djpg', '-S891,630');
print('03-precision-Newmark.pdf', '-dpdf', '-S891,630');



%
%               droite des moindres carres
alpha=codre(m,x,y) 





