No BSD License

Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

Howard Wilson (view profile)

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