function table6_6

%  Integrate  f(x)  over  [a, b]  using Gaussian quadrature (no subintervals) and composite simpson
%  Compares computed and exact value
%  f(x) given at end of file

global x

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

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

% composite simpson
ic=0;
for k=2:10
    % n =  number of subintervals
    n=2^k;
    if k==10
        n=1000;
    end
    ic=ic+1;
    
    % simpson
    evalsS(ic)=n+1;
    I_S(ic)=simp(a,b,n);
    errS(ic)=abs(exact-I_S(ic));
    
end

% gauss
ic=0;
for k=2:80
    ic=ic+1;
    G(ic)=gauss(a,b,k);
    evalsG(ic)=k;
    errG(ic)=abs(exact-G(ic));
    
    if k<16
        figure(1)
        clf
        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)
        % get(gcf)
        set(gcf,'Position', [6 1062 627 283])
        plot(x,zeros(length(x)),'xr','MarkerSize',12,'LineWidth',1.5)
        say=['m = ', int2str(k)];
        text(0.5*(a+b),0.4,say,'FontSize',16,'FontWeight','bold')
        axis([a b -1 1])
        xlabel('Location of Gauss Nodes')
        set(gca,'FontSize',16,'FontWeight','bold')
        pause
    end
end

% format shortE
% errG

figure(2)
clf
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)
% get(gcf)
set(gcf,'Position', [6 1062 627 283])
loglog(evalsS,errS,'--ob','MarkerSize',8,'LineWidth',1.3)
hold on
loglog(evalsG,errG,'-.+r','MarkerSize',9,'LineWidth',1.3)
grid on
xlabel('Number of f(x) Evaluations')
ylabel('Error')
axis([1 1.01e3 1e-16 1])
set(gca,'ytick',[1e-16 1e-11 1e-6 1e-1])
set(gca,'YMinorGrid','off')
legend({' Simpson',' Gauss'},'Location','SouthWest','FontSize',16,'FontWeight','bold')
set(gca,'FontSize',16,'FontWeight','bold')



%%% gaussian quadrature: n = number of points
function s=gauss(a,b,n)
global x
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;



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