%%
%% Script: ML_FEM_OrderOfAccuracy.m
%%
%% Compute order of accuracy using L2 error measure
%%

clear all, close all, clc

%
% Data obtained from running 
%  ML_FEM_1D_DiffusionVariableCoeff_MMS.m 
% with the following mesh sequence
%   no_elements = { 4, 8, 16, 32, 64, 128, 256, 512, 1024 }
%

data = [
5.0000e-01,  2.375693e+00;
2.5000e-01,  6.843966e-01;
1.2500e-01,  1.770845e-01; 
6.2500e-02,  4.460979e-02;
3.1250e-02,  1.117224e-02;
1.5625e-02,  2.794296e-03;
7.8125e-03,  6.986512e-04;
3.9062e-03,  1.746675e-04;
1.9531e-03,  4.366748e-05;
];

h = data(:,1);
E = data(:,2);

logh = log10(h);
logE = log10(E);

% plot log(h) vs log(E)
plot(logh,logE,'-ok', logh,logE,'-b'); 
xlabel('log10(h)'); ylabel('log10(E)'); grid on; 
legend('convergence',4); drawnow;

% perform linear regression and obtain the correlation coefficient
p = polyfit(logh,logE,1);
R = corrcoef(logh,logE);

fprintf(1,'Order of accuracy:       %1.4f \n', p(1) );
fprintf(1,'Correlation coefficient: %1.4f \n', R(1,2) );




