tel/TPmatlab/equation_nonlineaire/TP4g/aitken.m

   1 function [ait,deb_acc]=aitken(x0,exp_g,eps_seuil,n_max,nul_seuil)\r
   2 % produit  les accélérés d'Aitken d'une suite du point fixe et l'indice de début d'accélération.\r
   3 % La convergence de la suite initiale est supposée connue,\r
   4 % ainsi que celle des accélérés.\r
   5 %\r
   6 % variables d'entrée\r
   7 % x0 est la valeur initiale,\r
   8 % exp_g le nom de la fonction d'itération passé sous la forme 'g(x)'\r
   9 % eps_seuil désigne comme annoncé dans l'ouvrage le seuil de variation\r
  10 % relative des itérés qui déclenche le processus d'accélération.\r
  11 % n_max est la taille maximale du vecteur attendu en sortie;\r
  12 % on se limite à des valeurs d'itérés distinctes au sens de nul_seuil.\r
  13 % nul_seuil est un réel positif en-dessous duquel un réel est considéré\r
  14 % nul;\r
  15 % plus précisément pour un vecteur d'itérés on ne prend plus en compte\r
  16 % les valeurs distantes de la valeur finale de moins de nul_seuil.\r
  17 %\r
  18 % variables de sortie\r
  19 % ait désigne le vecteur constitué des premiers itérés ordinaires de x0\r
  20 % par la fonction g,\r
  21 % suivis dès que possible au sens des paramètres passés,\r
  22 % par les accélérés d'Aitken,\r
  23 % avec arret dès que les suivants sont identiques au dernier terme\r
  24 % au sens de nul_seuil choisi.\r
  25 % deb_acc est l'indice au sens matlab du premier accéléré produit.\r
  26 \r
  27 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
  28 %              Fonctions appelées                                                                    *\r
  29 %                                                  \r
  30 %                                           \r
  31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
  32 \r
  33 \r
  34 % spécifications\r
  35 format long e;\r
  36 \r
  37 \r
  38 % tests d'entrée\r
  39 if nargin~=5\r
  40     error('Nombre des arguments incorrect');\r
  41 end\r
  42 if n_max<2\r
  43     error('taille n_max insuffisante pour accélérer');\r
  44 end\r
  45 % Il convient sans doute de tester d'autres paramètres, leur type etc...\r
  46 % selon la dangerosité des utilisateurs!\r
  47 \r
  48 \r
  49 % Algorithme proprement dit\r
  50 \r
  51 % Génération des itérés simples par la fonction g;\r
  52 % on peut remplacer par une écriture matricielle ce bloc. Il sera moins\r
  53 % lisible.\r
  54 x=x0;val(1)=x;\r
  55 for i=1:n_max-1\r
  56     x=eval(exp_g);\r
  57     val(i+1)=x;\r
  58 end\r
  59 \r
  60 \r
  61 % Mise à l'abri de misères!\r
  62 % On évite des valeurs identiques à partir d'un certain rang.\r
  63 ind=find(abs(diff(val))>nul_seuil);val=val(ind);\r
  64 n_max=size(val,2);\r
  65 \r
  66 % recherche des ratios de variation relative des éléments de val;\r
  67 % voir cours.\r
  68 \r
  69 % Ici q(n) représentera ce qui dans le cours est q d'indice n\r
  70 % contrairement à l'ordinaire matlab.\r
  71 q=diff(val(1:n_max-1))./diff(val(2:n_max));\r
  72 \r
  73 % variation des qn\r
  74 rel=diff(q)./q(2:n_max-2);\r
  75 ind=find(abs(rel)<eps_seuil);\r
  76 sprintf('Début de l''accélération x d''indice:%3.0f ',ind(1));\r
  77 deb_acc=ind(1)+1;\r
  78 \r
  79 % détermination du vecteur ait\r
  80 if ind(1)<=n_max-2\r
  81     % calcul des accélérés sous aitken\r
  82     ait=zeros(1,n_max-1);\r
  83     ait(1:ind(1))=val(1:ind(1));\r
  84     \r
  85     % version matricielle\r
  86     ait(ind(1)+1:n_max-1)=val(ind(1)+2:n_max)+diff(val(ind(1)+1:n_max))./(q(ind(1):n_max-2)-1);\r
  87 end\r
  88 \r
  89 \r
  90 % on suprime les valeurs présentes dans ait égales au sens de nul_seuil\r
  91 % à la dernière calculée considérée comme la limite.\r
  92 ind=find(abs(ait-ait(size(ait,2)))>=nul_seuil);\r
  93 ait=ait(ind);\r
  94 \r
  95 \r
  96 % affichage des résultats à la demande de l'utilisateur\r
  97 % Ce bloc peut etre remplacé par le passage d'un champ complémentaire\r
  98 % dans la fonction aitken, par exemple dem_affichage,\r
  99 % qui serait un booléen égal à 1\r
 100 % ssi l'utlisateur veut afficher les résultats.\r
 101 % Ainsi l'utilisateur par une copie de la fenetre de commande\r
 102 % après exécution récupère et imprime l'appel de fonction\r
 103 % avec les paramètres passés, suivi des résultats obtenus.\r
 104 \r
 105 \r
 106 ch=input('Voulez-vous afficher les résultats? (Si oui, taper 1)','s');\r
 107 if str2num(ch)==1\r
 108     % à modifier selon l'impression souhaitée\r
 109     nb_lig=3;\r
 110     \r
 111     s=[];\r
 112     for i=1:size(ait,2)\r
 113         s=[s sprintf('ait%2.0f = %12.10f   ',i,ait(i))];\r
 114         % pourrait etre remplacé par un calibrage en fonction\r
 115         % de nul_seuil du type cal=-log10(2*nul_seuil)\r
 116         % à passer après transformation chaine à la place de 12.10f.\r
 117         if rem(i,nb_lig)==0|(i==size(ait,2))\r
 118             disp(s);sprintf('\n');\r
 119             s=[];\r
 120         end\r
 121     end\r
 122 end\r
 123 \r
 124 % fin de fonction aitken\r