aitkd2(F,x,epsx,w,maxit,ipr)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   aitkd2.m        1992-09-07  %   Accelerated Aitken's delta-square 
%       (c)         M. Balda    %          iteration process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        modified on 2006-05-15
%
%   Function aitkd2 is designated for acceleration of vector iteration
%   processes. The original algorithm (see Irons, Tuck) has been
%   complemented by a function fw for dynamic weighing of relaxation
%   factor for improving numerical stability of the iteration process. 
% 
% Forms of calls:
% ~~~~~~~~~~~~~~~
%  x       = aitkd2(F,x);
%  [x,cnt] = aitkd2(F,x);
%  [x,cnt] = aitkd2(F,x,epsx,w,maxit,ipr);
%     F       string; a name of the M-function of the form  x = F(x)
%             that produces one iteration step over a vector x
%     x       vector of iterated values
%     epsx    tolerance on the first norm of vector of differences of x 
%             in two consecutive iterations for finishing the process, 
%             default eps=1e-5
%     w       weighing coefficient for function fw, default w = 1
%             < 0     fw = -w     constant weight for relaxation factor R
%             = 0     fw = 0      zero R = no acceleration
%             < 1     fw = k*w/((k-1)*w+1+eps)
%             = 1     fw = 1      full R
%     maxit   max number of iterations,  default maxit=99
%     ipr     = 0     no print of intermediate results (default)
%             = p>0   print every p-th iteration
%     cnt     >0  -> reached number of iterations to the solution
%             -maxit -> solution not reached in maxit iterations 
%
%  Bruce M. Irons, Robert C. Tuck:
%  A version of the Aitken accelerator for computer iteration
%  International Journal for Numerical Methods in Engineering
%  Volume 1, Issue 3, Date: July/September 1969, Pages: 275-277
%
% An example:
% ~~~~~~~~~~~
% See the demAitkd2 script

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,cnt] = aitkd2(F,x,epsx,w,maxit,ipr)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

R  = 1;
e1 = eps+1;
d1 = 1e99*ones(length(x),1);

if nargin<6, ipr   = 0;    end
if nargin<5, maxit = 99;   end
if nargin<4, w     = 1;    end
if nargin<3, epsx  = 1e-5; end

if ipr>0
    disp(' ')
    disp('     iter        R       x(i)  ->')
end

for k = 0:maxit
    y = feval(F,x);
    if ipr>0 && rem(k,ipr)==0, disp([k,R,y.']); end
    do = d1;
    d1 = x-y;
    d2 = do-d1;
    if w>0
       fw = k*w/((k-1)*w+e1);
    else
       fw = -w;
    end
    R = (R+(R-1)*(d2.'*d1)/(d2.'*d2))*fw;
    x = R*d1+y;
    if norm(d1)<epsx, k=-k; break, end
end
k = -k;

if nargout==2, cnt=k; end