function [stef,deb_acc]=steffensen(x0,exp_g,eps_seuil,n_max,nul_seuil)
% Produit accélérés de Steffensen d'une suite du point fixe et indice de
% début d'accélération.
% La convergence de la suite initiale est supposée connue,
% ainsi que celle des accélérés.
% Ce fichier bénéficie de la connaissance préalable de la fonction ait.
%
% variables d'entrée
% x0 est la valeur initiale,
% exp_g le nom de la fonction d'itération passé sous la forme 'g(x)'
% eps_seuil désigne comme annoncé dans l'ouvrage le seuil de variation
% relative des itérés,
% qui déclenche le processus d'accélération.
% n_max est la taille maximale du vecteur attendu en sortie;
% on se limite à des valeurs d'itérés distinctes au sens de nul_seuil.
% nul_seuil est un réel positif en-dessous duquel
% un réel est considéré nul.
%
% variables de sortie
% stef désigne le vecteur constitué des premiers itérés ordinaires de x0
% par la fonction g,
% suivis dès que possible au sens des paramètres passés,
% par les accélérés de Steffensen,
% avec arret dès que les suivants sont identiques au sens de nul_seuil
% choisi.
% deb_acc est l'indice au sens matlab du premier accéléré produit.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              Fonctions appelées                                                                    *
%                                                  
%                                           
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% spécifications
format long e;


% Tests d'entrée
if nargin~=5
    error('Nombre des arguments incorrect');
end
if n_max<2
    error('taille n_max insuffisante pour accélérer');
end
% Il convient sans doute de tester d'autres paramètres, leur type etc...
% selon la dangerosité des utilisateurs!
% Test complémentaire sur donnée
if n_max<2
    error('taille n_max insuffisante pour accélérer');
end

% Algorithme proprement dit

% Initialisation des vecteurs stef et util

stef=zeros(1,n_max);stef(1)=x0;
% util vecteur de longueur 4, contient initialement les itérés lents
% ordinaires.
x=x0;util(1)=x;
for i=1:3
    x=eval(exp_g);
    util(i+1)=x;
end

% Initialisations
% numéro du  premier candidat à l'accélération
% au sens des indices matlab
compt=2;deb_acc=2;

% Initialisations avant l'entrée en boucle
% en pratique toujours parcourue une fois...

q_util=diff(util(1:3))./diff(util(2:4));
rel=diff(q_util)/q_util(2);
% définition du booléen d'autorisation d'accélération: voir cours.
b_acc=(abs(rel)<eps_seuil);

% Puis on recherche les ratios relatifs à util,
% vecteur des seules valeurs utiles à l'accélération.
% Attention ici q_util(n) représente q_util d'indice n contrairement
% à l'ordinaire matlab.

% Attente d'accélération (comme en Aitken)
while (b_acc==0)&(compt<=n_max)
    
    % construction de stef avant accélération
    stef(compt)=util(2);
    
    % mise à jour de util
    x=util(4);
    util(1:3)=util(2:4);
    util(4)=eval(exp_g);
    
    % test pour accélération
    q_util=diff(util(1:3))./diff(util(2:4));
    rel=diff(q_util)/q_util(2);
    b_acc=(abs(rel)<eps_seuil);
    
    % mise à jour de compt et deb_acc
    compt=compt+1;deb_acc=deb_acc+1;
end


% Après une zone tampon, l'accélération systématique:
% on va changer de cheval!
% Le vecteur acc contiendra trois accélérés successifs;
% on calcule le premier sous ce bloc if.

% initialisations
acc=zeros(1,3);acc(1:2)=util(1:2);

% zone tampon
if compt<=n_max
    % premier style de calcul:utilisation de eval
    g_acc=[exp_g,'+(' exp_g '-acc(2))/((acc(2)-acc(1))/(' exp_g '-acc(2))-1)'];
    x=acc(2);val=eval(g_acc);
    stef(compt)=val;acc(3)=val;
    % préparation de la suite: glissement des accélérés
    % et mise à jour de compt
    acc(1:2)=acc(2:3);
    x=acc(2);acc(3)=eval(exp_g);
    compt=compt+1;
end

% Calcul des accélérés avec ré-entrance
while (compt<=n_max)&(abs(diff(acc(1:2)))>nul_seuil)
    
    % calcul de l'accéléré suivant
    % ici par formule du delta2 d'Aitken!
    val=acc(3)-(diff(acc(2:3))^2)/diff(diff(acc));
    stef(compt)=val;acc(3)=val;
    
    % Il convient de noter que les accélérés obtenus pourraient etre
    % améliorés en considérant les deux accélérés obtenus à partir de 
    % trois itérés successifs; ceci éviterait le phénomène qui apparait
    % sous la fonction entrelacs... appliquée aux steffensen.
    
    % préparation du tour suivant
    acc(1:2)=acc(2:3);
    x=acc(2);acc(3)=eval(exp_g);
    compt=compt+1;
end

% suppression des valeurs rendues inutiles après l'arret du calcul
stef;ind=find(stef~=0);stef=stef(ind);


% affichage des résultats à la demande de l'utilisateur
% Ce bloc peut etre remplacé par le passage d'un champ complémentaire
%dans la fonction steffensen par exemple dem_affichage,
% qui serait un booléen égal à 1
% ssi l'utilisateur veut afficher les résultats.
% Ainsi l'utilisateur par une copie de la fenetre de commande
% après exécution récupère
% et imprime l'appel de fonction avec les paramètres passés
% suivi des résultats obtenus.

ch=input('Voulez-vous afficher les résultats? (Si oui, taper 1)','s');
if str2num(ch)==1
    % à modifier selon l'impression souhaitée
    nb_lig=3;
    
    s=[];
    for i=1:size(stef,2)
        s=[s sprintf('stef%2.0f = %12.10f   ',i,stef(i))];
        % pourrait etre remplacé par un calibrage en fonction
        % de nul_seuil du type cal=-log10(2*nul_seuil)
        % à passer après transformation chaine à la place de 12.10f.
        if rem(i,nb_lig)==0|(i==size(stef,2))
            disp(s);sprintf('\n');
            s=[];
        end
    end
end
   
% fin de fonction steffensen

