function [p, q] = eoc(h,e)
% The function eoc returns an "average" experimental order of convergence for
% the sequences given by h and e.
%
% This is useful if you want to verify a theoretical relation like
%         e_i <= C * h_i^p
% by a numerical experiment. eoc gives you ONE estimate of the convergence
% order p for given sequences h and e.
%
% In contrast to the often used "sequential" eoc
%                   log(e_i)
%                ------------
%                log(e_{i+1})
%     p_i =  ----------------------
%                   log(h_i)
%                ------------
%                log(h_{i+1})
%
% it has a "average" character, since the least squares problem
%   min 1/2 || A*x - y ||^2_2
%
% where A and y are given by
% A = [ log(h_1) ,  1  ; 
%         ...    , ... ;
%       log(h_n) ,  1  ];
%
% y = [ log(e_1) ; 
%         ...    ;
%       log(e_n) ];
%
% is solved. x(1) is a estimate of the convergence order p and q = x(2) = log(C).
% The computed p corresponds to the average slope in a loglog-plot.
%
% INPUT: Two column vectors h and e of same length.
%
% EXAMPLE: p = eoc( [1, 1/2, 1/4, 1/8]' , [1, 1/4, 1/16, 1/64]' )
%          gives p=2 (quadratic order of convergence)
  

% Build matrix A and vector y
A = [ log(h) , ones(size(h)) ];
y = log(e);

% Solve the least squares problem
x = (A'*A)\(A'*y);

% The first component is an estimate for the convergence rate p
p = x(1);
% The second component is an estimate for log(C) 
q = x(2);
