PROGRAM Recurrences;

{
Epreuve informatique de l'Ecole Polytechnique                   93.161


1°  Dans un espace euclidien de dimension n désigné par E,
associe à tout couple de points A et B un vecteur noté AB et une
distance notée AB. Le produit scalaire de deux vecteurs u et v se
note u.v .

Ecrire un programme Pascal pour étudier la
convergence éventuelle de la suite de points Mn définis par
récurrence par l'équation Mk+3= f(Mk, Mk+1, Mk+2 ) où les trois
premiers points M0, M1,M2 sont fixés. Discuter la convergence en
fonction des paramètres Ó, ß, þ dans le cas où :



 ------------------------------------------------------------+
¦              Imprimer tous les résultats                   ¦
¦     en indiquant chaque fois à quoi ils correspondent      ¦
+------------------------------------------------------------+

                            -=-=-=
 }



uses
  crt, graph, modubase,matric;


Const
     maxd=3;
     maxk=9;

TYPE
    Angle = real;
    Points = record
              x : array[1..maxd] of real;
              t  : string[1]
            end;

    Vecteur = record
              x : array[1..maxd] of real;
              end;

    Matrice = array[1..maxk,1..maxk] of real;



Const
     maxn=50;
     reduction = 1;
VAR
     Car         : Char ;
     Question    : Char;
     Xmin, Xmax,Ymin,Ymax, yv : real;
     Cas         : integer;
     i,k,n,d       : integer;
     W           : Points;
     a : array[1..maxk] of real;
     O :   vecteur;

Type
    Figure = array [1..maxk+1] of points;


Procedure Alea_P(VAR p:points);
var
   i:integer;
begin
     for i:=1 to d do
         p.x[i] := random*2-1;
end;


Procedure Vectorise(VAR v:Vecteur;P1,P2:Points);
begin
     for i:=1 to d do
         v.x[i]:=P2.x[i]-P1.x[i];
end;

Procedure Translate(v:Vecteur;P1:Points;Var P2:Points);
var
   i:integer;
begin
     for i:=1 to d do
         P2.x[i]:=P1.x[i]+v.x[i];
end;

Procedure Homothetie(lambda:real;v1:Vecteur;Var v2:Vecteur);
var
   i:integer;
begin
     for i:=1 to d do
         v2.x[i]:=v1.x[i]*lambda;
end;

Function Scalaire(u,v:Vecteur):Real;
var
   i:integer;
   x:real;
begin
     x:=0;
     for i:=1 to d do
         x:=x+u.x[i]*v.x[i];
     Scalaire := x;

end;

Function Module2(V:Vecteur):Real;
var
   i:integer;
   x:real;
begin
     Module2 := scalaire(v,v);
end;

Function Module(v:Vecteur):Real;
begin
     Module := SQRT(Module2(v));
end;

Function Distance(A,B:Points):Real;
VAR
   v:vecteur;
BEGIN
     Vectorise(V,A,B);
     Distance := Module(v);

END;


PROCEDURE Normalise(VAR v:vecteur);
VAR
   M : Real;
begin
     M := Module(v);
     for i:=1 to d do
         v.x[i] := v.x[i]/M;
end;

Function ProjX(p:points):real;
begin
     ProjX := (p.x[2]-p.x[1]*0.15)*reduction;
end;

Function ProjY(p:points):real;
begin
     ProjY := (p.x[3]-p.x[1]*0.32)*reduction;
end;

Procedure TracePoint(p:points);
VAR
   x,y : real;
begin
     x := ProjX(p);
     y := ProjY(p);
     croix(x,y);
     Ecris(p.t);
end ;


Procedure TraceTrait(p1,p2:Points);
BEGIN
      Deplace(ProjX(p1),ProjY(p1));
      Trace(ProjX(p2),ProjY(p2));
END;


Procedure TraceAxe(x,y,z:real;t:string);
VAR
   p : points;
   i:integer;
begin
     p.x[1]:=x;
     p.x[2]:=y;
     p.x[3]:=z;
     p.t:=t;
     TraceTrait(W,p);
     Croix(ProjX(p),ProjY(p));
     Ecris(t);
end;




Procedure Add_vec(V,W:Vecteur;VAR U:Vecteur);
Var
   i:integer;
begin
     for i:=1 to d do
         u.x[i] := V.x[i]+W.x[i];
end;



