% fil2.m
% fd, 21 novembre 2006, 03 mars 2025. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Newton sur le probleme de point selle 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear ; format long 

   
hold off 
%%% title ('Chainette numerique. Algorithme de Newton sur les equations \
%%% de point selle.') ;

jdif = 0
  
if jdif==0 ; longth = 1.6 ; niterNewton = 9 , ro = 1. , miny = -.21 ; end ; 
if jdif==1 ; longth = 4.831525 ; niterNewton = 20 , ro = .7 ; miny = -1.9 ; end 

%%  longth = 3.5  
%%  longth = 4.831525  %% 1.6  , %  longth = 4.831525 
init = 0 ;  if (init==0) ; 
n=200 ; dx=1/n ;  x = zeros(1,n+1) ;  y=x ;  
for j = 1 : n+1 ;    %  initialisation droite 
xx = (j-1)*dx ; x(1,j:j) = xx ; y(1,j:j) = 1 - xx - 7.5*xx*(1-xx)  ; end ; 
plot ( x ,  y , '*b;condition initiale : parabole;' , x,y,'r;;' ) ;
title ('chainette numerique ; algorithme de Newton pour les equations de point selle') ;
axis ( [ -.01 , 1.01 , -1.55  , 1.05 ] ) ; 
print('02a-chainette-numerique-Newton', '-djpg','-S891,630'); pause (.5)   ;  

lambda = 0.45   ;    %  ok pour longth = 1.67 
lambda = 1.26   ;    %  ok pour longth = 3.5 
lambda = 1.75   ;    %  ok pour longth = 4.5
lambda = 1.9
  
%  Tactique. tel que partir de longth = 3. 
%  Puis avancer la valeur de longth de .5 par continuation
 else ;  y=zz ; end ; 

%%%%%%%%%
for jj = 1 : niterNewton ;   %  iterations de Newton
cout  = 0 ;  long = 0 ; for j = 1 : n ; 
yy = y(1,j:j) ; yyp1 = y(1,j+1:j+1) ; 
delta = yyp1-yy ; pp = .5*(yy+yyp1) + lambda ; 
rr = sqrt ( 1 + (delta/dx)**2 ) ; 
cout  = cout  + pp*rr ; long = long + rr ;  end ; 
cout = dx * cout  ;   long = dx * long ;  
%                                   gradients
gr  = zeros(1,n+1) ; grco  = zeros(1,n+1) ;  
for j = 1 : n ;  yy = y(1,j:j) ; yyp1 = y(1,j+1:j+1) ; 
delta = yyp1-yy ; pp = .5*(yy+yyp1) + lambda ;  
rr = sqrt ( 1 + (delta/dx)**2 ) ; 
gr (1,j:j)       = gr (1,j:j)       + .5*rr*dx - pp*delta / (rr*dx) ; 
gr (1,j+1:j+1)   = gr (1,j+1:j+1)   + .5*rr*dx + pp*delta / (rr*dx) ; 
grco (1,j:j)     = grco (1,j:j)     - delta / (rr*dx) ; 
grco (1,j+1:j+1) = grco (1,j+1:j+1) + delta / (rr*dx) ; 
end ; 
%                                   Hessienne
matelt = zeros(2,2,n) ;        %    matrices elementaires 
for j = 1 : n ;  yy = y(1,j:j) ; yyp1 = y(1,j+1:j+1) ; 
delta    = yyp1-yy ; 
pp       = .5*(yy+yyp1)  + lambda ; 
rr       = sqrt ( 1 + (delta/dx)**2 ) ; 
dshr     = delta / (dx*rr) ; 
pshr     =  pp / (rr*dx) ; 
pd2sh3r3 =  pp*delta**2 / ((dx**3)*(rr**3)) ; 
matelt(1:1,1:1,j:j) =  - dshr + pshr - pd2sh3r3 ; 
matelt(1:1,2:2,j:j) =         - pshr + pd2sh3r3 ; 
matelt(2:2,1:1,j:j) =    matelt(1:1,2:2,j:j) ; 
matelt(2:2,2:2,j:j) =    dshr + pshr - pd2sh3r3 ; 
end ; 
%                                   assemblage  
aa = zeros(1,n-1) ; bb = aa ; cc = aa ; 
aa(1,1:1) = matelt(2:2,2:2,1:1) ;   
for j = 2 : n-1 ; 
aa(1,j-1:j-1) = aa(1,j-1:j-1) +  matelt(1:1,1:1,j:j) ;  
bb(1,j:j)     = bb(1,j:j)     +  matelt(1:1,2:2,j:j) ; 
cc(1,j:j)     = cc(1,j:j)     +  matelt(2:2,1:1,j:j) ;   
aa(1,j:j)     = aa(1,j:j)     +  matelt(2:2,2:2,j:j) ;    
end ;       
aa(1,n-1:n-1) = aa(1,n-1:n-1) +  matelt(1:1,1:1,n:n) ; 
%                            resolution des systemes lineaires 
sm = aa ;  
for j = 1 : n-1 ;  sm(1,j:j) = gr (1,j+1:j+1) ; end ; 
%                                   resolution du systeme lineaire 
xx=tridia(n-1,aa,bb,cc,sm) ;       
for j = 1 : n-1 ;  sm(1,j:j) = grco (1,j+1:j+1) ; end ; 
xxco=tridia(n-1,aa,bb,cc,sm) ;      
%                                   produits scalaires
prod = 0 ;  ppss = 0 ;  
for j = 1 : n-1 ; 
prod = prod +  grco (1,j+1:j+1)*xx(1,j:j) ; 
ppss = ppss +  grco (1,j+1:j+1)*xxco(1,j:j) ; end ;
jj , dlam = ( long - longth - prod ) / ppss 
%                                         incrementation 
for j = 2 : n ; 
y(1,j:j) = y(1,j:j) - ro* (xx (1,j-1:j-1) + dlam*xxco(1,j-1:j-1)) ; end ;

