+function [stef,deb_acc]=steffensen(x0,exp_g,eps_seuil,n_max,nul_seuil)\r
+% Produit accélérés de Steffensen d'une suite du point fixe et indice de\r
+% début d'accélération.\r
+% La convergence de la suite initiale est supposée connue,\r
+% ainsi que celle des accélérés.\r
+% Ce fichier bénéficie de la connaissance préalable de la fonction ait.\r
+%\r
+% variables d'entrée\r
+% x0 est la valeur initiale,\r
+% exp_g le nom de la fonction d'itération passé sous la forme 'g(x)'\r
+% eps_seuil désigne comme annoncé dans l'ouvrage le seuil de variation\r
+% relative des itérés,\r
+% qui déclenche le processus d'accélération.\r
+% n_max est la taille maximale du vecteur attendu en sortie;\r
+% on se limite à des valeurs d'itérés distinctes au sens de nul_seuil.\r
+% nul_seuil est un réel positif en-dessous duquel\r
+% un réel est considéré nul.\r
+%\r
+% variables de sortie\r
+% stef désigne le vecteur constitué des premiers itérés ordinaires de x0\r
+% par la fonction g,\r
+% suivis dès que possible au sens des paramètres passés,\r
+% par les accélérés de Steffensen,\r
+% avec arret dès que les suivants sont identiques au sens de nul_seuil\r
+% choisi.\r
+% deb_acc est l'indice au sens matlab du premier accéléré produit.\r
+\r
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
+% Fonctions appelées *\r
+% \r
+% \r
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
+\r
+\r
+% spécifications\r
+format long e;\r
+\r
+\r
+% Tests d'entrée\r
+if nargin~=5\r
+ error('Nombre des arguments incorrect');\r
+end\r
+if n_max<2\r
+ error('taille n_max insuffisante pour accélérer');\r
+end\r
+% Il convient sans doute de tester d'autres paramètres, leur type etc...\r
+% selon la dangerosité des utilisateurs!\r
+% Test complémentaire sur donnée\r
+if n_max<2\r
+ error('taille n_max insuffisante pour accélérer');\r
+end\r
+\r
+% Algorithme proprement dit\r
+\r
+% Initialisation des vecteurs stef et util\r
+\r
+stef=zeros(1,n_max);stef(1)=x0;\r
+% util vecteur de longueur 4, contient initialement les itérés lents\r
+% ordinaires.\r
+x=x0;util(1)=x;\r
+for i=1:3\r
+ x=eval(exp_g);\r
+ util(i+1)=x;\r
+end\r
+\r
+% Initialisations\r
+% numéro du premier candidat à l'accélération\r
+% au sens des indices matlab\r
+compt=2;deb_acc=2;\r
+\r
+% Initialisations avant l'entrée en boucle\r
+% en pratique toujours parcourue une fois...\r
+\r
+q_util=diff(util(1:3))./diff(util(2:4));\r
+rel=diff(q_util)/q_util(2);\r
+% définition du booléen d'autorisation d'accélération: voir cours.\r
+b_acc=(abs(rel)<eps_seuil);\r
+\r
+% Puis on recherche les ratios relatifs à util,\r
+% vecteur des seules valeurs utiles à l'accélération.\r
+% Attention ici q_util(n) représente q_util d'indice n contrairement\r
+% à l'ordinaire matlab.\r
+\r
+% Attente d'accélération (comme en Aitken)\r
+while (b_acc==0)&(compt<=n_max)\r
+ \r
+ % construction de stef avant accélération\r
+ stef(compt)=util(2);\r
+ \r
+ % mise à jour de util\r
+ x=util(4);\r
+ util(1:3)=util(2:4);\r
+ util(4)=eval(exp_g);\r
+ \r
+ % test pour accélération\r
+ q_util=diff(util(1:3))./diff(util(2:4));\r
+ rel=diff(q_util)/q_util(2);\r
+ b_acc=(abs(rel)<eps_seuil);\r
+ \r
+ % mise à jour de compt et deb_acc\r
+ compt=compt+1;deb_acc=deb_acc+1;\r
+end\r
+\r
+\r
+% Après une zone tampon, l'accélération systématique:\r
+% on va changer de cheval!\r
+% Le vecteur acc contiendra trois accélérés successifs;\r
+% on calcule le premier sous ce bloc if.\r
+\r
+% initialisations\r
+acc=zeros(1,3);acc(1:2)=util(1:2);\r
+\r
+% zone tampon\r
+if compt<=n_max\r
+ % premier style de calcul:utilisation de eval\r
+ g_acc=[exp_g,'+(' exp_g '-acc(2))/((acc(2)-acc(1))/(' exp_g '-acc(2))-1)'];\r
+ x=acc(2);val=eval(g_acc);\r
+ stef(compt)=val;acc(3)=val;\r
+ % préparation de la suite: glissement des accélérés\r
+ % et mise à jour de compt\r
+ acc(1:2)=acc(2:3);\r
+ x=acc(2);acc(3)=eval(exp_g);\r
+ compt=compt+1;\r
+end\r
+\r
+% Calcul des accélérés avec ré-entrance\r
+while (compt<=n_max)&(abs(diff(acc(1:2)))>nul_seuil)\r
+ \r
+ % calcul de l'accéléré suivant\r
+ % ici par formule du delta2 d'Aitken!\r
+ val=acc(3)-(diff(acc(2:3))^2)/diff(diff(acc));\r
+ stef(compt)=val;acc(3)=val;\r
+ \r
+ % Il convient de noter que les accélérés obtenus pourraient etre\r
+ % améliorés en considérant les deux accélérés obtenus à partir de \r
+ % trois itérés successifs; ceci éviterait le phénomène qui apparait\r
+ % sous la fonction entrelacs... appliquée aux steffensen.\r
+ \r
+ % préparation du tour suivant\r
+ acc(1:2)=acc(2:3);\r
+ x=acc(2);acc(3)=eval(exp_g);\r
+ compt=compt+1;\r
+end\r
+\r
+% suppression des valeurs rendues inutiles après l'arret du calcul\r
+stef;ind=find(stef~=0);stef=stef(ind);\r
+\r
+\r
+% affichage des résultats à la demande de l'utilisateur\r
+% Ce bloc peut etre remplacé par le passage d'un champ complémentaire\r
+%dans la fonction steffensen par exemple dem_affichage,\r
+% qui serait un booléen égal à 1\r
+% ssi l'utilisateur veut afficher les résultats.\r
+% Ainsi l'utilisateur par une copie de la fenetre de commande\r
+% après exécution récupère\r
+% et imprime l'appel de fonction avec les paramètres passés\r
+% suivi des résultats obtenus.\r
+\r
+ch=input('Voulez-vous afficher les résultats? (Si oui, taper 1)','s');\r
+if str2num(ch)==1\r
+ % à modifier selon l'impression souhaitée\r
+ nb_lig=3;\r
+ \r
+ s=[];\r
+ for i=1:size(stef,2)\r
+ s=[s sprintf('stef%2.0f = %12.10f ',i,stef(i))];\r
+ % pourrait etre remplacé par un calibrage en fonction\r
+ % de nul_seuil du type cal=-log10(2*nul_seuil)\r
+ % à passer après transformation chaine à la place de 12.10f.\r
+ if rem(i,nb_lig)==0|(i==size(stef,2))\r
+ disp(s);sprintf('\n');\r
+ s=[];\r
+ end\r
+ end\r
+end\r
+ \r
+% fin de fonction steffensen\r
+\r