# Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

### Howard Wilson (view profile)

14 Oct 2002 (Updated )

Companion Software (amamhlib)

```function [vg,tg,vL,tL,pctdiff]=sqrtquadtest
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function compares the accuracy and
% computation time for functions quadgsqrt
% where radical(x) is sqrt(x-a), sqrt(b-x), or
% sqrt((x-a)*(b-x))

%----------------------------------
% Program Output

% EVALUATING INTEGRALS WITH SQUARE ROOT TYPE
%     SINGULARITIES AT THE END POINTS

% Function integrated:
% ftest(x,u,v)=exp(u*x).*cos(v*x)

% a = 1   b = 4
% u = 3   v = 10

% 4.836504484e+003 -8.060993912e+003 -4.264510048e+003
% Computation time = 0.0159 sec.

% 4.836504484e+003 -8.060993912e+003 -4.264510048e+003
% Computation time = 7.03 sec.

% Percent difference for the two methods
% -3.6669e-012 -1.5344e-012 1.4929e-012
% >>

%-------------------------------------

% The test function
ftest=inline('exp(u*x).*cos(v*x)','x','u','v');

% Limits and function parameters
a=1; b=4; u=3; v=10;

nloop=100; tic;
for j=1:nloop
end
vg=[v1g,v2g,v3g]; tg=toc/nloop;
disp(' ')
disp('EVALUATING INTEGRALS WITH SQUARE ROOT TYPE')
disp('     SINGULARITIES AT THE END POINTS')
disp(' ')
disp('Function integrated:')
disp('ftest(x,u,v)=exp(u*x).*cos(v*x)')
disp(' ')
disp(['a = ',num2str(a),'   b = ',num2str(b)])
disp(['u = ',num2str(u),'   v = ',num2str(v)])
disp(' ')
fprintf('%17.9e %17.9e %17.9e\n',vg)
disp(['Computation time = ',num2str(tg),' sec.'])

tol=1e-12; tic;
vL=[v1L,v2L,v3L]; tL=toc;

disp(' ')
fprintf('%17.9e %17.9e %17.9e\n',vL)
disp(['Computation time = ',num2str(tL),' sec.'])

pctdiff=100*(vg-vL)./vL; disp(' ')
disp('Percent difference for the two methods')
fprintf('%13.4e %12.4e %12.4e\n',pctdiff)

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

func,type,a,b,norder,nsegs,varargin)
%
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function evaluates an integral having a
% square root type singularity at one or both ends
% of the integration interval a<x<b. Composite
% Gauss integration is used with func(x) treated
% as a polynomial of degree norder.
% The integrand has the form:
% func(x)/sqrt(x-a) if type==1.
% func(x)/sqrt(b-x) if type==2.
% func(x)/sqrt((x-a)*(b-x)) if type==3.
% The integration interval is subdivided into
% nsegs subintervals of equal length.
%
% func   - a character string or function handle
%          naming a function continuous in the
%          interval from x=a to x=b
% type   - 1 if the integrand is singular at x=a
%          2 if the integrand is singular at x=b
%          3 if the integrand is singular at both
%            x=a and x=b.
% a,b    - integration limits with b>a
% norder - polynomial interpolation order within
%          each interval. Lowest norder is 20.
% nsegs  - number of integration subintervals
%
% User m functions called: gcquad
%
% Reference: Abromowitz and Stegun, 'Handbook of
%            Mathematical Functions', Chapter 25
% -------------------------------------------

if nargin<6, nsegs=1; end;
if nargin<5, norder=50; end
switch type
case 1  % Singularity at the left end.
'',0,sqrt(b-a),norder+1,nsegs);
t=a+bp.^2; y=feval(func,t,varargin{:});
v=wf(:)'*y(:)*2;
case 2  % Singularity at the right end.
'',0,sqrt(b-a),norder+1,nsegs);
t=b-bp.^2; y=feval(func,t,varargin{:});
v=wf(:)'*y(:)*2;
case 3  % Singularity at both ends.
% Use Chebyshev integration
n=norder; bp=cos(pi/(2*n+2)*(1:2:2*n+1));
c1=(b+a)/2; c2=(b-a)/2; t=c1+c2*bp;
y=feval(func,t,varargin{:});
v=pi/(n+1)*sum(y);
end

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

%
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function uses the MATLAB integrator quadl
% to evaluate integrals having square root type
% singularities at one or both ends of the
% integration interval a < x < b.
% The integrand has the form:
% func(x)/sqrt(x-a) if type==1.
% func(x)/sqrt(b-x) if type==2.
% func(x)/sqrt((x-a)*(b-x)) if type==3.
%
% func   - the handle for a function continuous
%          from x=a to x=b
% type   - 1 if the integrand is singular at x=a
%          2 if the integrand is singular at x=b
%          3 if the integrand is singular at both
%            x=a and x=b.
% a,b    - integration limits with b > a

if nargin<6 | isempty(trace), trace=0; end
if nargin<5 | isempty(tol), tol=1e-8; end
if nargin<7
varargin{1}=type; varargin{2}=[a,b];
varargin{3}=fname;
else
n=length(varargin); c=[a,b]; varargin{n+1}=type;
varargin{n+2}=c; varargin{n+3}=fname;
end

if type==1 | type==2
tol,trace,varargin{:});
else
end

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

function u=fshift(x,varargin)
% u=fshift(x,varargin)
% This function shifts arguments to produce
% a nonsingular integrand called by quadl
N=length(varargin); fname=varargin{N};
c=varargin{N-1}; type=varargin{N-2};
a=c(1); b=c(2); c1=(b+a)/2; c2=(b-a)/2;

switch type
case 1, t=a+x.^2; case 2, t=b-x.^2;
case 3, t=c1+c2*cos(x);
end

if N>3, u=feval(fname,t,varargin{1:N-3});
else, u=feval(fname,t); end

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