program ep95_246;

{
Epreuve informatique de l'Ecole Polytechnique           95.246
--------------------------------------------------------------



	Soit Y=  [ y11  y12  ]
                 [ y12  y22  ]
          une matrice carrée d'ordre 2 à coefficients réels.

On cherche à définir une autre matrice X qui soit sa racine carrée, de manière à
donner un sens à la notation  X=sqrt(Y)  par analogie avec le raisonnement
usuel dans R.


1° 	En notant X²  le produit matriciel de la matrice X par elle-même,
on considère l'équation :


                      X²=Y

Discuter l'existence et le nombre de solutions X en fonction de Y.
Ecrire un programme Pascal qui recherche et calcule automatiquement ces 
solutions, en distinguant bien les différents cas possibles. Vérifier
le bon fonctionnement du programme dans chacun des quatre cas suivants :


1) Y=  [ 1 0 ]   2) Y = [ 4 3 ]   3)   Y = [ 4 -3 ]  4)   Y = [0  0 ]
       [ 0 1 ]          [ 3 4 ]            [ 3  4 ]           [0  2 ]


2°  	On cherche à résoudre l'équation matricielle suivante  :

           AX²+BX+C=0

dans le cas où A, B et C sont des matrices carrées à  2x2  coefficient réels, 
et où X est une matrice inconnue.

Que peut-on dire a priori des solutions de cette équation ?

Ecrire un programme Pascal qui traite le problème en fonction des matrices A, B et C.


Application au cas suivant  :

		
   A=  [ 2  0 ]        B = [ -1  -2 ]   C  =  [ 0  1 ]
       [ 0  2 ]            [ -1  -1 ]         [ 1  0 ]



:--------------------------------------------------------------:
:    Imprimer tous les résultats                               :
:     en indiquant chaque fois à quoi ils correspondent        :
:--------------------------------------------------------------:
                            -=-=-=-
 }



uses
  crt, modubase;



TYPE
    MyReal=Double;
    matrice=array[1..2,1..2] of Double;

Const
     reduction = 0.3;
     n=2;




VAR
     Car         : Char ;
     Question    : Char;

Function determinant(m:matrice):Myreal;
begin
     determinant:=m[1,1]*m[2,2]-m[1,2]*m[2,1];
end;


Function Norme(m:matrice):MyReal;
begin
     Norme:= abs(m[1,1])
            +abs(m[1,2])
            +abs(m[2,1])
            +abs(m[2,2]);
end;


function acos(x:real):real;
begin
     if x=1 then
        acos:=0
     else if x=-1 then
     acos:=pi
     else
       if abs(x)<1 then
            acos:=atan(sqrt(1-sqr(x))/x)
       else
          begin
               Write('erreur dans acos');pause;
          end;
end;

procedure Gen_Mat(VAR x:matrice;x11,x12,x21,x22:Myreal);
begin
     x[1,1]:=x11;
     x[1,2]:=x12;
     x[2,1]:=x21;
     x[2,2]:=x22;
end;


procedure Inv_mat(m:matrice;var m1:matrice);
var
   d:real;
begin
     d:=determinant(m);
     Gen_mat(M1,m[2,2]/d,-m[1,2]/d,-m[2,1]/d,m[1,1]/d);
end;

procedure Alea_Mat(VAR m:matrice);
Var
   i,j:integer;
begin
     Randomize;
     for i:=1 to 2 do
      for j:=1 to 2 do
          m[i,j]:=2*random-1;
end;

procedure Aff_Mat(t:string;m:matrice);
Var
   i,j:integer;
begin
     WriteLN('[',t,']');
     for i:=1 to n do
     begin
      for j:=1 to n do
       begin
           write('  [',i,',',j,']=',m[i,j]:12:5);
       end;
      WriteLN;
      end;
      {WriteLN('Déterminant =', determinant(m):15:8);}
end;

procedure scal(var a:matrice;s:Myreal;b:matrice);
Var
   i,j : byte;
begin
     for i:=1 to n do
         for j:=1 to n do
             a[i,j]:=s*b[i,j];
end;

procedure add_mat(var a:matrice;a1,a2:matrice);
Var
   i,j : byte;