if jj==1 ; plot ( x ,  y , '*b;iteration 1;' , x,y,'r;;' ) ;
title ('chainette numerique ; algorithme de Newton pour les equations de point selle') ;
axis ( [ -.01 , 1.01 , miny  , 1.05 ] ) ; 
print('02b-chainette-numerique-Newton', '-djpg','-S891,630'); pause  (.5) ; end 
if jj==2 ; plot ( x ,  y , '*b;iteration 2;' , x,y,'r;;' ) ;
title ('chainette numerique ; algorithme de Newton pour les equations de point selle') ;
axis ( [ -.01 , 1.01 , miny  , 1.05 ] ) ;  
print('02c-chainette-numerique-Newton', '-djpg','-S891,630'); pause  (.5) ; end
if jj==3 ; plot ( x ,  y , '*b;iteration 3;' , x,y,'r;;' ) ; 
title ('chainette numerique ; algorithme de Newton pour les equations de point selle') ;
axis ( [ -.01 , 1.01 , miny  , 1.05 ] ) ;  
print('02d-chainette-numerique-Newton', '-djpg','-S891,630'); pause  (.5) ; end  
if jj==4 ; plot ( x ,  y , '*b;iteration 4;' , x,y,'r;;' ) ; 
title ('chainette numerique ; algorithme de Newton pour les equations de point selle') ;
axis ( [ -.01 , 1.01 , miny  , 1.05 ] ) ;  
print('02e-chainette-numerique-Newton', '-djpg','-S891,630'); pause  (.5) ; end     
if jj==5 ; plot ( x ,  y , '*b;iteration 5;' , x,y,'r;;' ) ; 
title ('chainette numerique ; algorithme de Newton pour les equations de point selle') ;
axis ( [ -.01 , 1.01 , miny  , 1.05 ] ) ; 
print('02f-chainette-numerique-Newton', '-djpg','-S891,630'); pause  (.5) ; end
if jj==6 ; plot ( x ,  y , '*b;iteration 6;' , x,y,'r;;' ) ; 
title ('chainette numerique ; algorithme de Newton pour les equations de point selle') ;
axis ( [ -.01 , 1.01 , miny  , 1.05 ] ) ; 
print('02g-chainette-numerique-Newton', '-djpg','-S891,630'); pause  (.5) ; end
if jj>=7 ; plot ( x ,  y , '*b;iterations;' , x,y,'r;;' ) ; 
title ('chainette numerique ; algorithme de Newton pour les equations de point selle') ;
axis ( [ -.01 , 1.01 , miny  , 1.05 ] ) ; pause  (.001) ; end



lambda = lambda + ro*dlam ; 
long ;  lambda ; pause (.5)  , end ; 
%                                   calcul de la valeur de theta  
for j = 2 : n ; 
yym1 = y(1,j-1:j-1) ; yy = y(1,j:j) ; yyp1 = y(1,j+1:j+1) ;  
if ((yy<=yym1)&&(yy<=yyp1)) jth = j ; end ; end ; 
theta = (jth-1) * dx 
%
zz=y ;  

long
lambda

plot ( x ,  y , '*b;;' , x,y,'r;;' ) ; 
%%% axis ( [ -.01 , 1.01 , -.21  , 1.05 ] ) ;
axis ( [ -.01 , 1.01 , miny  , 1.05 ] ) ; 
title ('chainette numerique ; algorithme de Newton sur les equations de point selle') ;
print('02h-chainette-numerique-Newton-finale', '-djpg','-S891,630'); 
