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)

v=quadgsqrt(...
function v=quadgsqrt(...
              func,type,a,b,norder,nsegs,varargin)
%
% v=quadgsqrt(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.
          % Use Gauss quadrature
    [dumy,bp,wf]=gcquad(...
      '',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.
          % Use Gauss quadrature
    [dumy,bp,wf]=gcquad(...
      '',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

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

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

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function integrates a general function using 
% a composite Gauss formula of arbitrary order. The 
% integral value is returned along with base points
% and weight factors obtained by an eigenvalue based 
% method. The integration interval is divided into
% mparts subintervals of equal length and integration
% over each part is performed with a Gauss formula
% making nquad function evaluations. Results are 
% exact for polynomials of degree up to 2*nquad-1.   
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% func         - name of a function to be integrated
%                having an argument list of the form
%                func(x,p1,p2,...) where any auxiliary
%                parameters p1,p2,.. are passed through
%                variable varargin. Use [ ] for the
%                function name if only the base points
%                and weight factors are needed.
% xlow,xhigh   - integration limits
% nquad        - order of Gauss formula chosen
% mparts       - number of subintervals selected in 
%                the composite integration
% varargin     - variable length parameter used to
%                pass additional arguments needed in
%                the integrand func
% val          - numerical value of the integral
% bp,wf        - vectors containing base points and 
%                weight factors in the composite
%                integral formula
%
% A typical calculation such as:
% Fun=inline('(sin(w*t).^2).*exp(c*t)','t','w','c');
% A=0; B=12; nquad=21; mparts=10; w=10; c=8;  
% [value,pcterr]=integrate(Fun,A,B,nquad,mparts,w,c);
% gives value = 1.935685556078172e+040 which is 
% accurate within an error of 1.9e-13 percent.
%
% User m functions called:  the function name passed 
%                           in the argument list

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

 if isempty(nquad),  nquad=10; end
 if isempty(mparts), mparts=1; end

% Compute base points and weight factors 
% for the single interval [-1,1]. (Ref:
% 'Methods of Numerical Integration' by
% P. Davis and P. Rabinowitz, page 93)

u=(1:nquad-1)./sqrt((2*(1:nquad-1)).^2-1);
[vc,bp]=eig(diag(u,-1)+diag(u,1));
[bp,k]=sort(diag(bp)); wf=2*vc(1,k)'.^2;

% Modify the base points and weight factors
% to apply for a  composite interval
d=(xhigh-xlow)/mparts;  d1=d/2;
dbp=d1*bp(:); dwf=d1*wf(:);  dr=d*(1:mparts);
cbp=dbp(:,ones(1,mparts))+ ...
dr(ones(nquad,1),:)+(xlow-d1);
cwf=dwf(:,ones(1,mparts)); wf=cwf(:); bp=cbp(:);

% Compute the integral
if isempty(func)
  val=[];
else
  f=feval(func,bp,varargin{:}); val=wf'*f(:);
end

Contact us