begin
     for i:=1 to n do
         for j:=1 to n do
             a[i,j]:=a1[i,j]+a2[i,j];
end;

procedure Sub_mat(var a:matrice;a1,a2:matrice);
begin
     Scal(a2,-1,a2);
     Add_mat(a,a1,a2);
end;

procedure mul_mat(var a:matrice;a1,a2:matrice);
Var
   i,j,k : byte;
   z: Myreal;
begin
     for i:=1 to n do
         for j:=1 to n do
           begin
             z:=0;
             for k:=1 to n do
                 z:=z+a1[i,k]*a2[k,j];
             a[i,j]:=z;
           end;
end;

procedure Calc_f(Var z:matrice;x,a,b,c:matrice);  {z=ax²+bx+c}
begin
     Mul_mat(z,a,x);
     {z:=ax+b}
     Add_mat(z,z,b);
     {z:=z*x}
     Mul_mat(z,z,x);
     {z:=z+c}
     Add_mat(z,z,c);
end;

procedure Resoud(Var x:matrice;a,b,c:matrice);
Var
   h,d,dmin     : Myreal;
   i1,i2,i3,i4  : integer;
   j1,j2,j3,j4  : integer;
   y,ymin,z:matrice;
begin
     h:=1;
     Gen_mat(x,0,0,0,0);
     repeat
      dmin:=1.0E10;
      for i1:=-1 to 1 do
         for i2:=-1 to 1 do
             for i3:=-1 to 1 do
                 for i4:=-1 to 1 do
                          begin
{                               WriteLN('i1,i2,i3,i4=',i1:4,i2:4,i3:4,i4:4);}
                               Gen_Mat(y,x[1,1]+i1*h,
                                         x[1,2]+i2*h,
                                         x[2,1]+i3*h,
                                         x[2,2]+i4*h);
                 {              Aff_mat('essai=',y);Pause;}
                               Calc_f(z,y,a,b,c);
                               d:=norme(z);
                               if d<dmin then
                                  begin
                                       j1:=i1;
                                       j2:=i2;
                                       j3:=i3;
                                       j4:=i4;
                                       ymin:=y;
                                       dmin:=d;
                                  end;
             end;
         x:=ymin;
         if (j1=0) and (j2=0) and (j3=0) and (j4=0) then
                 h:=h*0.8;
{         WriteLN('h=',h:8:5,' Norme=',dmin:8:5);}
          GotoXY(10,16);
          Aff_mat('x=',x);
          WriteLN('h=',h:15:12,' ³f(ax²+bx+c)³=',dmin:12:8);
      until dmin<0.00001;
      WriteLN;
      {Aff_mat(' ymin=',ymin);}
end;

procedure Iteration(Y:matrice);
const
   eps=1E-8;
var
   h:real;
   x,x1,U:matrice;
begin
     if Determinant(Y)<=0 then
        begin
             WriteLN('déterminant négatif ou nul');
        end
     else
      begin
        Gen_mat(x,1,0,0,1);
       repeat
               Efface;
               Aff_mat('Recherche de la racine de Y=',Y);
               Inv_mat(x,x1);
              Mul_mat(u,Y,x1);
              Add_mat(U,U,X);
              Scal(X,0.5,U);
              Aff_mat('X=',X);
              Mul_mat(U,X,X);
              Sub_mat(U,U,Y);
              h:=norme(U);
              writeLN('h=',h:12:8);
              Delay(100);
       until h<eps;
       end;
    Pause;
end;

procedure Iteration2(A,B,C:matrice);
const
   eps=1E-8;
var
   h:real;
   x,x1,U,V,I,A1,B1,C1:matrice;
begin
     Inv_mat(A,A1);
     Mul_mat(B1,A1,B);
     Mul_mat(C1,A1,C);
     Gen_mat(x,1,0,0,1);
       repeat
              Efface;
              Inv_mat(x,x1);
              Mul_mat(u,C1,x1);
              Add_mat(U,U,B1);
              Sub_mat(X,X,U);
              Scal(X,0.5,X);
              Aff_mat('X=',X);
              Delay(100);
       until false;
     Pause;
end;

procedure RacineCarree(Y:matrice);
var
   M,M1,S,R,X:matrice;
   det,disc:real;
   beta:real;
   rho,theta:real;
   Lambda1,Lambda2:real;
   Alpha:real;
   rho1:real;
