| [f,a,b,epsabs,epsrel,result,abserr,neval,ier]=qng(f,a,b,epsabs,epsrel,result,abserr,neval,ier); |
function [f,a,b,epsabs,epsrel,result,abserr,neval,ier]=qng(f,a,b,epsabs,epsrel,result,abserr,neval,ier);
persistent absc centr dhlgth epmach fcentr firstCall fv1 fv2 fv3 fv4 fval fval1 fval2 hlgth ipx k l res10 res21 res43 res87 resabs resasc reskh savfun uflow w10 w21a w21b w43a w43b w87a w87b x1 x2 x3 x4 ; if isempty(firstCall),firstCall=1;end;
if isempty(w21b), w21b=zeros(1,6); end;
%***BEGIN PROLOGUE QNG
%***PURPOSE The routine calculates an approximation result to a
% given definite integral I = integral of F over (A,B),
% hopefully satisfying following claim for accuracy
% ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
%***LIBRARY SLATEC (QUADPACK)
%***CATEGORY H2A1A1
%***TYPE SINGLE PRECISION (QNG-S, DQNG-D)
%***KEYWORDS AUTOMATIC INTEGRATOR, GAUSS-KRONROD(PATTERSON) RULES,
% NONADAPTIVE, QUADPACK, QUADRATURE, SMOOTH INTEGRAND
%***AUTHOR Piessens, Robert
% Applied Mathematics and Programming Division
% K. U. Leuven
% de Doncker, Elise
% Applied Mathematics and Programming Division
% K. U. Leuven
%***DESCRIPTION
%
% NON-ADAPTIVE INTEGRATION
% STANDARD FORTRAN SUBROUTINE
% REAL VERSION
%
% F - Real version
% 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 - Real version
% Lower limit of integration
%
% B - Real version
% Upper limit of integration
%
% EPSABS - Real
% Absolute accuracy requested
% EPSREL - Real
% 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.
%
% ON RETURN
% RESULT - Real
% Approximation to the integral I
% Result is obtained by applying the 21-POINT
% GAUSS-KRONROD RULE (RES21) obtained by optimal
% addition of abscissae to the 10-POINT GAUSS RULE
% (RES10), or by applying the 43-POINT RULE (RES43)
% obtained by optimal addition of abscissae to the
% 21-POINT GAUSS-KRONROD RULE, or by applying the
% 87-POINT RULE (RES87) obtained by optimal addition
% of abscissae to the 43-POINT RULE.
%
% ABSERR - Real
% Estimate of the modulus of the absolute error,
% which should EQUAL or EXCEED ABS(I-RESULT)
%
% NEVAL - Integer
% Number of integrand evaluations
%
% IER - 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. It is
% assumed that the requested accuracy has
% not been achieved.
% ERROR MESSAGES
% IER = 1 The maximum number of steps has been
% executed. The integral is probably too
% difficult to be calculated by DQNG.
% = 6 The input is invalid, because
% EPSABS.LE.0 AND
% EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28).
% RESULT, ABSERR and NEVAL are set to zero.
%
%***REFERENCES (NONE)
%***ROUTINES CALLED R1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 800101 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890531 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
%***end PROLOGUE QNG
%
if isempty(absc), absc=0; end;
if isempty(centr), centr=0; end;
if isempty(dhlgth), dhlgth=0; end;
if isempty(epmach), epmach=0; end;
if isempty(fcentr), fcentr=0; end;
if isempty(fval), fval=0; end;
if isempty(fval1), fval1=0; end;
if isempty(fval2), fval2=0; end;
if isempty(fv1), fv1=zeros(1,5); end;
if isempty(fv2), fv2=zeros(1,5); end;
if isempty(fv3), fv3=zeros(1,5); end;
if isempty(fv4), fv4=zeros(1,5); end;
if isempty(hlgth), hlgth=0; end;
if isempty(res10), res10=0; end;
if isempty(res21), res21=0; end;
if isempty(res43), res43=0; end;
if isempty(res87), res87=0; end;
if isempty(resabs), resabs=0; end;
if isempty(resasc), resasc=0; end;
if isempty(reskh), reskh=0; end;
if isempty(savfun), savfun=zeros(1,21); end;
if isempty(uflow), uflow=0; end;
if isempty(w10), w10=zeros(1,5); end;
if isempty(w21a), w21a=zeros(1,5); end;
if isempty(w43a), w43a=zeros(1,10); end;
if isempty(w43b), w43b=zeros(1,12); end;
if isempty(w87a), w87a=zeros(1,21); end;
if isempty(w87b), w87b=zeros(1,23); end;
if isempty(x1), x1=zeros(1,5); end;
if isempty(x2), x2=zeros(1,5); end;
if isempty(x3), x3=zeros(1,11); end;
if isempty(x4), x4=zeros(1,22); end;
if isempty(ipx), ipx=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
%
%
% THE FOLLOWING DATA STATEMENTS CONTAIN THE
% ABSCISSAE AND WEIGHTS OF THE INTEGRATION RULES USED.
%
% X1 ABSCISSAE COMMON TO THE 10-, 21-, 43-
% AND 87-POINT RULE
% X2 ABSCISSAE COMMON TO THE 21-, 43- AND
% 87-POINT RULE
% X3 ABSCISSAE COMMON TO THE 43- AND 87-POINT
% RULE
% X4 ABSCISSAE OF THE 87-POINT RULE
% W10 WEIGHTS OF THE 10-POINT FORMULA
% W21A WEIGHTS OF THE 21-POINT FORMULA FOR
% ABSCISSAE X1
% W21B WEIGHTS OF THE 21-POINT FORMULA FOR
% ABSCISSAE X2
% W43A WEIGHTS OF THE 43-POINT FORMULA FOR
% ABSCISSAE X1, X3
% W43B WEIGHTS OF THE 43-POINT FORMULA FOR
% ABSCISSAE X3
% W87A WEIGHTS OF THE 87-POINT FORMULA FOR
% ABSCISSAE X1, X2, X3
% W87B WEIGHTS OF THE 87-POINT FORMULA FOR
% ABSCISSAE X4
%
if firstCall, x1(1) =[0.9739065285171717e+00]; end;
if firstCall, x1(2) =[0.8650633666889845e+00]; end;
if firstCall, x1(3) =[0.6794095682990244e+00]; end;
if firstCall, x1(4) =[0.4333953941292472e+00]; end;
if firstCall, x1(5)=[0.1488743389816312e+00]; end;
if firstCall, x2(1) =[0.9956571630258081e+00]; end;
if firstCall, x2(2) =[0.9301574913557082e+00]; end;
if firstCall, x2(3) =[0.7808177265864169e+00]; end;
if firstCall, x2(4) =[0.5627571346686047e+00]; end;
if firstCall, x2(5)=[0.2943928627014602e+00]; end;
if firstCall, x3(1) =[0.9993333609019321e+00]; end;
if firstCall, x3(2) =[0.9874334029080889e+00]; end;
if firstCall, x3(3) =[0.9548079348142663e+00]; end;
if firstCall, x3(4) =[0.9001486957483283e+00]; end;
if firstCall, x3(5) =[0.8251983149831142e+00]; end;
if firstCall, x3(6) =[0.7321483889893050e+00]; end;
if firstCall, x3(7) =[0.6228479705377252e+00]; end;
if firstCall, x3(8) =[0.4994795740710565e+00]; end;
if firstCall, x3(9) =[0.3649016613465808e+00]; end;
if firstCall, x3(10) =[0.2222549197766013e+00]; end;
if firstCall, x3(11)=[0.7465061746138332e-01]; end;
if firstCall, x4(1) =[0.9999029772627292e+00]; end;
if firstCall, x4(2) =[0.9979898959866787e+00]; end;
if firstCall, x4(3) =[0.9921754978606872e+00]; end;
if firstCall, x4(4) =[0.9813581635727128e+00]; end;
if firstCall, x4(5) =[0.9650576238583846e+00]; end;
if firstCall, x4(6) =[0.9431676131336706e+00]; end;
if firstCall, x4(7) =[0.9158064146855072e+00]; end;
if firstCall, x4(8) =[0.8832216577713165e+00]; end;
if firstCall, x4(9) =[0.8457107484624157e+00]; end;
if firstCall, x4(10) =[0.8035576580352310e+00]; end;
if firstCall, x4(11) =[0.7570057306854956e+00]; end;
if firstCall, x4(12) =[0.7062732097873218e+00]; end;
if firstCall, x4(13) =[0.6515894665011779e+00]; end;
if firstCall, x4(14) =[0.5932233740579611e+00]; end;
if firstCall, x4(15) =[0.5314936059708319e+00]; end;
if firstCall, x4(16) =[0.4667636230420228e+00]; end;
if firstCall, x4(17) =[0.3994248478592188e+00]; end;
if firstCall, x4(18) =[0.3298748771061883e+00]; end;
if firstCall, x4(19) =[0.2585035592021616e+00]; end;
if firstCall, x4(20) =[0.1856953965683467e+00]; end;
if firstCall, x4(21)=[0.1118422131799075e+00]; end;
if firstCall, x4(22)=[0.3735212339461987e-01]; end;
if firstCall, w10(1) =[0.6667134430868814e-01]; end;
if firstCall, w10(2) =[0.1494513491505806e+00]; end;
if firstCall, w10(3) =[0.2190863625159820e+00]; end;
if firstCall, w10(4) =[0.2692667193099964e+00]; end;
if firstCall, w10(5)=[0.2955242247147529e+00]; end;
if firstCall, w21a(1) =[0.3255816230796473e-01]; end;
if firstCall, w21a(2) =[0.7503967481091995e-01]; end;
if firstCall, w21a(3) =[0.1093871588022976e+00]; end;
if firstCall, w21a(4) =[0.1347092173114733e+00]; end;
if firstCall, w21a(5)=[0.1477391049013385e+00]; end;
if firstCall, w21b(1) =[0.1169463886737187e-01]; end;
if firstCall, w21b(2) =[0.5475589657435200e-01]; end;
if firstCall, w21b(3) =[0.9312545458369761e-01]; end;
if firstCall, w21b(4) =[0.1234919762620659e+00]; end;
if firstCall, w21b(5) =[0.1427759385770601e+00]; end;
if firstCall, w21b(6)=[0.1494455540029169e+00]; end;
if firstCall, w43a(1) =[0.1629673428966656e-01]; end;
if firstCall, w43a(2) =[0.3752287612086950e-01]; end;
if firstCall, w43a(3) =[0.5469490205825544e-01]; end;
if firstCall, w43a(4) =[0.6735541460947809e-01]; end;
if firstCall, w43a(5) =[0.7387019963239395e-01]; end;
if firstCall, w43a(6) =[0.5768556059769796e-02]; end;
if firstCall, w43a(7) =[0.2737189059324884e-01]; end;
if firstCall, w43a(8) =[0.4656082691042883e-01]; end;
if firstCall, w43a(9) =[0.6174499520144256e-01]; end;
if firstCall, w43a(10)=[0.7138726726869340e-01]; end;
if firstCall, w43b(1) =[0.1844477640212414e-02]; end;
if firstCall, w43b(2) =[0.1079868958589165e-01]; end;
if firstCall, w43b(3) =[0.2189536386779543e-01]; end;
if firstCall, w43b(4) =[0.3259746397534569e-01]; end;
if firstCall, w43b(5) =[0.4216313793519181e-01]; end;
if firstCall, w43b(6) =[0.5074193960018458e-01]; end;
if firstCall, w43b(7) =[0.5837939554261925e-01]; end;
if firstCall, w43b(8) =[0.6474640495144589e-01]; end;
if firstCall, w43b(9) =[0.6956619791235648e-01]; end;
if firstCall, w43b(10) =[0.7282444147183321e-01]; end;
if firstCall, w43b(11) =[0.7450775101417512e-01]; end;
if firstCall, w43b(12)=[0.7472214751740301e-01]; end;
if firstCall, w87a(1) =[0.8148377384149173e-02]; end;
if firstCall, w87a(2) =[0.1876143820156282e-01]; end;
if firstCall, w87a(3) =[0.2734745105005229e-01]; end;
if firstCall, w87a(4) =[0.3367770731163793e-01]; end;
if firstCall, w87a(5) =[0.3693509982042791e-01]; end;
if firstCall, w87a(6) =[0.2884872430211531e-02]; end;
if firstCall, w87a(7) =[0.1368594602271270e-01]; end;
if firstCall, w87a(8) =[0.2328041350288831e-01]; end;
if firstCall, w87a(9) =[0.3087249761171336e-01]; end;
if firstCall, w87a(10) =[0.3569363363941877e-01]; end;
if firstCall, w87a(11) =[0.9152833452022414e-03]; end;
if firstCall, w87a(12)=[0.5399280219300471e-02]; end;
if firstCall, w87a(13) =[0.1094767960111893e-01]; end;
if firstCall, w87a(14) =[0.1629873169678734e-01]; end;
if firstCall, w87a(15) =[0.2108156888920384e-01]; end;
if firstCall, w87a(16) =[0.2537096976925383e-01]; end;
if firstCall, w87a(17) =[0.2918969775647575e-01]; end;
if firstCall, w87a(18) =[0.3237320246720279e-01]; end;
if firstCall, w87a(19) =[0.3478309895036514e-01]; end;
if firstCall, w87a(20) =[0.3641222073135179e-01]; end;
if firstCall, w87a(21)=[0.3725387550304771e-01]; end;
if firstCall, w87b(1) =[0.2741455637620724e-03]; end;
if firstCall, w87b(2) =[0.1807124155057943e-02]; end;
if firstCall, w87b(3) =[0.4096869282759165e-02]; end;
if firstCall, w87b(4) =[0.6758290051847379e-02]; end;
if firstCall, w87b(5) =[0.9549957672201647e-02]; end;
if firstCall, w87b(6) =[0.1232944765224485e-01]; end;
if firstCall, w87b(7) =[0.1501044734638895e-01]; end;
if firstCall, w87b(8) =[0.1754896798624319e-01]; end;
if firstCall, w87b(9) =[0.1993803778644089e-01]; end;
if firstCall, w87b(10) =[0.2219493596101229e-01]; end;
if firstCall, w87b(11) =[0.2433914712600081e-01]; end;
if firstCall, w87b(12)=[0.2637450541483921e-01]; end;
if firstCall, w87b(13) =[0.2828691078877120e-01]; end;
if firstCall, w87b(14) =[0.3005258112809270e-01]; end;
if firstCall, w87b(15) =[0.3164675137143993e-01]; end;
if firstCall, w87b(16) =[0.3305041341997850e-01]; end;
if firstCall, w87b(17) =[0.3425509970422606e-01]; end;
if firstCall, w87b(18) =[0.3526241266015668e-01]; end;
if firstCall, w87b(19) =[0.3607698962288870e-01]; end;
if firstCall, w87b(20) =[0.3669860449845609e-01]; end;
if firstCall, w87b(21) =[0.3712054926983258e-01]; end;
if firstCall, w87b(22) =[0.3733422875193504e-01]; end;
if firstCall, w87b(23)=[0.3736107376267902e-01]; end;
firstCall=0;
%
% LIST OF MAJOR VARIABLES
% -----------------------
%
% CENTR - MID POINT OF THE INTEGRATION INTERVAL
% HLGTH - HALF-LENGTH OF THE INTEGRATION INTERVAL
% FCENTR - FUNCTION VALUE AT MID POINT
% ABSC - ABSCISSA
% FVAL - FUNCTION VALUE
% SAVFUN - ARRAY OF FUNCTION VALUES WHICH
% HAVE ALREADY BEEN COMPUTED
% RES10 - 10-POINT GAUSS RESULT
% RES21 - 21-POINT KRONROD RESULT
% RES43 - 43-POINT RESULT
% RES87 - 87-POINT RESULT
% RESABS - APPROXIMATION TO THE INTEGRAL OF ABS(F)
% RESASC - APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A))
%
% MACHINE DEPENDENT CONSTANTS
% ---------------------------
%
% EPMACH IS THE LARGEST RELATIVE SPACING.
% UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
%
%***FIRST EXECUTABLE STATEMENT QNG
[epmach ]=r1mach(4);
[uflow ]=r1mach(1);
%
% TEST ON VALIDITY OF PARAMETERS
% ------------------------------
%
result = 0.0e+00;
abserr = 0.0e+00;
neval = 0;
ier = 6;
if( epsabs>0.0e+00 || epsrel>=max(0.5e-14,0.5e+02.*epmach) )
hlgth = 0.5e+00.*(b-a);
dhlgth = abs(hlgth);
centr = 0.5e+00.*(b+a);
fcentr = f(centr);
neval = 21;
ier = 1;
%
% COMPUTE THE INTEGRAL USING THE 10- AND 21-POINT FORMULA.
%
for l = 1 : 3;
if( l==2 )
%
% COMPUTE THE INTEGRAL USING THE 43-POINT FORMULA.
%
res43 = w43b(12).*fcentr;
neval = 43;
for k = 1 : 10;
res43 = res43 + savfun(k).*w43a(k);
end; k = fix(10+1);
for k = 1 : 11;
ipx = fix(ipx + 1);
absc = hlgth.*x3(k);
fval = f(absc+centr) + f(centr-absc);
res43 = res43 + fval.*w43b(k);
savfun(ipx) = fval;
end; k = fix(11+1);
%
% TEST FOR CONVERGENCE.
%
result = res43.*hlgth;
abserr = abs((res43-res21).*hlgth);
elseif( l==3 ) ;
%
% COMPUTE THE INTEGRAL USING THE 87-POINT FORMULA.
%
res87 = w87b(23).*fcentr;
neval = 87;
for k = 1 : 21;
res87 = res87 + savfun(k).*w87a(k);
end; k = fix(21+1);
for k = 1 : 22;
absc = hlgth.*x4(k);
res87 = res87 + w87b(k).*(f(absc+centr)+f(centr-absc));
end; k = fix(22+1);
result = res87.*hlgth;
abserr = abs((res87-res43).*hlgth);
else;
res10 = 0.0e+00;
res21 = w21b(6).*fcentr;
resabs = w21b(6).*abs(fcentr);
for k = 1 : 5;
absc = hlgth.*x1(k);
fval1 = f(centr+absc);
fval2 = f(centr-absc);
fval = fval1 + fval2;
res10 = res10 + w10(k).*fval;
res21 = res21 + w21a(k).*fval;
resabs = resabs + w21a(k).*(abs(fval1)+abs(fval2));
savfun(k) = fval;
fv1(k) = fval1;
fv2(k) = fval2;
end; k = fix(5+1);
ipx = 5;
for k = 1 : 5;
ipx = fix(ipx + 1);
absc = hlgth.*x2(k);
fval1 = f(centr+absc);
fval2 = f(centr-absc);
fval = fval1 + fval2;
res21 = res21 + w21b(k).*fval;
resabs = resabs + w21b(k).*(abs(fval1)+abs(fval2));
savfun(ipx) = fval;
fv3(k) = fval1;
fv4(k) = fval2;
end; k = fix(5+1);
%
% TEST FOR CONVERGENCE.
%
result = res21.*hlgth;
resabs = resabs.*dhlgth;
reskh = 0.5e+00.*res21;
resasc = w21b(6).*abs(fcentr-reskh);
for k = 1 : 5;
resasc = resasc + w21a(k).*(abs(fv1(k)-reskh)+abs(fv2(k)-reskh)) + w21b(k).*(abs(fv3(k)-reskh)+abs(fv4(k)-reskh));
end; k = fix(5+1);
abserr = abs((res21-res10).*hlgth);
resasc = resasc.*dhlgth;
end;
if( resasc~=0.0e+00 && abserr~=0.0e+00 )
abserr = resasc.*min(0.1e+01,(0.2e+03.*abserr./resasc).^1.5e+00);
end;
if( resabs>uflow./(0.5e+02.*epmach) )
abserr = max((epmach.*0.5e+02).*resabs,abserr);
end;
if( abserr<=max(epsabs,epsrel.*abs(result)) )
ier = 0;
end;
% ***JUMP OUT OF DO-LOOP
if( ier==0 )
return;
end;
end; l = fix(3+1);
end;
[dumvar1,dumvar2,dumvar3,ier]=xermsg('SLATEC','QNG','ABNORMAL RETURN',ier,0);
end
%DECK QPDOC
|
|
