quadg.m

Numerical quadrature using open methods.

   %  End of the "real work".  The code below is
   %  for the graphical and optional output.
   
   %  Check verbosity and output storage
   if OutputOn
      %  Store the Nodes and Integrand Values
      AllNodes =  [AllNodes; Nodes(:)];
      AllIntegrand =  [AllIntegrand; Integrand(:)];
      if verbosity
          % Update the listbox
          updatefigure(Listbox,Axis,Iter,Nodes,GoodRegions,Integrand)          
      end   
   end  %   End of verbose output section

end 
%  End of the main while loop
%--------------------------------------------------------------

if Iter>=MaxIter
   warning('Reached maximum iterations; result may be inaccurate')
end 

%______________________________________________________________
%  Define output structure and update display if necessary
%--------------------------------------------------------------
if OutputOn
   ElapsedTime              = cputime - t0;
   Output.Iterations        = Iter;
   Output.FunctionSamples   = length(AllNodes);
   [Output.Nodes idx]       = sort(AllNodes);
   Output.Integrand         = AllIntegrand(idx);
   Output.Time              = ElapsedTime;
   Output.basicNodes        = Options.nodes;
   Output.basicWeights      = Options.weights;
   if verbosity
      axes(Axis)
         %  Plot the points of the integrand used to compute the result
         AllNodes = Output.Nodes; AllIntegrand = Output.Integrand;
         plot(AllNodes,AllIntegrand,'Linewidth',1.5,'Color','r')

         %  Update the title
         Title = get(Axis,'title');
         title([get(Title,'String') ' = ' num2str(Result)])
         set(Title,'Fontweight','bold')

         %  Plot a histogram of the samples
      axes(HistogramAxis)
         if round(length(AllNodes)/15)>1
            bins = min(round(length(AllNodes)/15),length(AllNodes));
         else
            bins = length(AllNodes);
         end
         hist(AllNodes,linspace(min(AllNodes),max(AllNodes),bins))
         xlim([min(AllNodes) max(AllNodes)])
         ylabel('Count')
         Tstring = [num2str(length(AllIntegrand)) ' fcn samples, ',...
                    num2str(ElapsedTime) ' seconds.'];
         title(Tstring,'Fontweight','Bold','Fontname','FixedWidth');

      setptr(Fig,'arrow')
      drawnow
   end
end
%--------------------------------------------------------------


%*********************************************************
%  End of main function, subfunctions are below
%  First the computationally relevant subfunctions, 
%  followed by the graphical and interface related 
%  functions.
%*********************************************************


%*********************************************************
function ScaledX = scalex(X,a,b);
%  This function scales the abscissae to the proper values, 
%  in a vectorized fashion.
%
%   I used KRON, but this should probably be rewritten to 
%   avoid this.
a = a(:); b = b(:);
Node_Scale  = (b-a)/2;
ScaledX = kron(Node_Scale.',X) + repmat((Node_Scale+a).',length(X),1);


