clear;
close();
// Minimisation de l'Energie de Ginzburg Landau d une fonction radiale
//  u = rho( r ) (cos(d theta), sin(d theta)).
//  L'Energie est definie par
//
//             / 1
//   K(rho) = |     [ 1/2 rho^2 d^2/r  + 1/2 (rho')^2 + 1/(4*eps^2) (1 - rho^2)^2 ] dr   
//            / 0
//


d = 2;    // entier qui caracterise la donnee au bord
eps = .05; // parametre de penalisation
N = 1000;   // Entier caracterisant la discretisation
itermax = 1000; // Nbre d'Iteration maximal pour le gradient
h = 1/N;  // pas de la discretisation
r = [0:N-1]*h + h/2; // liste des milieux des mailles
rr = [1:N]'*h; // liste des abscisses des degres de liberte. 

p = 10*eps^2;

////////////////    ***************      //////////////     
//                                                         /1                                        
// ASSEMBLAGE de la matrice A telle que (A R_1, R_2) =  |  [ d^2 *rho_1 rho_2 / r + rho_1' rho_2' r ] dr  = (rho_1, rho_1)_H
//                                                        /0
//Calcul des coefficients diagonaux
	c_diag = zeros(1,N);
	c_diag(1:N-1) = d^2*h/4*( r(1:N-1).^(-1) + r(2:N).^(-1) ) + 1/h*(r(1:N-1) + r(2:N));
	c_diag (N) = d^2*h/(4*r(N)) + r(N)/h;
//Calcul des coefficients sous et sur diagonaux
	c_pasdiag = zeros(1,N-1);
	c_pasdiag(1:N-1) = (d^2*h/4)*r(2:N).^(-1) - 1/h*r(2:N);
//Construction des matrices creuses A et A0 definie par 
//           A(i,i) = c_diag(i),    A(i,i+1) = A(i+1,i) = c_pasdiag(i),    A(i,j) = 0  si |j-i|>1. 
	A = sparse([1:N ; 1:N]', c_diag, [N,N]) + sparse([1:N-1 ; 2:N]', c_pasdiag, [N,N]) + sparse([2:N ; 1:N-1]', c_pasdiag, [N,N]);
	A0 = sparse([1:N-1 ; 1:N-1]', c_diag(1:N-1), [N-1,N-1]) + sparse([1:N-2 ; 2:N-1]', c_pasdiag(1:N-2), [N-1,N-1]) + sparse([2:N-1 ; 1:N-2]', c_pasdiag(1:N-2), [N-1,N-1]);
///////////////    *****************    //////////////////
// 
// INITIALISATION 
R = rr.^2;
iter = 0; 
err = 1;
ancienK = 1/2*(R'*A*R) + h/(4*eps^2)*sum(((1-(.5*R(2:N) + .5*R(1:N-1)).^2).^2).*r(2:N)');
Koo = ancienK;
B=zeros(N-1,1);


///////////////    *****************    //////////////////
// 
// BOUCLE DU GRADIENT
while (abs(err) > 1e-5) & (iter< itermax)
	iter = iter +1;
	// calcul du gradient T (calcul de B et A*R)
	B(1) =  .5*(1- (.5*R(1))^2)*R(1)*r(1);
	B(2:N-1) =  (1- (.5*R(2:N-1) + .5*R(1:N-2)).^2).*(.5*R(2:N-1) + .5*R(1:N-2)).*r(2:N-1)';
	B = B + (1- (.5*R(2:N) + .5*R(1:N-1) ).^2).*(.5*R(2:N) + .5*R(1:N-1)).*r(2:N)';
	B =  (-.5*h/eps^2)*B;
        AR = A*R;
	T = A0\(AR(1:N-1) + B);
	R(1:N-1) = R(1:N-1) - p*T;
	K = 1/2*(R'*A*R) + h/(4*eps^2)*sum(((1-(.5*R(2:N) + .5*R(1:N-1)).^2).^2).*r(2:N)')
        err = (ancienK - K)/Koo; ancienK = K;
end;
//iter
//K
plot(rr,R);