tel/TPmatlab/equation_nonlineaire/TP4g/steffensen.m

   1 function [stef,deb_acc]=steffensen(x0,exp_g,eps_seuil,n_max,nul_seuil)\r
   2 % Produit accélérés de Steffensen d'une suite du point fixe et indice de\r
   3 % début d'accélération.\r
   4 % La convergence de la suite initiale est supposée connue,\r
   5 % ainsi que celle des accélérés.\r
   6 % Ce fichier bénéficie de la connaissance préalable de la fonction ait.\r
   7 %\r
   8 % variables d'entrée\r
   9 % x0 est la valeur initiale,\r
  10 % exp_g le nom de la fonction d'itération passé sous la forme 'g(x)'\r
  11 % eps_seuil désigne comme annoncé dans l'ouvrage le seuil de variation\r
  12 % relative des itérés,\r
  13 % qui déclenche le processus d'accélération.\r
  14 % n_max est la taille maximale du vecteur attendu en sortie;\r
  15 % on se limite à des valeurs d'itérés distinctes au sens de nul_seuil.\r
  16 % nul_seuil est un réel positif en-dessous duquel\r
  17 % un réel est considéré nul.\r
  18 %\r
  19 % variables de sortie\r
  20 % stef désigne le vecteur constitué des premiers itérés ordinaires de x0\r
  21 % par la fonction g,\r
  22 % suivis dès que possible au sens des paramètres passés,\r
  23 % par les accélérés de Steffensen,\r
  24 % avec arret dès que les suivants sont identiques au sens de nul_seuil\r
  25 % choisi.\r
  26 % deb_acc est l'indice au sens matlab du premier accéléré produit.\r
  27 \r
  28 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
  29 %              Fonctions appelées                                                                    *\r
  30 %                                                  \r
  31 %                                           \r
  32 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
  33 \r
  34 \r
  35 % spécifications\r
  36 format long e;\r
  37 \r
  38 \r
  39 % Tests d'entrée\r
  40 if nargin~=5\r
  41     error('Nombre des arguments incorrect');\r
  42 end\r
  43 if n_max<2\r
  44     error('taille n_max insuffisante pour accélérer');\r
  45 end\r
  46 % Il convient sans doute de tester d'autres paramètres, leur type etc...\r
  47 % selon la dangerosité des utilisateurs!\r
  48 % Test complémentaire sur donnée\r
  49 if n_max<2\r
  50     error('taille n_max insuffisante pour accélérer');\r
  51 end\r
  52 \r
  53 % Algorithme proprement dit\r
  54 \r
  55 % Initialisation des vecteurs stef et util\r
  56 \r
  57 stef=zeros(1,n_max);stef(1)=x0;\r
  58 % util vecteur de longueur 4, contient initialement les itérés lents\r
  59 % ordinaires.\r
  60 x=x0;util(1)=x;\r
  61 for i=1:3\r
  62     x=eval(exp_g);\r
  63     util(i+1)=x;\r
  64 end\r
  65 \r
  66 % Initialisations\r
  67 % numéro du  premier candidat à l'accélération\r
  68 % au sens des indices matlab\r
  69 compt=2;deb_acc=2;\r
  70 \r
  71 % Initialisations avant l'entrée en boucle\r
  72 % en pratique toujours parcourue une fois...\r
  73 \r
  74 q_util=diff(util(1:3))./diff(util(2:4));\r
  75 rel=diff(q_util)/q_util(2);\r
  76 % définition du booléen d'autorisation d'accélération: voir cours.\r
  77 b_acc=(abs(rel)<eps_seuil);\r
  78 \r
  79 % Puis on recherche les ratios relatifs à util,\r
  80 % vecteur des seules valeurs utiles à l'accélération.\r
  81 % Attention ici q_util(n) représente q_util d'indice n contrairement\r
  82 % à l'ordinaire matlab.\r
  83 \r
  84 % Attente d'accélération (comme en Aitken)\r
  85 while (b_acc==0)&(compt<=n_max)\r
  86     \r
  87     % construction de stef avant accélération\r
  88     stef(compt)=util(2);\r
  89     \r
  90     % mise à jour de util\r
  91     x=util(4);\r
  92     util(1:3)=util(2:4);\r
  93     util(4)=eval(exp_g);\r
  94     \r
  95     % test pour accélération\r
  96     q_util=diff(util(1:3))./diff(util(2:4));\r
  97     rel=diff(q_util)/q_util(2);\r
  98     b_acc=(abs(rel)<eps_seuil);\r
  99     \r
 100     % mise à jour de compt et deb_acc\r
 101     compt=compt+1;deb_acc=deb_acc+1;\r
 102 end\r
 103 \r
 104 \r
 105 % Après une zone tampon, l'accélération systématique:\r
 106 % on va changer de cheval!\r
 107 % Le vecteur acc contiendra trois accélérés successifs;\r
 108 % on calcule le premier sous ce bloc if.\r
 109 \r
 110 % initialisations\r
 111 acc=zeros(1,3);acc(1:2)=util(1:2);\r
 112 \r
 113 % zone tampon\r
 114 if compt<=n_max\r
 115     % premier style de calcul:utilisation de eval\r
 116     g_acc=[exp_g,'+(' exp_g '-acc(2))/((acc(2)-acc(1))/(' exp_g '-acc(2))-1)'];\r
 117     x=acc(2);val=eval(g_acc);\r
 118     stef(compt)=val;acc(3)=val;\r
 119     % préparation de la suite: glissement des accélérés\r
 120     % et mise à jour de compt\r
 121     acc(1:2)=acc(2:3);\r
 122     x=acc(2);acc(3)=eval(exp_g);\r
 123     compt=compt+1;\r
 124 end\r
 125 \r
 126 % Calcul des accélérés avec ré-entrance\r
 127 while (compt<=n_max)&(abs(diff(acc(1:2)))>nul_seuil)\r
 128     \r
 129     % calcul de l'accéléré suivant\r
 130     % ici par formule du delta2 d'Aitken!\r
 131     val=acc(3)-(diff(acc(2:3))^2)/diff(diff(acc));\r
 132     stef(compt)=val;acc(3)=val;\r
 133     \r
 134     % Il convient de noter que les accélérés obtenus pourraient etre\r
 135     % améliorés en considérant les deux accélérés obtenus à partir de \r
 136     % trois itérés successifs; ceci éviterait le phénomène qui apparait\r
 137     % sous la fonction entrelacs... appliquée aux steffensen.\r
 138     \r
 139     % préparation du tour suivant\r
 140     acc(1:2)=acc(2:3);\r
 141     x=acc(2);acc(3)=eval(exp_g);\r
 142     compt=compt+1;\r
 143 end\r
 144 \r
 145 % suppression des valeurs rendues inutiles après l'arret du calcul\r
 146 stef;ind=find(stef~=0);stef=stef(ind);\r
 147 \r
 148 \r
 149 % affichage des résultats à la demande de l'utilisateur\r
 150 % Ce bloc peut etre remplacé par le passage d'un champ complémentaire\r
 151 %dans la fonction steffensen par exemple dem_affichage,\r
 152 % qui serait un booléen égal à 1\r
 153 % ssi l'utilisateur veut afficher les résultats.\r
 154 % Ainsi l'utilisateur par une copie de la fenetre de commande\r
 155 % après exécution récupère\r
 156 % et imprime l'appel de fonction avec les paramètres passés\r
 157 % suivi des résultats obtenus.\r
 158 \r
 159 ch=input('Voulez-vous afficher les résultats? (Si oui, taper 1)','s');\r
 160 if str2num(ch)==1\r
 161     % à modifier selon l'impression souhaitée\r
 162     nb_lig=3;\r
 163     \r
 164     s=[];\r
 165     for i=1:size(stef,2)\r
 166         s=[s sprintf('stef%2.0f = %12.10f   ',i,stef(i))];\r
 167         % pourrait etre remplacé par un calibrage en fonction\r
 168         % de nul_seuil du type cal=-log10(2*nul_seuil)\r
 169         % à passer après transformation chaine à la place de 12.10f.\r
 170         if rem(i,nb_lig)==0|(i==size(stef,2))\r
 171             disp(s);sprintf('\n');\r
 172             s=[];\r
 173         end\r
 174     end\r
 175 end\r
 176    \r
 177 % fin de fonction steffensen\r
 178 \r