function fig6_17

%  Integrate  f(x)  over  [a, b]  using midpoint and simpson's rule
%  with compensated summation
%

a=0; b=1;
exact=(exp(3)-1)/3;

% a=0; b=3*pi;
% exact=3*pi;

% midpoint
fprintf('\n Error Using Midpint Rule \n')
fprintf('\n n  \t I_M \t\t I_M + CS \n')
nk=20;
n=120;
for k=1:nk
    intervM(k)=2*n;
    n=intervM(k);
    xd=linspace(a,b,n+1);
    h=xd(2)-xd(1);
    sumMM=f(xd(1)+0.5*h);
    err=0;
    sumM=f(xd(1)+0.5*h);
    for j=2:n
        z=f(xd(j)+0.5*h)+err;
        q=sumMM;
        sumMM=q+z;
        err=q-sumMM+z;
        sumM=sumM+f(xd(j)+0.5*h);
    end
    MM(k)=h*sumMM;
    M(k)=h*sumM;
    errM(k)=abs(exact-M(k));
    errMM(k)=abs(exact-MM(k));
    fprintf('%i \t %7.1se \t %7.1e \n',intervM(k),errM(k),errMM(k));
end
fprintf('\n')

%  simpson
fprintf('\n Error Using Simpsons Rule \n')
fprintf('\n n  \t I_S \t\t I_S + CS \n')
nk=14;
n=20;
% nk=10;
% n=10;
for k=1:nk
    intervS(k)=2*n;
    n=intervS(k);
    xd=linspace(a,b,n+1);
    h=xd(2)-xd(1);
    sumS=f(xd(1));
    err=0;
    sumSS=f(xd(1));
    for j=2:2:n
        ff2=f(xd(j+1));
        z=4*f(xd(j))+2*ff2+err;
        q=sumSS;
        sumSS=q+z;
        err=q-sumSS+z;
        sumS=sumS+4*f(xd(j))+2*ff2;
    end
    S(k)=h*(sumS-ff2)/3;
    SS(k)=h*(sumSS-ff2)/3;
    errS(k)=abs(exact-S(k));
    errSS(k)=abs(exact-SS(k));
    fprintf('%i \t %7.1e \t %7.1e \n',intervS(k),errS(k),errSS(k));
end
fprintf('\n')

% plot error using loglog plot
clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.06,'P',0.04)
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

loglog(intervM,errM,'-or','MarkerSize',8,'LineWidth',1.6)
hold on
loglog(intervM,errMM,'-*b','MarkerSize',8,'LineWidth',1.4)
loglog(intervS,errS,'--or','MarkerSize',8,'LineWidth',1.6)
loglog(intervS,errSS,'--*b','MarkerSize',8,'LineWidth',1.4)
grid on
xlabel('n-axis')
ylabel('Error')
axis([1 1e8 1e-16 1e-4]);
set(gca,'ytick',[1e-16 1e-12 1e-8 1e-4])
legend({' Original',' Comp Sum'},'Location','NorthEast','FontSize',14,'FontWeight','bold')

set(gca,'FontSize',16,'FontWeight','bold')

%exportgraphics(gcf,'/Users/mark/Desktop/comsum.eps')

function y=f(x)
y=exp(3*x);
% y=x*sin(x);

% function y=df(x)
% y=exp(x);















