function fig6_12

%  Integrate  f(x)  over  [a, b]  using simp-rom, and gauss

%%% pick your example
% exp(3x)
a=0; b=1;
exact=(exp(3)-1)/3;

% tanh(100*x-30)
% a=0; b=1;
% r1=1+exp(-140); r2=1+exp(-60);
% exact=0.4+log(r1/r2);

% 1/(25*x^2 + 1)
% a=-1; b=1;
% exact=(2*atan(5))/5;

%  composite simpson-romberg
n=2;
i=1;
S0=simp(a,b,n);
fnS(1)=3;
R0=S0;
errS(1)=abs(R0-exact)/exact;
while n<16
    n=2*n;
    i=i+1;
    fnS(i)=n+1;
    S=simp(a,b,n);
    R=(16*S-S0)/15;
    errS(i)=abs(R-exact)/exact;
    S0=S;
    R0=R;
end

% gauss
for k=1:16
    G(k)=gauss(a,b,k);
    fnG(k)=k;
    errG(k)=abs(exact-G(k))/exact;
end

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)

semilogy(fnS,errS,'--*r','MarkerSize',8,'lineWidth',1.6)
hold on
box on
grid on
axis([3 17 1e-16 1])
set(gca,'xtick',[3 5 7 9 11 13 15 17 19 23])
set(gca,'ytick',[1e-16 1e-12 1e-8 1e-4 1])
xlabel('Number of Points Used')
ylabel('Error')
semilogy(fnG,errG,'--ob','MarkerSize',8,'lineWidth',1.6)

legend({' Simp-Romb',' Gaussian'},'Location','NorthEast','FontSize',16,'FontWeight','bold')

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

%exportgraphics(gca,'/Users/mark/Desktop/gauss.eps')


function y=f(x)
y=exp(3*x);
% y=tanh(100*x-30);
%  y=1/(25*x^2 + 1);

%%% gaussian quadrature: n = number of points
function s=gauss(a,b,n)
beta=0.5./sqrt(1-(2*(1:n-1)).^(-2));
A=diag(beta,1)+diag(beta,-1);
[V,D]=eig(A);
x=diag(D);
x=a+0.5*(b-a)*(1+x);
w=2*V(1,:).^2;
sum=0;
for i=1:length(x)
    sum=sum+w(i)*f(x(i));
end
s=sum*(b-a)/2;

%%% composite simpson: n = number of subintervals
function s=simp(a,b,n)
xd=linspace(a,b,n+1);
h=xd(2)-xd(1);
s=f(xd(1));
for j=2:2:n
    ff2=f(xd(j+1));
    s=s+4*f(xd(j))+2*ff2;
end
s=h*(s-ff2)/3;
































