No BSD License  

Highlights from
Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

image thumbnail

Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

by

Howard Wilson

 

14 Oct 2002 (Updated )

Companion Software (amamhlib)

[L,G,names]=quadtest(secs)
function [L,G,names]=quadtest(secs)
%
% [L,G,names]=quadtest(secs)
% ~~~~~~~~~~~~~~~~~~~~~~
% This program compares the accuracy and
% computation times for several integrals
% evaluated using quadl and gcquad
%
% secs  - the number of seconds each integration
%         is repeated to get accurate timing. The
%         default value is 60 seconds.
% L,G   - matrices with columns containing
%         results from quadl and from gcquad.
%         The matrices are structured as:
%         [IntegralValue,PercentError,...
%         FunctionEvaluations,ComputationSeconds]
% names - character matrix with rows
%         describing the functions 
%         which were integrated
%
% User functions called: ftest, gcquad
%---------------------------------------------

global nvals

if nargin==0, secs=60; end

fprintf('\nPRESS RETURN TO BEGIN COMPUTATION > ')
pause

% Summary of the five integrands used
names=strvcat('sqrt(x)','log(x)','humps(x)',...
              'exp(10*x).*cos(10*pi*x)',...
              'cos(20*pi*x-20*sin(pi*x))');
fprintf(['\n\nINTEGRATION TEST COMPARING',...
      ' FUNCTIONS QUADL AND GCQUAD\n'])
fprintf('\nThe functions being integrated are:\n')
disp(names)          
          
% Compute exact values of integrals          
exact=[2/3; -1; quadl(@humps,0,1,1e-12);
      real((exp(10+10*pi*i)-1)/(10+10*pi*i));
      besselj(20,20)];    
  
% Find time to make a loop and call the clock
nmax=5000; nclock=0; t0=clock;
while nclock<nmax
    nclock=nclock+1; tclock=etime(clock,t0);
end
tclock=tclock/nclock; 
          
% Evaluate each integral individually. Repeat
% the integrations for secs seconds to get
% accurate timing.  Save results in array L.
L=zeros(5,4); tol=1e-6; e=exact; warning off;
for k=1:5 
  nquad=0; tim=0; t0=clock;
  while tim<secs
    [v,nfuns]=quadl(@ftest,0,1,tol,[],k);
    nquad=nquad+1; tim=etime(clock,t0);
  end
  tim=tim/nquad-tclock; pe=100*(v/e(k)-1);
  L(k,:)=[v,nfuns,pe,tim];
end
warning on;

% Obtain time to compute base points and weight
% factors for a Gauss formula of order 100
nloop=100; t0=clock;
for j=1:nloop
  [dumy,bp,wf]=gcquad([],0,1,100,1);
end
tbpwf=etime(clock,t0)/nloop;

% Perform the Gauss integration using a
% vector integrand. Save results in array G
          
ngquad=0; tim=0; t0=clock; 
while tim<secs
  v=ftest(bp,6)*wf;
  ngquad=ngquad+1; tim=etime(clock,t0);
end
tim=tim/ngquad+tbpwf-tclock; pe=100*(v./e-1);
G=[v,100*ones(5,1),pe,tim/5*ones(5,1)];

format short e
disp(' ')
disp('            Results Using Function quadl')
disp(...
'    Integral     Function      Percent   Computation')
disp(...
'     values     evaluations     error     seconds')
disp(L)
disp('            Results Using Function gcquad')
disp(...
'   Integral     Function      Percent    Computation')
disp(...
'    values     evaluations     error       seconds')
disp(G)
format short
disp(['(Total time using quadl)/',...
      '(Total time using gcquad)'])
disp(['equals ',...
       num2str(sum(L(:,end))/sum(G(:,end)))])
disp(' ')
 
%=============================================          
          
function y=ftest(x,n)
% Integrands used by function quadl
global nvals
switch n
case 1, y=sqrt(x); case 2, y=log(x); 
case 3, y=humps(x); 
case 4, y=exp(10*x).*cos(10*pi*x);    
case 5, y=cos(20*pi*x-20*sin(pi*x));
otherwise
  x=x(:)'; y=[sqrt(x);log(x);humps(x);
              exp(10*x).*cos(10*pi*x);
              cos(20*pi*x-20*sin(pi*x))];
end
if n<6, nvals=nvals+length(x);  
else, nvals=nvals+5*length(x); end

%============================================= 

% function [val,bp,wf]=gcquad(func,xlow,...
%                    xhigh,nquad,mparts,varargin)
% See Appendix B

Contact us