| [f,a,b,c,epsabs,epsrel,limit,result,abserr,neval,ier,alist,blist,rlist,elist,iord,last]=dqawce(f,a,b,c,epsabs,epsrel,limit,result,abserr,neval,ier,alist,blist,rlist,elist,iord,last); |
function [f,a,b,c,epsabs,epsrel,limit,result,abserr,neval,ier,alist,blist,rlist,elist,iord,last]=dqawce(f,a,b,c,epsabs,epsrel,limit,result,abserr,neval,ier,alist,blist,rlist,elist,iord,last);
%***BEGIN PROLOGUE DQAWCE
%***PURPOSE The routine calculates an approximation result to a
% CAUCHY PRINCIPAL VALUE I = Integral of F*W over (A,B)
% (W(X) = 1/(X-C), (C.NE.A, C.NE.B), hopefully satisfying
% following claim for accuracy
% ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I))
%***LIBRARY SLATEC (QUADPACK)
%***CATEGORY H2A2A1, J4
%***TYPE doubleprecision (QAWCE-S, DQAWCE-D)
%***KEYWORDS AUTOMATIC INTEGRATOR, CAUCHY PRINCIPAL VALUE,
% CLENSHAW-CURTIS METHOD, QUADPACK, QUADRATURE,
% SPECIAL-PURPOSE
%***AUTHOR Piessens, Robert
% Applied Mathematics and Programming Division
% K. U. Leuven
% de Doncker, Elise
% Applied Mathematics and Programming Division
% K. U. Leuven
%***DESCRIPTION
%
% Computation of a CAUCHY PRINCIPAL VALUE
% Standard fortran subroutine
% doubleprecision version
%
% PARAMETERS
% ON ENTRY
% F - doubleprecision
% function subprogram defining the integrand
% function F(X). The actual name for F needs to be
% declared E X T E R N A L in the driver program.
%
% A - doubleprecision
% Lower limit of integration
%
% B - doubleprecision
% Upper limit of integration
%
% C - doubleprecision
% Parameter in the WEIGHT function, C.NE.A, C.NE.B
% If C = A OR C = B, the routine will end with
% IER = 6.
%
% EPSABS - doubleprecision
% Absolute accuracy requested
% EPSREL - doubleprecision
% Relative accuracy requested
% If EPSABS.LE.0
% and EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
% the routine will end with IER = 6.
%
% LIMIT - Integer
% Gives an upper bound on the number of subintervals
% in the partition of (A,B), LIMIT.GE.1
%
% ON RETURN
% RESULT - doubleprecision
% Approximation to the integral
%
% ABSERR - doubleprecision
% Estimate of the modulus of the absolute error,
% which should equal or exceed ABS(I-RESULT)
%
% NEVAL - Integer
% Number of integrand evaluations
%
% IER - Integer
% IER = 0 Normal and reliable termination of the
% routine. It is assumed that the requested
% accuracy has been achieved.
% IER.GT.0 Abnormal termination of the routine
% the estimates for integral and error are
% less reliable. It is assumed that the
% requested accuracy has not been achieved.
% ERROR MESSAGES
% IER = 1 Maximum number of subdivisions allowed
% has been achieved. One can allow more sub-
% divisions by increasing the value of
% LIMIT. However, if this yields no
% improvement it is advised to analyze the
% the integrand, in order to determine the
% the integration difficulties. If the
% position of a local difficulty can be
% determined (e.g. SINGULARITY,
% DISCONTINUITY within the interval) one
% will probably gain from splitting up the
% interval at this point and calling
% appropriate integrators on the subranges.
% = 2 The occurrence of roundoff error is detec-
% ted, which prevents the requested
% tolerance from being achieved.
% = 3 Extremely bad integrand behaviour
% occurs at some interior points of
% the integration interval.
% = 6 The input is invalid, because
% C = A or C = B or
% (EPSABS.LE.0 and
% EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
% or LIMIT.LT.1.
% RESULT, ABSERR, NEVAL, RLIST(1), ELIST(1),
% IORD(1) and LAST are set to zero. ALIST(1)
% and BLIST(1) are set to A and B
% respectively.
%
% ALIST - doubleprecision
% Vector of dimension at least LIMIT, the first
% LAST elements of which are the left
% end points of the subintervals in the partition
% of the given integration range (A,B)
%
% BLIST - doubleprecision
% Vector of dimension at least LIMIT, the first
% LAST elements of which are the right
% end points of the subintervals in the partition
% of the given integration range (A,B)
%
% RLIST - doubleprecision
% Vector of dimension at least LIMIT, the first
% LAST elements of which are the integral
% approximations on the subintervals
%
% ELIST - doubleprecision
% Vector of dimension LIMIT, the first LAST
% elements of which are the moduli of the absolute
% error estimates on the subintervals
%
% IORD - Integer
% Vector of dimension at least LIMIT, the first K
% elements of which are pointers to the error
% estimates over the subintervals, so that
% ELIST(IORD(1)), ..., ELIST(IORD(K)) with K = LAST
% If LAST.LE.(LIMIT/2+2), and K = LIMIT+1-LAST
% otherwise, form a decreasing sequence
%
% LAST - Integer
% Number of subintervals actually produced in
% the subdivision process
%
%***REFERENCES (NONE)
%***ROUTINES CALLED D1MACH, DQC25C, DQPSRT
%***REVISION HISTORY (YYMMDD)
% 800101 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 890831 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
%***end PROLOGUE DQAWCE
%
persistent a1 a2 aa area area1 area12 area2 b1 b2 bb epmach errbnd errmax erro12 error1 error2 errsum iroff1 iroff2 k krule maxerr nev nrmax uflow ;
if isempty(aa), aa=0; end;
if isempty(area), area=0; end;
if isempty(area1), area1=0; end;
if isempty(area12), area12=0; end;
if isempty(area2), area2=0; end;
if isempty(a1), a1=0; end;
if isempty(a2), a2=0; end;
if isempty(bb), bb=0; end;
if isempty(b1), b1=0; end;
if isempty(b2), b2=0; end;
if isempty(epmach), epmach=0; end;
if isempty(errbnd), errbnd=0; end;
if isempty(errmax), errmax=0; end;
if isempty(error1), error1=0; end;
if isempty(erro12), erro12=0; end;
if isempty(error2), error2=0; end;
if isempty(errsum), errsum=0; end;
if isempty(uflow), uflow=0; end;
if isempty(iroff1), iroff1=0; end;
if isempty(iroff2), iroff2=0; end;
if isempty(k), k=0; end;
if isempty(krule), krule=0; end;
if isempty(maxerr), maxerr=0; end;
if isempty(nev), nev=0; end;
if isempty(nrmax), nrmax=0; end;
%
alist_shape=size(alist);alist=reshape(alist,1,[]);
blist_shape=size(blist);blist=reshape(blist,1,[]);
rlist_shape=size(rlist);rlist=reshape(rlist,1,[]);
elist_shape=size(elist);elist=reshape(elist,1,[]);
iord_shape=size(iord);iord=reshape(iord,1,[]);
%
%
% LIST OF MAJOR VARIABLES
% -----------------------
%
% ALIST - LIST OF LEFT END POINTS OF ALL SUBINTERVALS
% CONSIDERED UP TO NOW
% BLIST - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS
% CONSIDERED UP TO NOW
% RLIST(I) - APPROXIMATION TO THE INTEGRAL OVER
% (ALIST(I),BLIST(I))
% ELIST(I) - ERROR ESTIMATE APPLYING TO RLIST(I)
% MAXERR - POINTER TO THE INTERVAL WITH LARGEST
% ERROR ESTIMATE
% ERRMAX - ELIST(MAXERR)
% AREA - SUM OF THE INTEGRALS OVER THE SUBINTERVALS
% ERRSUM - SUM OF THE ERRORS OVER THE SUBINTERVALS
% ERRBND - REQUESTED ACCURACY MAX(EPSABS,EPSREL*
% ABS(RESULT))
% *****1 - VARIABLE FOR THE LEFT SUBINTERVAL
% *****2 - VARIABLE FOR THE RIGHT SUBINTERVAL
% LAST - INDEX FOR SUBDIVISION
%
%
% MACHINE DEPENDENT CONSTANTS
% ---------------------------
%
% EPMACH IS THE LARGEST RELATIVE SPACING.
% UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
%
%***FIRST EXECUTABLE STATEMENT DQAWCE
[epmach ]=d1mach(4);
[uflow ]=d1mach(1);
%
%
% TEST ON VALIDITY OF PARAMETERS
% ------------------------------
%
ier = 6;
neval = 0;
last = 0;
alist(1) = a;
blist(1) = b;
rlist(1) = 0.0d+00;
elist(1) = 0.0d+00;
iord(1) = 0;
result = 0.0d+00;
abserr = 0.0d+00;
if( ~(c==a || c==b ||(epsabs<=0.0d+00 && epsrel<max(0.5d+02.*epmach,0.5d-28))) )
%
% FIRST APPROXIMATION TO THE INTEGRAL
% -----------------------------------
%
aa = a;
bb = b;
if( a>b )
aa = b;
bb = a;
end;
ier = 0;
krule = 1;
[f,aa,bb,c,result,abserr,krule,neval]=dqc25c(f,aa,bb,c,result,abserr,krule,neval);
last = 1;
rlist(1) = result;
elist(1) = abserr;
iord(1) = 1;
alist(1) = a;
blist(1) = b;
%
% TEST ON ACCURACY
%
errbnd = max(epsabs,epsrel.*abs(result));
if( limit==1 )
ier = 1;
end;
if( abserr>=min(0.1d-01.*abs(result),errbnd) && ier~=1 )
%
% INITIALIZATION
% --------------
%
alist(1) = aa;
blist(1) = bb;
rlist(1) = result;
errmax = abserr;
maxerr = 1;
area = result;
errsum = abserr;
nrmax = 1;
iroff1 = 0;
iroff2 = 0;
%
% MAIN DO-LOOP
% ------------
%
for last = 2 : limit;
%
% BISECT THE SUBINTERVAL WITH NRMAX-TH LARGEST
% ERROR ESTIMATE.
%
a1 = alist(maxerr);
b1 = 0.5d+00.*(alist(maxerr)+blist(maxerr));
b2 = blist(maxerr);
if( c<=b1 && c>a1 )
b1 = 0.5d+00.*(c+b2);
end;
if( c>b1 && c<b2 )
b1 = 0.5d+00.*(a1+c);
end;
a2 = b1;
krule = 2;
[f,a1,b1,c,area1,error1,krule,nev]=dqc25c(f,a1,b1,c,area1,error1,krule,nev);
neval = fix(neval + nev);
[f,a2,b2,c,area2,error2,krule,nev]=dqc25c(f,a2,b2,c,area2,error2,krule,nev);
neval = fix(neval + nev);
%
% IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
% AND ERROR AND TEST FOR ACCURACY.
%
area12 = area1 + area2;
erro12 = error1 + error2;
errsum = errsum + erro12 - errmax;
area = area + area12 - rlist(maxerr);
if( abs(rlist(maxerr)-area12)<0.1d-04.*abs(area12) &&erro12>=0.99d+00.*errmax && krule==0 )
iroff1 = fix(iroff1 +1);
end;
if( last>10 && erro12>errmax && krule==0 )
iroff2 = fix(iroff2 + 1);
end;
rlist(maxerr) = area1;
rlist(last) = area2;
errbnd = max(epsabs,epsrel.*abs(area));
if( errsum>errbnd )
%
% TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
%
if( iroff1>=6 && iroff2>20 )
ier = 2;
end;
%
% SET ERROR FLAG IN THE CASE THAT NUMBER OF INTERVAL
% BISECTIONS EXCEEDS LIMIT.
%
if( last==limit )
ier = 1;
end;
%
% SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
% AT A POINT OF THE INTEGRATION RANGE.
%
if( max(abs(a1),abs(b2))<=(0.1d+01+0.1d+03.*epmach).*(abs(a2)+0.1d+04.*uflow) )
ier = 3;
end;
end;
%
% APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
%
if( error2>error1 )
alist(maxerr) = a2;
alist(last) = a1;
blist(last) = b1;
rlist(maxerr) = area2;
rlist(last) = area1;
elist(maxerr) = error2;
elist(last) = error1;
else;
alist(last) = a2;
blist(maxerr) = b1;
blist(last) = b2;
elist(maxerr) = error1;
elist(last) = error2;
end;
%
% CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
% IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
% WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
%
[limit,last,maxerr,errmax,elist,iord,nrmax]=dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax);
% ***JUMP OUT OF DO-LOOP
if( ier~=0 || errsum<=errbnd )
break;
end;
end;
%
% COMPUTE FINAL RESULT.
% ---------------------
%
result = 0.0d+00;
for k = 1 : last;
result = result + rlist(k);
end; k = fix(last+1);
abserr = errsum;
end;
if( aa==b )
result = -result;
end;
end;
alist_shape=zeros(alist_shape);alist_shape(:)=alist(1:numel(alist_shape));alist=alist_shape;
blist_shape=zeros(blist_shape);blist_shape(:)=blist(1:numel(blist_shape));blist=blist_shape;
rlist_shape=zeros(rlist_shape);rlist_shape(:)=rlist(1:numel(rlist_shape));rlist=rlist_shape;
elist_shape=zeros(elist_shape);elist_shape(:)=elist(1:numel(elist_shape));elist=elist_shape;
iord_shape=zeros(iord_shape);iord_shape(:)=iord(1:numel(iord_shape));iord=iord_shape;
end
%DECK DQAWC
|
|
