| [f,a,b,epsabs,epsrel,result,abserr,neval,ier]=dqng(f,a,b,epsabs,epsrel,result,abserr,neval,ier); |
function [f,a,b,epsabs,epsrel,result,abserr,neval,ier]=dqng(f,a,b,epsabs,epsrel,result,abserr,neval,ier);
%***BEGIN PROLOGUE DQNG
%***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 doubleprecision (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
% doubleprecision VERSION
%
% 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
%
% 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.
%
% ON RETURN
% RESULT - doubleprecision
% 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 - doubleprecision
% 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 D1MACH, 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 DQNG
%
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(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(w21b), w21b=zeros(1,6); 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
%
%
% GAUSS-KRONROD-PATTERSON QUADRATURE COEFFICIENTS FOR USE IN
% QUADPACK ROUTINE QNG. THESE COEFFICIENTS WERE CALCULATED WITH
% 101 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON, BELL LABS, NOV 1981.
%
if firstCall, x1(1)=[0.973906528517171720077964012084452d0]; end;
if firstCall, x1(2)=[0.865063366688984510732096688423493d0]; end;
if firstCall, x1(3)=[0.679409568299024406234327365114874d0]; end;
if firstCall, x1(4)=[0.433395394129247190799265943165784d0]; end;
if firstCall, x1(5)=[0.148874338981631210884826001129720d0]; end;
if firstCall, w10(1)=[0.066671344308688137593568809893332d0]; end;
if firstCall, w10(2)=[0.149451349150580593145776339657697d0]; end;
if firstCall, w10(3)=[0.219086362515982043995534934228163d0]; end;
if firstCall, w10(4)=[0.269266719309996355091226921569469d0]; end;
if firstCall, w10(5)=[0.295524224714752870173892994651338d0]; end;
%
if firstCall, x2(1)=[0.995657163025808080735527280689003d0]; end;
if firstCall, x2(2)=[0.930157491355708226001207180059508d0]; end;
if firstCall, x2(3)=[0.780817726586416897063717578345042d0]; end;
if firstCall, x2(4)=[0.562757134668604683339000099272694d0]; end;
if firstCall, x2(5)=[0.294392862701460198131126603103866d0]; end;
if firstCall, w21a(1)=[0.032558162307964727478818972459390d0]; end;
if firstCall, w21a(2)=[0.075039674810919952767043140916190d0]; end;
if firstCall, w21a(3)=[0.109387158802297641899210590325805d0]; end;
if firstCall, w21a(4)=[0.134709217311473325928054001771707d0]; end;
if firstCall, w21a(5)=[0.147739104901338491374841515972068d0]; end;
if firstCall, w21b(1)=[0.011694638867371874278064396062192d0]; end;
if firstCall, w21b(2)=[0.054755896574351996031381300244580d0]; end;
if firstCall, w21b(3)=[0.093125454583697605535065465083366d0]; end;
if firstCall, w21b(4)=[0.123491976262065851077958109831074d0]; end;
if firstCall, w21b(5)=[0.142775938577060080797094273138717d0]; end;
if firstCall, w21b(6)=[0.149445554002916905664936468389821d0]; end;
%
if firstCall, x3(1)=[0.999333360901932081394099323919911d0]; end;
if firstCall, x3(2)=[0.987433402908088869795961478381209d0]; end;
if firstCall, x3(3)=[0.954807934814266299257919200290473d0]; end;
if firstCall, x3(4)=[0.900148695748328293625099494069092d0]; end;
if firstCall, x3(5)=[0.825198314983114150847066732588520d0]; end;
if firstCall, x3(6)=[0.732148388989304982612354848755461d0]; end;
if firstCall, x3(7)=[0.622847970537725238641159120344323d0]; end;
if firstCall, x3(8)=[0.499479574071056499952214885499755d0]; end;
if firstCall, x3(9)=[0.364901661346580768043989548502644d0]; end;
if firstCall, x3(10)=[0.222254919776601296498260928066212d0]; end;
if firstCall, x3(11)=[0.074650617461383322043914435796506d0]; end;
if firstCall, w43a(1)=[0.016296734289666564924281974617663d0]; end;
if firstCall, w43a(2)=[0.037522876120869501461613795898115d0]; end;
if firstCall, w43a(3)=[0.054694902058255442147212685465005d0]; end;
if firstCall, w43a(4)=[0.067355414609478086075553166302174d0]; end;
if firstCall, w43a(5)=[0.073870199632393953432140695251367d0]; end;
if firstCall, w43a(6)=[0.005768556059769796184184327908655d0]; end;
if firstCall, w43a(7)=[0.027371890593248842081276069289151d0]; end;
if firstCall, w43a(8)=[0.046560826910428830743339154433824d0]; end;
if firstCall, w43a(9)=[0.061744995201442564496240336030883d0]; end;
if firstCall, w43a(10)=[0.071387267268693397768559114425516d0]; end;
if firstCall, w43b(1)=[0.001844477640212414100389106552965d0]; end;
if firstCall, w43b(2)=[0.010798689585891651740465406741293d0]; end;
if firstCall, w43b(3)=[0.021895363867795428102523123075149d0]; end;
if firstCall, w43b(4)=[0.032597463975345689443882222526137d0]; end;
if firstCall, w43b(5)=[0.042163137935191811847627924327955d0]; end;
if firstCall, w43b(6)=[0.050741939600184577780189020092084d0]; end;
if firstCall, w43b(7)=[0.058379395542619248375475369330206d0]; end;
if firstCall, w43b(8)=[0.064746404951445885544689259517511d0]; end;
if firstCall, w43b(9)=[0.069566197912356484528633315038405d0]; end;
if firstCall, w43b(10)=[0.072824441471833208150939535192842d0]; end;
if firstCall, w43b(11)=[0.074507751014175118273571813842889d0]; end;
if firstCall, w43b(12)=[0.074722147517403005594425168280423d0]; end;
%
if firstCall, x4(1)=[0.999902977262729234490529830591582d0]; end;
if firstCall, x4(2)=[0.997989895986678745427496322365960d0]; end;
if firstCall, x4(3)=[0.992175497860687222808523352251425d0]; end;
if firstCall, x4(4)=[0.981358163572712773571916941623894d0]; end;
if firstCall, x4(5)=[0.965057623858384619128284110607926d0]; end;
if firstCall, x4(6)=[0.943167613133670596816416634507426d0]; end;
if firstCall, x4(7)=[0.915806414685507209591826430720050d0]; end;
if firstCall, x4(8)=[0.883221657771316501372117548744163d0]; end;
if firstCall, x4(9)=[0.845710748462415666605902011504855d0]; end;
if firstCall, x4(10)=[0.803557658035230982788739474980964d0]; end;
if firstCall, x4(11)=[0.757005730685495558328942793432020d0]; end;
if firstCall, x4(12)=[0.706273209787321819824094274740840d0]; end;
if firstCall, x4(13)=[0.651589466501177922534422205016736d0]; end;
if firstCall, x4(14)=[0.593223374057961088875273770349144d0]; end;
if firstCall, x4(15)=[0.531493605970831932285268948562671d0]; end;
if firstCall, x4(16)=[0.466763623042022844871966781659270d0]; end;
if firstCall, x4(17)=[0.399424847859218804732101665817923d0]; end;
if firstCall, x4(18)=[0.329874877106188288265053371824597d0]; end;
if firstCall, x4(19)=[0.258503559202161551802280975429025d0]; end;
if firstCall, x4(20)=[0.185695396568346652015917141167606d0]; end;
if firstCall, x4(21)=[0.111842213179907468172398359241362d0]; end;
if firstCall, x4(22)=[0.037352123394619870814998165437704d0]; end;
if firstCall, w87a(1)=[0.008148377384149172900002878448190d0]; end;
if firstCall, w87a(2)=[0.018761438201562822243935059003794d0]; end;
if firstCall, w87a(3)=[0.027347451050052286161582829741283d0]; end;
if firstCall, w87a(4)=[0.033677707311637930046581056957588d0]; end;
if firstCall, w87a(5)=[0.036935099820427907614589586742499d0]; end;
if firstCall, w87a(6)=[0.002884872430211530501334156248695d0]; end;
if firstCall, w87a(7)=[0.013685946022712701888950035273128d0]; end;
if firstCall, w87a(8)=[0.023280413502888311123409291030404d0]; end;
if firstCall, w87a(9)=[0.030872497611713358675466394126442d0]; end;
if firstCall, w87a(10)=[0.035693633639418770719351355457044d0]; end;
if firstCall, w87a(11)=[0.000915283345202241360843392549948d0]; end;
if firstCall, w87a(12)=[0.005399280219300471367738743391053d0]; end;
if firstCall, w87a(13)=[0.010947679601118931134327826856808d0]; end;
if firstCall, w87a(14)=[0.016298731696787335262665703223280d0]; end;
if firstCall, w87a(15)=[0.021081568889203835112433060188190d0]; end;
if firstCall, w87a(16)=[0.025370969769253827243467999831710d0]; end;
if firstCall, w87a(17)=[0.029189697756475752501446154084920d0]; end;
if firstCall, w87a(18)=[0.032373202467202789685788194889595d0]; end;
if firstCall, w87a(19)=[0.034783098950365142750781997949596d0]; end;
if firstCall, w87a(20)=[0.036412220731351787562801163687577d0]; end;
if firstCall, w87a(21)=[0.037253875503047708539592001191226d0]; end;
if firstCall, w87b(1)=[0.000274145563762072350016527092881d0]; end;
if firstCall, w87b(2)=[0.001807124155057942948341311753254d0]; end;
if firstCall, w87b(3)=[0.004096869282759164864458070683480d0]; end;
if firstCall, w87b(4)=[0.006758290051847378699816577897424d0]; end;
if firstCall, w87b(5)=[0.009549957672201646536053581325377d0]; end;
if firstCall, w87b(6)=[0.012329447652244853694626639963780d0]; end;
if firstCall, w87b(7)=[0.015010447346388952376697286041943d0]; end;
if firstCall, w87b(8)=[0.017548967986243191099665352925900d0]; end;
if firstCall, w87b(9)=[0.019938037786440888202278192730714d0]; end;
if firstCall, w87b(10)=[0.022194935961012286796332102959499d0]; end;
if firstCall, w87b(11)=[0.024339147126000805470360647041454d0]; end;
if firstCall, w87b(12)=[0.026374505414839207241503786552615d0]; end;
if firstCall, w87b(13)=[0.028286910788771200659968002987960d0]; end;
if firstCall, w87b(14)=[0.030052581128092695322521110347341d0]; end;
if firstCall, w87b(15)=[0.031646751371439929404586051078883d0]; end;
if firstCall, w87b(16)=[0.033050413419978503290785944862689d0]; end;
if firstCall, w87b(17)=[0.034255099704226061787082821046821d0]; end;
if firstCall, w87b(18)=[0.035262412660156681033782717998428d0]; end;
if firstCall, w87b(19)=[0.036076989622888701185500318003895d0]; end;
if firstCall, w87b(20)=[0.036698604498456094498018047441094d0]; end;
if firstCall, w87b(21)=[0.037120549269832576114119958413599d0]; end;
if firstCall, w87b(22)=[0.037334228751935040321235449094698d0]; end;
if firstCall, w87b(23)=[0.037361073762679023410321241766599d0]; 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 DQNG
[epmach ]=d1mach(4);
[uflow ]=d1mach(1);
%
% TEST ON VALIDITY OF PARAMETERS
% ------------------------------
%
result = 0.0d+00;
abserr = 0.0d+00;
neval = 0;
ier = 6;
if( epsabs>0.0d+00 || epsrel>=max(0.5d+02.*epmach,0.5d-28) )
hlgth = 0.5d+00.*(b-a);
dhlgth = abs(hlgth);
centr = 0.5d+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.0d+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.5d+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.0d+00 && abserr~=0.0d+00 )
abserr = resasc.*min(0.1d+01,(0.2d+03.*abserr./resasc).^1.5d+00);
end;
if( resabs>uflow./(0.5d+02.*epmach) )
abserr = max((epmach.*0.5d+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','DQNG','ABNORMAL RETURN',ier,0);
end
%DECK DQPSRT
|
|