Procedure Mat_Vec(m:matrice;u:vecteur;Var v:vecteur);
Var
   i,j:integer;
   z:real;
begin
     for i:=1 to d do
         begin
              z:=0;
              for j:=1 to d do
                  z:=z+m[i,j]*u.x[j];
              v.x[i]:=z;
         end;
end;

Function Discussion:string;  {  cas  A1+a2+a3=1}

         {      Racines = ½(a1-1±¹(sqr(1-a1)-4a3))
                                                                 }
Va
   x,det,re,im:real;
   d,s:string;
   j:integer;
begin
     det:=sqr(1-a[1])/4-a[3];
     re:=(a[1]-1)/2;
     if det<0 then
        begin
             Im:=sqrt(-det);
             Str(re:12:8,s);
             d:='Racines complexes : Re+i*Im='+s;
             Str(Im:12:8,s);
             d:=d+'  ±i'+s;
        end
     else
         begin
              im :=sqrt(det);
              Str(re-im:12:8,s);
              d:='Racines réelles='; d:=d+' '+s;
              Str(re+im:12:8,s);
              d:=d+' & '+s;
         end;
    if det>0 then
       begin
            if abs(re-im)>abs(re+im) then
               x:=re-im
            else
               x:=re+im;
            Str(x:12:8,s);
            if abs(x)>1 then
               d:=d+' Div'
                  else
               d:=d+' Conv';
            if x>0 then
               d:=d+' stable '+s
            else
               d:=d+' altern '+s;
       end
    else
       begin
            if (abs(Re)>1) or (abs(Im)>1) then
               d:=d+' Puls div'
            else
               d:=d+' Puls conv';
            if abs(Re)>Im then
               begin
                    x:=Re;
                    d:=d+ ' Re='
               end
            else
                begin
                    if Re>0 then
                       x:=Im
                    else
                        x:=-Im;
                    d:=d+ ' Im='
                end;
            Str(x:12:8,s);
            d:=d+s;
       end
    Discussion:=d;

end;

Procedure Recurrence(M:figure;Var p:Points;k:integer);
Var
   i:integer;
   u:vecteur;
begin
     p:=W;
     for i:=1 to k do
         begin
              Vectorise(u,W,M[k+1-i]);
              homothetie(a[i],u,u)
              Translate(u,p,p);
         end;
end;

procedure init(titre:string);
begin
     Modegraphique;
     Efface;
     Xmin := -3;
     Xmax := +4;
     Ymin := -2;
     Ymax := 2;
     Fenetre(Xmin,Xmax,Ymin,Ymax);
     yv := 1.2;
     Deplacf(Xmin,Ymin*0.9);
     Ecris(Titre);

     SetBkColor(Blue);

     Couleur(Rouge);
     for i:=1 to d do
         W.x[i]:=0;
     W.t:='W';
     TracePoint(W);
     TraceAxe(1.2,0,0,'x');
     TraceAxe(0,1.2,0,'y');
     TraceAxe(0,0,1.2,'z');

     Couleur(-Brillant);
end;                     (* init *)





Procedure Question1;
Var
   i:integer;
   k:integer;
   M: Figure;
   x,z,zz:real;
   s:string;
