%Composite Newton cotes rules demo
clc
clear all
close all
%composite rules
%cos function - composite simpson rule
%Range is 0 to 3
fname = 'humps';
a=0;
b=3;
m=3; % Simpson rule
c= WNC(m); %Quadrature weights
niter=10;  %number of cycles
% n is the number of points
%ni is number of intervals
%nseg = number of integration segments 
%using a 3-pt rule, two intervals are in each segment
%initialize xnew, xold
%
errold = 30.;
nold=2;
xold=[ a b];  %initial full integral range
fxold = feval(fname,xold);  %function at xold points
integExact =  2.39680798400577e+001;  %known answer
fprintf ('F(x) = %7s, a = %6.3f, b=%6.3f \n', fname, a, b)
fprintf ('Newton Cotes rule order m= %4.0f \n',m)
fprintf('Exact integral = %20.16f \n',integExact)
% on each interation, double the number of panels used
for in = 1:niter
    n = 2^in + 1;       %number of points in [a b]
    ni = n-1;           %number of intervals in range
    h=(b-a)/ni;         %size of each interval
    fx = zeros(1,n);
    x = linspace(a,b,n);  %nodes
    xnew = x(2:2:n-1);    %xnew is set of new points in this cycle
    %  the object function is evaluated only for the new points
    %  being introduced on this iteration
    % evaluation of function may be largest part of work
    fxnew = feval(fname,xnew);  
    %  fx is a merged set of the f(x) at old and new points
    fx(1:2:n) = fxold;  %resuses old points
    fx(2:2:n-1) = fxnew; %intersperses the new intermediates
    %fx is the entire set of function evaluations
    %set up vectors for entire set of regions
    index1=1;
    index2=1;
    nseg = ni/2;  %number of integration panels
                  % two intervals/panel in m=3 case.
                  % there are three weights per panel
                  % ccomp will hold them in a column
    ccomp=zeros(nseg*3,1);
                  % there are also three values of the function
                  % per panel -- fxcomp holds them in a row
    fxcomp=zeros(1,nseg*3);
    % Loop over the panels, assembling the two vectors
    for ns = 1:nseg
        fxcomp(index1:index1+2)=fx(index2:index2+2);  %three f(x)'s / panel
        ccomp(index1:index1+2)=c;  %weights are all the same
        index1=index1+3;
        index2=index2+2;   %incremented by two because end point f(x)s
                           % are shared between adjacent panels.
    end
    %calaculate the approximate integral with a dot product over the 
    % entire range, summing them for each panel.
    % width of a panel is 2h or 2(b-a)/ni
    % composite rule integral
    integ=2*h*(fxcomp*ccomp);
    %  error
    errnew = integExact - integ;
    disp (sprintf('N = %2.0f panels Int= %20.16f  Error = %10.5e',...
        nseg,integ,errnew))
     fxold = fx;
end