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

mm=1 ; cc=.1 ; kk=1 ;


mmcible = mm
cccible = cc
kkcible = kk 
  
tmax=4*pi 
n=256
dt=tmax/n ; 
%
%   parametres initiaux            
% 
masse=1 ; amort=1  ; raideurini=.1 ; raideur = raideurini ;


mminitiale = masse
ccinitiale = amort
kkinitiale = raideurini

%
%                   nombre d iterations du gradient  
nitmax = 33 ; %% 120 ; 
%                  parametre de sous-relaxation
roro   = .01 ; 
%
t=zeros(n+1,1); 
%
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 ) ; 
ga = -mm \ ( kk*xxx + cc*yyy )  ;
xt = xxx + dt*yyy + .25*dt*dt*ga ;
yt = yyy + .5*dt * ga ; 
ga = ss \ ( - kk*xt - cc*yt ) ;
xn = xt + .25*dt*dt*ga ; 
yn = yt + .5*dt*ga ; 
x(i+1:i+1,1)=xn ; 
y(i+1:i+1,1)=yn ; 
xxx=xn ; 
yyy=yn ; 
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;position exacte;' , t(:,1) , y(:,1), '+m;vitesse  exacte;') ; 
axis ( [ -.1   , tmax+.1  , -1.1 , 1.1 ]) ;
title ('oscillateur harmonique amorti ; solution cible') ;
print('05a-etat-dynamique-cible', '-djpg', '-S891,630');
pause (1) 



%                                         boucle en les parametres 
%
kkas = zeros (nitmax+1,1) ; ccas = kkas ; zlogerr = kkas ; 
jiter = kkas ; 
% 
lll=0 ; 
jiter  (lll+1:lll+1,1)  = lll ; 
kkas (lll+1:lll+1,1)    = raideur ; 
ccas (lll+1:lll+1,1)    = amort ; 
zlogerr (lll+1:lll+1,1) = 1. ; 
%
for lll = 1:nitmax ;            %%%  iterations de gradient 
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;                      %%%  iterations en temps du probleme direct  
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 ; 

if (lll==1) ; 
plot( t(:,1) , xx(:,1), '*;x   position;' ,  t(:,1) , yy(:,1), '+b;y : vitesse;'  ) ; 
axis ( [ -.1   , tmax+.1  , -.15 , 1.1 ]) ;
title ('oscillateur harmonique amorti ; etat dynamique amorti au debut des iterations de gradient') ;
print('05b-etat-debut-iterations-gradient', '-djpg', '-S891,630');
pause (1) 
end  ; 

if (lll==nitmax) ; 
plot(t(:,1) , x(:,1), ';position exacte;' , t(:,1) , y(:,1), ';vitesse  exacte;' , ...
     t(:,1) , xx(:,1), '*;position approchee;' ,  t(:,1) , yy(:,1), '+b;vitesse approchee;'  ) ; 
axis ( [ -.1   , tmax+.1  , -1.1 , 1.1 ]) ;
title ('oscillateur harmonique amorti ; etat dynamique a la fin des iterations de gradient') ;
print('05e-etat-fin-iterations-gradient', '-djpg', '-S891,630');
pause (1) 
end  ; 


err = xx - x ; 
zer = 0. ;
for i=1:n+1;
zer = max (zer , abs (err(i:i,1))) ; 
end ; 
zer = 1.2 * zer ;
%  axis ( [ -.1   , tmax+.1   , -zer , zer ]) ;
%  plot(t(:,1) , err(:,1), '*;erreur;') ; 
%  pause 

%                                                 equation adjointe
ss = masse + .5*dt*amort + .25*dt*dt*raideur ; 
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;           %%%  iterations en temps du probleme adjoint   
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) + raideur*ro ; 
qt    = qq(jj+1:jj+1,1) + amort*ro ; 
rot   = - masse \ (.25*dt*dt*pt + .5*dt*qt ) ;
pp(jj:jj,1) = pt + dt*err(jj+1:jj+1,1) + raideur*rot ; 
qq(jj:jj,1) = qt + dt*pt               + amort*rot ; 
ymax = max (ymax,abs(pp(jj:jj,1))) ; 
ymax = max (ymax,abs(qq(jj:jj,1))) ; 
end ;