begin
     ModeTexte;
     Efface;
     Writeln('+---------------------- Question n°1 -------------------------------+');
     Writeln('¦                                                                   ¦');
     Writeln('¦     (1)  Convergence dans Espace affine euclidien R3              ¦');
     Writeln('¦                                                                   ¦');
     Writeln('¦      Mn=a1*Mn-1+a2*Mn-2+...+ak*Mn-k    M1,M2...Mk fixés           ¦');
     Writeln('¦                                                                   ¦');
     Writeln('¦                                                                   ¦');
     Writeln('+-------------------------------------------------------------------+');
{     Write('Entrez la valeur de k=');
     ReadLN(k); }
     k:=3;
     Randomize;
     repeat
          Init('Convergence Mn=a1*Mn-1+a2*Mn-2+...+ak*Mn-k');
          begin
               repeat
                     for i:=1 to k do
                         begin
                              a[i]:=2*random-1;
                              Alea_P(M[i]);
                         end;
                    a[3]:=1;
                    a[2]:=2*Random-1;
                    a[1]:=1-a[2]-a[3];
                     x:=0;
                     for i:=1 to k do
                              x:=x+a[i];
               until x<>0;
               for i:=1 to k do
                   a[i]:=a[i]/x;

           Couleur(-brillant);
           for i:=1 to k-1 do
             TraceTrait(M[i],M[i+1]);
           n:=k;
           repeat
               Deplace(xmin,ymax*0.9);
               Ecris('a[1..k]=');
               for i:=1 to k do
                   Ecrisreel(a[i]);
                 Inc(n);
                 Deplace(xmin,ymax*0.85);
                 Ecris('n=');Ecrisentier(n);
                 Recurrence(M,M[k+1],k);
                 TraceTrait(M[k+1],M[k]);
                 Deplace(xmin,ymax*0.8);
                 Ecris('x[1..d]=');
                 for i:=1 to d do
                   Ecrisreel(M[k+1].x[i]);
                 Deplace(xmin,ymax*0.75);
                 Ecris('MnMn-1=');
                 Ecrisreel(distance(M[k],M[k+1]));
                 Deplace(xmin,ymax*0.70);
                 Ecris('  %=');
                 for i:=1 to k-1 do
                     Ecrisreel(distance(M[k+2-i],M[k+1-i])/distance(M[k+1-i],M[k-i]));
                 Str(-Exp(Ln(distance(M[k+1],M[k]))/n):12:8,s);
                 Ecris(' ==>        '+s);
                 Deplace(xmin,ymax*0.65);
                 Ecris(Discussion);
                 for i:=1 to k do
                   M[i]:=M[i+1];
                 Delay(100);
           until  KeyPressed or (Abs(Log(distance(M[k-1],M[k])))>12);
           Car:=ReadKey;
          end;
    until Car=#27;

end;

Procedure Question2;
Var
   i,k:integer;
   M:array[0..3] of Points;
   x,y,z:real;
begin
     Efface;
     Writeln('+---------------------- Question n°2 -------------------------------+');
     Writeln('¦                                                                   ¦');
     Writeln('¦   Calcul itératif de Mk+3=f(Mk+2,Mk+1,Mk)   dans R2               ¦');
     Writeln('¦                      D=f(A,B,C)                                   ¦');
     Writeln('¦                                                                   ¦');
     Writeln('¦   AD=phi(r,s,t)AB+psi(r,s,t)AC                                    ¦');
     Writeln('¦                                                                   ¦');
     Writeln('¦       où  r=AC/AB s=signe(AB/\AC) et t=cos(BAC)                   ¦');
     Writeln('¦                                                                   ¦');
     Writeln('¦   (1) D symétrique de C par rapport à la droite AB                ¦');
     Writeln('¦   (2) D est le centre du cercle inscrit dans le triangle ABC      ¦');
     Writeln('+-------------------------------------------------------------------+');
     repeat
           write('Cas choisi ? ');
           readln(cas);
     until cas in [1..2];

     Randomize;
     repeat
     until car=#27;
end;


Procedure Question3;
Var
   i,k:integer;
   M:array[0..3] of Points;
   x,y,z:real;
begin
     Efface;
     Writeln('+---------------------- Question n°3 -------------------------------+');
     Writeln('¦                                                                   ¦');
     Writeln('¦   Calcul itératif de Mk+3=f(Mk+2,Mk+1,Mk)   dans R3               ¦');
     Writeln('¦                      D=f(A,B,C)                                   ¦');
     Writeln('¦                                                                   ¦');
     Writeln('¦   (1) Barycentre                                                  ¦');
     Writeln('¦   (2) AD=(Ó/BC)(AB/\AC)+ßAB+þAC                                   ¦');
     Writeln('+-------------------------------------------------------------------+');
     repeat
           write('Cas choisi ? ');
           readln(cas);
     until cas in [1..2];

     Randomize;
end;





Procedure Presentation;
begin
    ModeTexte;
    Efface;
    WriteLN('      Suite récurrente dans Rd         ');
    WriteLN('                                       ');
    WriteLN('                                       ');
    WriteLN(' (1)  Convergence dans R3              ');
    WriteLN(' (2)  ?????????????????????????????    ');
    WriteLN(' (3)  ??????                           ');
    WriteLN('                                       ');
end;



begin
  d:=3;
  Initgraphique;
  while true do
  begin
     Question1;
     Presentation;
     write('Question choisie ? ');
     readln(question);
     Efface;
     writeln('======== Question n°',question,'===============');
     writeln('');

     case question of
          '1':  Question1;
          '2':  Question2;
          '3' : Question3;
     end;
  end;

end.