%*********************************************************
function [X, W] = gausslegendre(order)
%   This function generates the nodes and weights
%   for a Gauss Legendre quadrature over [-1,1].
%   Note that the weights and nodes for both the lower
%   order and higher order accuracy check are created, so:
%   
%   X(:,1:?) = nodes for lower order quadrature
%   X(:,?:end) = nodes for higher order quadrature
%       and
%   W(1,:) = weights for lower order quadrature
%   W(2,:) = weights for higher order quadrature
%
%       and
%
%   Q = W*feval(X.')*span_of_integration
%   where Q(1,:) = quadrature estimate from lower order method
%   and   Q(2,:) = quadrature estimate from higher order method.

%  Generate the weights and the abscissa for the first pass
vector =(1:order-1)./sqrt((2*(1:order-1)).^2-1);
[w1,xi1]=eig(diag(vector,-1)+diag(vector,1));
xi1 = diag(xi1);
w1=w1(1,:)'.^2;

%  Now generate the weights and abscissae for the second
%  pass.  Append that to the first pass values, so that
%  W*X gives [first_pass_quadrature; second_pass_quadrature];
vector =(1:order)./sqrt((2*(1:order)).^2-1);
[w2,xi2]=eig(diag(vector,-1)+diag(vector,1));
xi2 = diag(xi2);
w2=w2(1,:)'.^2;

X = [xi1; xi2];

First_w  = [w1(:).' zeros(1,order+1)];
Second_w = [zeros(1,order) w2(:).'];
W = [First_w; Second_w];

%  Normalize the weights - depending on how you scale various values,
%  this may not be necessary.
%  I'm doing this just to be careful about things...
W = diag(1./sum(W,2))*W;
%*********************************************************


%*********************************************************
function [X, W] = gausskronrod(order)
%   This function generates the nodes and weights
%   for a Gauss Kronrod quadrature over [-1,1].
%
%   See gausslegendre for an explanation of how the
%   outputs are specified.

if order~=7
    error('Only 7(15) order available for Gauss-Kronrod')
    return
end

%  First get the Gauss-Legendre weights
[GX GW] = gausslegendre(order);

%  Get the weights and abscissae I need and
%  sort them to make later work easier
GW = GW(1,1:order); GX = GX(1:order,1);
[GX idx] = sort(GX);
GW = GW(idx);
%  Now add on the new abscissae and form the 
%  new weights to make up the Kronrod extension.
Positive_Kronrod_Additional_Nodes = [...
        0.9914553711208126392068546; 0.8648644233597690727897127;
        0.5860872354676911302941448; 0.2077849550078984676006894];
Kronrod_Weights = [...
        0.0229353220105292249637320, 0.0630920926299785532907006,...
        0.1047900103222501838398763, 0.1406532597155259187451895,...
        0.1690047266392679028265834, 0.1903505780647854099132564,...
        0.2044329400752988924141619, 0.2094821410847278280129991];
Kronrod_Weights_for_Original_Nodes = Kronrod_Weights(2:2:end);
Kronrod_Weights_for_Additional_Nodes = Kronrod_Weights(1:2:end);

%  Put all the nodes together and sort them
X = [   GX(:);...
        Positive_Kronrod_Additional_Nodes;...
        -Positive_Kronrod_Additional_Nodes];
Kronrod_Weights = Kronrod_Weights(idx);

W = [   GW(:).' zeros(1,order+1);
        Kronrod_Weights_for_Original_Nodes,...
        Kronrod_Weights_for_Original_Nodes((end-1):-1:1),...
        Kronrod_Weights_for_Additional_Nodes,...
        Kronrod_Weights_for_Additional_Nodes];

%  Normalize the weights
W = diag(1./sum(W,2))*W;
%*********************************************************






%*********************************************************
%  Graphics and interface subroutines below.
%*********************************************************
function updatefigure(Listbox,Axis,Iter,Nodes,GoodRegions,Integrand)
%  Update the listbox
Str = sprintf('Iteration %2.1i %4.1i (%4.3g%%) %8.1i %9.1i',...
    Iter,...
    full(sum(GoodRegions)),...
    full(100*sum(GoodRegions)/length(GoodRegions)),...
    full(sum(~GoodRegions)),...
    prod(size(Integrand)));

Current_String = get(Listbox,'String');
set(Listbox,'String',[Current_String(1);{Str};Current_String(2:end)])

%  Update the plot
axes(Axis)
stem(Nodes(:),Integrand(:));
drawnow

%  Pause to so that the animation can be seen
pause(.25)
%*********************************************************



%*********************************************************
function [Fig, Axis, Listbox, HistogramAxis, Spaces] = ...
         initializefigure(Options,fun,a,b);
   %  Initialize figure window display
   if isa(fun,'function_handle');
       fun = func2str(fun);
   else
       fun = char(fun);
   end
   Fig = figure('Name',[' Integration of ' fun],...
                       'Doublebuffer','on',...
                       'NumberTitle','off',...
                       'Units','normalized','Position',[.15 .2 .6 .7]);
         setptr(Fig,'watch')
   Axis = subplot(3,1,1);
      set(Axis,'units','normalized','xlim',[min(a,b) max(a,b)])
      AxPos = get(Axis,'position');
   TexTitle = ['_{' num2str(a) '}^{' num2str(b) '}' texlabel(fun) ];

   Title = title(['\int' TexTitle]);
      ylabel('Integrand')
   hold on

   Spaces = blanks(4);
   S = [    'Iteration #' Spaces 'Successes' Spaces ...
            'Failures' Spaces 'Samples'];

   Listbox = uicontrol('Style','Listbox','Units','Normalized',...
            'Position',[AxPos(1) .05 AxPos(3) .3],'Back',[1 1 1],...
            'Fore',[0 0 0],'String',{S},'Fontname','FixedWidth');

   HistogramAxis = subplot(3,1,2);
      set(HistogramAxis,'units','normalized','xlim',[min(a,b) max(a,b)])
%*********************************************************


%*********************************************************
function Options = optionscheck(opts)
%  This function checks options and
%  returns a structure of options along with the 
%  extra parameters the user wants to pass to the 
%  integrand function.

%  Set the defaults here:
Options.tol = 1e-6;
Options.display = 'off';
[Options.nodes Options.weights] = feval(@gausskronrod,7);

if isempty(opts); return; end

%  Loop through the fields of Options and
%  set the defaults as appropriate
Fields   = {'tol','display','nodes','weights'};
for i = 1:length(Fields)
    f = Fields{i};
    if isfield(opts,f) & ~isempty(getfield(opts,f))
        Options = setfield(Options,f,getfield(opts,f));
    end
end