ymax = ymax + .5 ; 
if (lll==1) ; 
plot( td(:,1) , pp(:,1), '*;etat adjoint;' , td(:,1) , qq(:,1), '+b;second etat adjoint;'  ) ; 
axis ( [ -.1   , tmax+.3  , -.15*ymax  , ymax  ] ) ;
title ('oscillateur harmonique amorti ; etat adjoint initial') ;
print('05c-etat-adjoint-initial', '-djpg', '-S891,630');
pause (1) 
end ; 

%                          gradient 
grkk  = 0 ; grcc = 0 ; 
for j = 1:n ;
xxx   = x(j:j,1) ; 
yyy   = y(j:j,1) ; 
ga    = -masse \ ( raideur*xxx + amort*yyy )  ;
xt    = xxx + dt*yyy + .25*dt*dt*ga ;
yt    = yyy + .5*dt * ga ; 
gap   = ss \ ( - raideur*xt - amort*yt ) ;
ro    = -ss \ (.25*dt*dt*pp(j:j,1) + .5*dt*qq(j:j,1) ) ; 
pt    = pp(j:j,1) + raideur*ro ; 
qt    = qq(j:j,1) + amort*ro ; 
% vieille version avec bug, fd le 14 mai 
% pt    = pp(jj+1:jj+1,1) + raideur*ro ; 
% qt    = qq(jj+1:jj+1,1) + amort*ro ; 
% vieille version avec bug, fd le 14 mai 
rot   = - masse \ (.25*dt*dt*pt + .5*dt*qt ) ;
grkk  = grkk + ro*(xt+dt*dt*gap/4) + rot* xxx ; 
grcc  = grcc + ro*(yt+dt*gap/2)    + rot* yyy ; 
end ; 

raideur = raideur - roro * grkk ; 
amort   = amort   - roro * grcc ; 

jiter  (lll+1:lll+1,1)  = lll ; 
kkas (lll+1:lll+1,1)    = raideur ; 
ccas (lll+1:lll+1,1)    = amort ; 
zlogerr (lll+1:lll+1,1) = log (zer) / log (10) ; 

end   %%%% fin des iterations de gradient




%%% ppf = zeros(n,1); qqf = ppf ; 
%%% for j = 1:n ;
%%% ppf(j:j,1) = pt(j:j,1) ; qqf(j:j,1) = qt(j:j,1) ; end ; 


hold off 
plot( td(:,1) , pp(:,1), '*;etat adjoint final;' , td(:,1) , qq(:,1), '+b;second etat adjoint final;'  ) ;
axis ( [ -.1   , tmax+.3  , -.010  , .010  ] ) ;
title ('oscillateur harmonique amorti ; etat adjoint final') ;
print('05d-etat-adjoint-final', '-djpg', '-S891,630');
pause (1) 



zzmm = nitmax +.5 ; 
plot ( jiter , zlogerr , '@m;log de l erreur;' )
if nitmax==120 ; axis ( [ -.5   , zzmm  , -14  , 1.1  ] ) ; else ;
axis ( [ -.5   , zzmm  , -3  , 1.1  ] ) ; end 
title ('oscillateur harmonique amorti ; erreur au cours des iterations de gradient') ;
print('05f-erreur-iterations-gradient.jpg', '-djpg', '-S891,630');
pause (1) 


zzmm = max (kk , raideurini) +.1 ; 
plot (jiter , kkas, '*;raideur;' , jiter , ccas , '+b;amortissement;' , ...
      [0 , nitmax] ,  [kk , kk]  , 'r;raideur cible;' , ...
      [0 , nitmax] ,  [cc , cc]  , 'b;amortissement cible;' ) 
axis ( [ -.5   , nitmax+.5  , -.05  , zzmm+.13  ] ) ;
title ('oscillateur harmonique amorti ; parametres du ressort au cours des iterations') ;
print('05g-parametres-iterations-gradient.jpg', '-djpg', '-S891,630');
