clear;
// Calcul 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 = 1;    // entier qui caracterise la donnee au bord
eps = .1; // parametre de penalisation
N = 16;   // Entier caracterisant la discretisation
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. 

////////////////    ***************      //////////////     
//                                                         /1                                        
// A] ASSEMBLAGE de la matrice A telle que (A R, S) =     |  [ rho sigma d^2/ r + rho' sigma' r ] dr  = (rho, sigma)_H
//                                                        /0
// ETAPE 1)  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;

// ETAPE 2) 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);

// ETAPE 3) Construction de la matrice creuse A 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]);

// Exemple : calcul du  produit scalaire (rho, sigma)_H pour rho( r) = r^3 et sigma2( r) = r^2 
	R = rr.^3; S = rr.^2;
        prodscal = R'*(A*S)  // (valeur exacte 7/5)

////////////////    ***************      /////////////// 
//                                             /1                                        
// B] Calcul du terme de penalisation    Int = |   (1 - rho^2)^2 r dr 
//                                            /0
// exemple dans le cas rho ( r ) = r^2
	R = rr^2;
	Integrales = zeros(N,1); // vecteur contenant la valeur de l'integrale sur l'intervalle [(i-1)h, i h]   
	Integrales(1) = h*(1-(R(1)/2)^2)^2*r(1);
	Integrales(2:N) = h*((1-(.5*R(2:N) + .5*R(1:N-1)).^2).^2).*r(2:N)';
	Int = sum (Integrales)   // (valeur exacte 4/15)

////////////////    ***************      /////////////// 
//           
// C] Calcul de l'energie pour rho ( r ) = r^2
//
	K = 1/2*(R'*(A*R)) + 1/(4*eps^2)*Int