begin
     Efface;
     Aff_mat('Y=',Y);
     det:=Determinant(y);
     WriteLN('Le déterminant de Y est ', det:12:8);
     if det<0 then
       begin
            WriteLN('Ce déterminant est négatif')
       end
     else if det=0 then
       begin
            WriteLN('Ce déterminant est nul')  ;
            WriteLN('Cherchons la valeur propre Lambda')  ;
       end
     else
       begin
            WriteLN('Ce déterminant est positif');
            WriteLN('Donc on va chercher les valeurs propres ');
            WriteLN('Elles ont solutions de l''équation caractéristique');
            WriteLN('z²-(y11+y22)z+d=0');
            beta:=y[1,1]+y[2,2];
            disc:=sqr(beta)-4*det;
            WriteLN('disc=',disc:12:8);
            if disc<=0 then
              begin
                   WriteLN('les solutions sont imaginaires de la forme');
                   WriteLN('rho(cos(theta)+isin(theta)');
                   WriteLN('rho(cos(theta)-isin(theta)');
                   rho:=sqrt(det);
                   WriteLN('beta=',beta:12:8, ' beta/2rho=',beta/(2*rho):12:8);
                   WriteLN('arccos=',acos(beta/(2*rho)):12:8);
                   theta:=acos(beta/(2*rho));
                   WriteLN('avec  rho=',rho:12:8,'  et theta=',theta:12:8);
                   WriteLN('Y=rho.MRM-1 où R est une matrice de rotation theta');
                   WriteLN('Donc la solution est de la forme');
                   WriteLN('MSM-1 où S est une similitude ');
                   rho1:=sqrt(rho);
                   WriteLN('de module  sqrt(rho)=',rho1:12:8);
                   alpha:=theta/2;
                   WriteLN('et d''angle alpha=',alpha:12:8);
                   if alpha=0 then
                   begin
                       WriteLN('alpha=0 ==> Il y a une infinité de solutions')
                   end
                   else
                   begin
                        Gen_mat(M,y[1,2]-rho*sin(theta),rho*cos(theta)-y[2,2],
                             rho*cos(theta)-y[1,1],rho*sin(theta)+y[2,1]);
                        Inv_mat(M,M1);
                         {  Mul_mat(X,M,M1); Aff_mat('vérification M*M1=',X);Pause;}
                        Gen_mat(S,rho1*cos(alpha),rho1*sin(alpha),
                            -rho1*sin(alpha),rho1*cos(alpha));
                        Mul_mat(X,S,M1);
                        Mul_mat(X,M,X);
                        Aff_mat('Rappel Y=',Y);
                        Aff_mat('Solution X= (au signe près)',X);
                        Mul_mat(X,X,X); Aff_mat('vérification X*X=',X);
                        Pause;
                   end
                 end
            else
              begin
                   WriteLN('Les deux valeurs propres sont ');
                   Lambda1:=(beta+sqrt(disc))/2;
                   Lambda2:=(beta-sqrt(disc))/2;
                   WriteLN('Lambda1=',Lambda1:12:8,' et  ','Lambda2=',Lambda2:12:8);
                   if (Lambda1>0)and(Lambda2>0) then
                      begin
                           WriteLN('les deux valeurs propres sont positives,');
                           WriteLN('donc les solutions sont de la forme :');
                           WriteLN('MSM-1 où S est diagonale ');
                           WriteLN('avec termes sqrt(Lambda1), sqrt(Lambda2)');
                           WriteLN('et M une matrice faite des deux vecteurs propres
associés');
                           Gen_mat(M,y[1,2],Lambda1-y[1,1],y[1,2],(Lambda2-y[1,1]));
                           Inv_mat(M,M1);
                        {   Mul_mat(X,M,M1); Aff_mat('verification M*M1=',X);Pause;}

                           Gen_mat(S,sqrt(Lambda1),0,0,sqrt(Lambda2));
                           Mul_mat(X,S,M1);
                           Mul_mat(X,M,X);
                           Aff_mat('Première solution X1= (au signe près)',X);

                          Mul_mat(X,X,X); Aff_mat('vérification X*X=',X);
                           Gen_mat(S,sqrt(Lambda1),0,0,-sqrt(Lambda2));
                           Mul_mat(X,S,M1);
                           Mul_mat(X,M,X);
                           Aff_mat('Deuxième solution X2= (au signe près)',X);
                           Mul_mat(X,X,X); Aff_mat('vérification X*X=',X);Pause;
                      end
              end;
       end;
     Pause;
end;

procedure Presentation;
begin;
      Efface;
      WriteLN(' ====== Equation du second degré de matrices =====');
      WriteLN('                                  ');
      WriteLN('                                                  ');
      WriteLN(' (1)  Solution de Xý=Y  par traitement direct     ');
      WriteLN(' (2)  Solution de Xý=Y  par itération X->(X+Y/X)/2');
      WriteLN(' (3)  Résolution par hypercubes                   ');
      WriteLN(' (4)  Résolution de   AXý+BX+C=0    par itération ');
      WriteLN('                              ');
      Write('   Tapez votre Choix  :  ');
end;

procedure Question1;
Var
   x,y,a,b,c:Matrice;
   d:real;
begin
     Efface;
     Writeln(':---------------------- Question n°1 -------------------------------:');
     Writeln(':                                                                   :');
     Writeln(':      Solution de X²=Y                                             :');
     Writeln(':-------------------------------------------------------------------:');
     Gen_mat(y,1,0,0,1);     RacineCarree(y);
     Gen_mat(y,4,3,3,4);     RacineCarree(y);
     Gen_mat(y,4,-3,3,4);    RacineCarree(y);
     Gen_mat(y,0,0,0,2);     RacineCarree(y);

end;

procedure Question2;
Var
   Y:Matrice;
begin
     Efface;
     Writeln(':---------------------- Question n°2 -------------------------------:');
     Writeln(':                                                                   :');
     Writeln(':      Solution de X²=Y  par itération X->(X+Y/X)/2                 :');
     Writeln(':-------------------------------------------------------------------:');
     Pause;
     Gen_mat(y,1,0,0,1);     Iteration(y);
     Gen_mat(y,4,3,3,4);     Iteration(y);
     Gen_mat(y,4,-3,3,4);    Iteration(y);
     Gen_mat(y,0,0,0,2);     Iteration(y);

     end;

procedure Question3;
Var
   x,a,b,c:Matrice;
begin
     Efface;
     Writeln(':---------------------- Question n°3 -------------------------------:');
     Writeln(':                                                                   :');
     Writeln(':      AX²+BX+C=0                                                   :');
     Writeln(':-------------------------------------------------------------------:');
     Pause;
     Gen_mat(a,2.8,0,0,2);
     Aff_mat('a=',a);
     Gen_mat(b,-1,-2,-2,-1);
     Aff_mat('b=',b);
     Gen_mat(c,0,1,1,0);
     Aff_mat('c=',c);
     Resoud(x,a,b,c);
     Aff_mat('Solution x=',x);
     Aff_mat('a=',a);
     Aff_mat('b=',b);
     Aff_mat('c=',c);
     Pause;
end;


procedure Question4;
Var
   x,a,b,c:Matrice;
begin
     Efface;
     Writeln(':---------------------- Question n°4 -------------------------------:');
     Writeln(':                                                                   :');
     Writeln(':      AX²+BX+C=0  par itération   X->  (A+I)X+C+C/X                :');
     Writeln(':-------------------------------------------------------------------:');
     Pause;
     Gen_mat(a,2.8,0,0,2);
     Aff_mat('a=',a);
     Gen_mat(b,-1,-2,-2,-1);
     Aff_mat('b=',b);
     Gen_mat(c,0,1,1,0);
     Aff_mat('c=',c);
     Iteration2(a,b,c);
     Aff_mat('Solution x=',x);
     Aff_mat('a=',a);
     Aff_mat('b=',b);
     Aff_mat('c=',c);
     Pause;
end;

begin
  Initgraphique;
  while true do
  begin
     Modetexte;
     Presentation;
     write('Question choisie ? ');
     readln(question);
     case question of
          '1':  Question1;
          '2':  Question2;
          '3' : Question3;
          '4' : Question4;
     end;
     Pause;
  end;
end.