% ressort04.m
% fd  09, 10, 14  mars 2006, 01 mars 2025 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% title ('Assimilation variationnelle pour un oscillateur harmonique amorti') ;

mm=1 ; cc=.01 ; kk=1 ; 
tmax=4*pi 

n=32
dt=tmax/n

t=zeros(n+1,1); 
uu=zeros(n+1,2) ; 


x=t ; y=x ; xx=x ; yy = xx ; td = td=zeros(n,1);
x0 = 1 ;
y0 = 0 ;
xxx = x0 ;
yyy = y0 ;
x(1:1,1)=xxx ; 
y(1:1,1)=yyy ; 
ss = mm + .5*dt*cc + .25*dt*dt*kk ; 
for i=1:n;
uu = direct ( mm , cc , kk , ss , xxx , yyy , dt ) ; 
xxx=uu(1:1,1) ;
yyy=uu(2:2,1) ; 
x(i+1:i+1,1)=xxx ; 
y(i+1:i+1,1)=yyy ; 
end ; 

for i=1:n;
t(i+1)=t(i)+dt;
td(i:i,1) = (i-.5)*dt ; 
end ;

hold off
plot(t(:,1) , x(:,1), 'ob;x exact;' , t(:,1) , x(:,1), 'b;;' , ...
     t(:,1) , y(:,1), '+m;y exact;', t(:,1) , y(:,1), 'm;;') ; 
axis ( [ -.1   , tmax+.1  , -1.1 , 1.3 ]) ;
title ('oscillateur harmonique amorti ; schema de Newmark') ;
print('04-probleme-direct.jpg', '-djpg','-S891,630');
%%% print('04-probleme-direct.pdf', '-dpdf');
pause (2) 
%
%               initialisation
masse=1 ; amort=.005 ; raideur=.5 ; 



%                              boucle en les parametres 
for lll = 1:10

xxx = x0 ;
yyy = y0 ;
xx(1:1,1)=xxx ; 
yy(1:1,1)=yyy ; 
ssiter = masse + .5*dt*amort + .25*dt*dt*raideur ; 
for i=1:n;
uu = direct ( masse , amort , raideur , ssiter , xxx , yyy , dt ) ; 
xxx=uu(1:1,1) ;
yyy=uu(2:2,1) ; 
xx(i+1:i+1,1)=xxx ; 
yy(i+1:i+1,1)=yyy ; 
end ; 

%% plot( t(:,1) , xx(:,1), '*;x approche;' ,  t(:,1) , yy(:,1), '+b;y approche;'  ) ; 
%% pause  (1) 

err = xx - x ; 
axis ( [ -.1   , tmax+ .1   , -2.1 , 2.1 ]) ;
plot(t(:,1) , err(:,1), '*;erreur;') ; 

%                          equation adjointe
%%%% cette ligne st un bug. fd le 14 mai 06
ss = mm + .5*dt*cc + .25*dt*dt*kk ; 
%%%% cette ligne st un bug. fd le 14 mai 06
pp = zeros(n,1); qq = pp ; 
pp(n:n,1) = .5 * dt * err(n+1:n+1,1) ; 
qq(n:n,1) = 0 ; 
ymax = abs(pp(n:n,1)) ; 
ymax = max (ymax,abs(qq(n:n,1))) ; 
for j = 1:n-1;
jj    = n-j ; 
ro    = -ss \ (.25*dt*dt*pp(jj+1:jj+1,1) + .5*dt*qq(jj+1:jj+1,1) ) ; 
pt    = pp(jj+1:jj+1,1) + kk*ro ; 
qt    = qq(jj+1:jj+1,1) + cc*ro ; 
rot   = - mm \ (.25*dt*dt*pt + .5*dt*qt ) ;
pp(jj:jj,1) = pt + dt*err(jj+1:jj+1,1) + kk*rot ; 
qq(jj:jj,1) = qt + dt*pt + cc*rot ; 
ymax = max (ymax,abs(pp(jj:jj,1))) ; 
ymax = max (ymax,abs(qq(jj:jj,1))) ; 
end ;

hold off
plot( td(:,1) , pp(:,1), '*;premier etat adjoint;' , td(:,1) , qq(:,1), '+b;second etat adjoint;'  ) ;
axis ( [ -.1   , tmax+.1  , -ymax-.5 , ymax+.5 ]) ;
pause(.5) 


%                          gradient 
gr    = 0 ; 
for j = 1:n ;
xxx   = x(j:j,1) ; 
yyy   = y(j:j,1) ; 
ga    = -masse \ ( kk*xxx + cc*yyy )  ;
xt    = xxx + dt*yyy + .25*dt*dt*ga ;
ro    = -ss \ (.25*dt*dt*pp(j:j,1) + .5*dt*qq(j:j,1) ) ; 
gr    = gr + ro*xt  ; 
end ; 

format long 
raideur = raideur - .03 * gr  

end 
