| [f,a,b,result,abserr,resabs,resasc]=qk61(f,a,b,result,abserr,resabs,resasc); |
function [f,a,b,result,abserr,resabs,resasc]=qk61(f,a,b,result,abserr,resabs,resasc);
%***BEGIN PROLOGUE QK61
%***PURPOSE To compute I = Integral of F over (A,B) with error
% estimate
% J = Integral of ABS(F) over (A,B)
%***LIBRARY SLATEC (QUADPACK)
%***CATEGORY H2A1A2
%***TYPE SINGLE PRECISION (QK61-S, DQK61-D)
%***KEYWORDS 61-POINT GAUSS-KRONROD RULES, QUADPACK, QUADRATURE
%***AUTHOR Piessens, Robert
% Applied Mathematics and Programming Division
% K. U. Leuven
% de Doncker, Elise
% Applied Mathematics and Programming Division
% K. U. Leuven
%***DESCRIPTION
%
% Integration rule
% Standard fortran subroutine
% Real version
%
%
% PARAMETERS
% ON ENTRY
% F - Real
% 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 calling program.
%
% A - Real
% Lower limit of integration
%
% B - Real
% Upper limit of integration
%
% ON RETURN
% RESULT - Real
% Approximation to the integral I
% RESULT is computed by applying the 61-point
% Kronrod rule (RESK) obtained by optimal addition of
% abscissae to the 30-point Gauss rule (RESG).
%
% ABSERR - Real
% Estimate of the modulus of the absolute error,
% which should equal or exceed ABS(I-RESULT)
%
% RESABS - Real
% Approximation to the integral J
%
% RESASC - Real
% Approximation to the integral of ABS(F-I/(B-A))
%
%***REFERENCES (NONE)
%***ROUTINES CALLED R1MACH
%***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)
%***end PROLOGUE QK61
%
persistent absc centr dhlgth epmach fc firstCall fsum fv1 fv2 fval1 fval2 hlgth j jtw jtwm1 resg resk reskh uflow wg wgk xgk ; 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(fc), fc=0; end;
if isempty(fsum), fsum=0; end;
if isempty(fval1), fval1=0; end;
if isempty(fval2), fval2=0; end;
if isempty(fv1), fv1=zeros(1,30); end;
if isempty(fv2), fv2=zeros(1,30); end;
if isempty(hlgth), hlgth=0; end;
if isempty(resg), resg=0; end;
if isempty(resk), resk=0; end;
if isempty(reskh), reskh=0; end;
if isempty(uflow), uflow=0; end;
if isempty(wg), wg=zeros(1,15); end;
if isempty(wgk), wgk=zeros(1,31); end;
if isempty(xgk), xgk=zeros(1,31); end;
if isempty(j), j=0; end;
if isempty(jtw), jtw=0; end;
if isempty(jtwm1), jtwm1=0; end;
%
%
% THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE
% INTERVAL (-1,1). BECAUSE OF SYMMETRY ONLY THE POSITIVE
% ABSCISSAE AND THEIR CORRESPONDING WEIGHTS ARE GIVEN.
%
% XGK - ABSCISSAE OF THE 61-POINT KRONROD RULE
% XGK(2), XGK(4) ... ABSCISSAE OF THE 30-POINT
% GAUSS RULE
% XGK(1), XGK(3) ... OPTIMALLY ADDED ABSCISSAE
% TO THE 30-POINT GAUSS RULE
%
% WGK - WEIGHTS OF THE 61-POINT KRONROD RULE
%
% WG - WEIGHTS OF THE 30-POINT GAUSS RULE
%
if firstCall, xgk(1) =[0.9994844100504906e+00]; end;
if firstCall, xgk(2) =[0.9968934840746495e+00]; end;
if firstCall, xgk(3) =[0.9916309968704046e+00]; end;
if firstCall, xgk(4) =[0.9836681232797472e+00]; end;
if firstCall, xgk(5) =[0.9731163225011263e+00]; end;
if firstCall, xgk(6) =[0.9600218649683075e+00]; end;
if firstCall, xgk(7)=[0.9443744447485600e+00]; end;
if firstCall, xgk(8) =[0.9262000474292743e+00]; end;
if firstCall, xgk(9) =[0.9055733076999078e+00]; end;
if firstCall, xgk(10)=[0.8825605357920527e+00]; end;
if firstCall, xgk(11) =[0.8572052335460611e+00]; end;
if firstCall, xgk(12) =[0.8295657623827684e+00]; end;
if firstCall, xgk(13) =[0.7997278358218391e+00]; end;
if firstCall, xgk(14) =[0.7677774321048262e+00]; end;
if firstCall, xgk(15) =[0.7337900624532268e+00]; end;
if firstCall, xgk(16) =[0.6978504947933158e+00]; end;
if firstCall, xgk(17) =[0.6600610641266270e+00]; end;
if firstCall, xgk(18) =[0.6205261829892429e+00]; end;
if firstCall, xgk(19) =[0.5793452358263617e+00]; end;
if firstCall, xgk(20)=[0.5366241481420199e+00]; end;
if firstCall, xgk(21) =[0.4924804678617786e+00]; end;
if firstCall, xgk(22) =[0.4470337695380892e+00]; end;
if firstCall, xgk(23) =[0.4004012548303944e+00]; end;
if firstCall, xgk(24) =[0.3527047255308781e+00]; end;
if firstCall, xgk(25) =[0.3040732022736251e+00]; end;
if firstCall, xgk(26) =[0.2546369261678898e+00]; end;
if firstCall, xgk(27) =[0.2045251166823099e+00]; end;
if firstCall, xgk(28) =[0.1538699136085835e+00]; end;
if firstCall, xgk(29) =[0.1028069379667370e+00]; end;
if firstCall, xgk(30) =[0.5147184255531770e-01]; end;
if firstCall, xgk(31)=[0.0e+00]; end;
if firstCall, wgk(1) =[0.1389013698677008e-02]; end;
if firstCall, wgk(2) =[0.3890461127099884e-02]; end;
if firstCall, wgk(3) =[0.6630703915931292e-02]; end;
if firstCall, wgk(4) =[0.9273279659517763e-02]; end;
if firstCall, wgk(5) =[0.1182301525349634e-01]; end;
if firstCall, wgk(6) =[0.1436972950704580e-01]; end;
if firstCall, wgk(7)=[0.1692088918905327e-01]; end;
if firstCall, wgk(8) =[0.1941414119394238e-01]; end;
if firstCall, wgk(9) =[0.2182803582160919e-01]; end;
if firstCall, wgk(10)=[0.2419116207808060e-01]; end;
if firstCall, wgk(11) =[0.2650995488233310e-01]; end;
if firstCall, wgk(12) =[0.2875404876504129e-01]; end;
if firstCall, wgk(13) =[0.3090725756238776e-01]; end;
if firstCall, wgk(14) =[0.3298144705748373e-01]; end;
if firstCall, wgk(15) =[0.3497933802806002e-01]; end;
if firstCall, wgk(16) =[0.3688236465182123e-01]; end;
if firstCall, wgk(17) =[0.3867894562472759e-01]; end;
if firstCall, wgk(18) =[0.4037453895153596e-01]; end;
if firstCall, wgk(19) =[0.4196981021516425e-01]; end;
if firstCall, wgk(20)=[0.4345253970135607e-01]; end;
if firstCall, wgk(21) =[0.4481480013316266e-01]; end;
if firstCall, wgk(22) =[0.4605923827100699e-01]; end;
if firstCall, wgk(23) =[0.4718554656929915e-01]; end;
if firstCall, wgk(24) =[0.4818586175708713e-01]; end;
if firstCall, wgk(25) =[0.4905543455502978e-01]; end;
if firstCall, wgk(26) =[0.4979568342707421e-01]; end;
if firstCall, wgk(27) =[0.5040592140278235e-01]; end;
if firstCall, wgk(28) =[0.5088179589874961e-01]; end;
if firstCall, wgk(29) =[0.5122154784925877e-01]; end;
if firstCall, wgk(30) =[0.5142612853745903e-01]; end;
if firstCall, wgk(31)=[0.5149472942945157e-01]; end;
if firstCall, wg(1) =[0.7968192496166606e-02]; end;
if firstCall, wg(2) =[0.1846646831109096e-01]; end;
if firstCall, wg(3) =[0.2878470788332337e-01]; end;
if firstCall, wg(4) =[0.3879919256962705e-01]; end;
if firstCall, wg(5) =[0.4840267283059405e-01]; end;
if firstCall, wg(6) =[0.5749315621761907e-01]; end;
if firstCall, wg(7) =[0.6597422988218050e-01]; end;
if firstCall, wg(8)=[0.7375597473770521e-01]; end;
if firstCall, wg(9) =[0.8075589522942022e-01]; end;
if firstCall, wg(10) =[0.8689978720108298e-01]; end;
if firstCall, wg(11) =[0.9212252223778613e-01]; end;
if firstCall, wg(12) =[0.9636873717464426e-01]; end;
if firstCall, wg(13) =[0.9959342058679527e-01]; end;
if firstCall, wg(14) =[0.1017623897484055e+00]; end;
if firstCall, wg(15)=[0.1028526528935588e+00]; end;
firstCall=0;
%
% LIST OF MAJOR VARIABLES
% -----------------------
%
% CENTR - MID POINT OF THE INTERVAL
% HLGTH - HALF-LENGTH OF THE INTERVAL
% ABSC - ABSCISSA
% FVAL* - FUNCTION VALUE
% RESG - RESULT OF THE 30-POINT GAUSS RULE
% RESK - RESULT OF THE 61-POINT KRONROD RULE
% RESKH - APPROXIMATION TO THE MEAN VALUE OF F
% OVER (A,B), I.E. TO I/(B-A)
%
% MACHINE DEPENDENT CONSTANTS
% ---------------------------
%
% EPMACH IS THE LARGEST RELATIVE SPACING.
% UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
%
%***FIRST EXECUTABLE STATEMENT QK61
[epmach ]=r1mach(4);
[uflow ]=r1mach(1);
%
centr = 0.5e+00.*(b+a);
hlgth = 0.5e+00.*(b-a);
dhlgth = abs(hlgth);
%
% COMPUTE THE 61-POINT KRONROD APPROXIMATION TO THE
% INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
%
resg = 0.0e+00;
fc = f(centr);
resk = wgk(31).*fc;
resabs = abs(resk);
for j = 1 : 15;
jtw = fix(j.*2);
absc = hlgth.*xgk(jtw);
fval1 = f(centr-absc);
fval2 = f(centr+absc);
fv1(jtw) = fval1;
fv2(jtw) = fval2;
fsum = fval1 + fval2;
resg = resg + wg(j).*fsum;
resk = resk + wgk(jtw).*fsum;
resabs = resabs + wgk(jtw).*(abs(fval1)+abs(fval2));
end; j = fix(15+1);
for j = 1 : 15;
jtwm1 = fix(j.*2 - 1);
absc = hlgth.*xgk(jtwm1);
fval1 = f(centr-absc);
fval2 = f(centr+absc);
fv1(jtwm1) = fval1;
fv2(jtwm1) = fval2;
fsum = fval1 + fval2;
resk = resk + wgk(jtwm1).*fsum;
resabs = resabs + wgk(jtwm1).*(abs(fval1)+abs(fval2));
end; j = fix(15+1);
reskh = resk.*0.5e+00;
resasc = wgk(31).*abs(fc-reskh);
for j = 1 : 30;
resasc = resasc + wgk(j).*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh));
end; j = fix(30+1);
result = resk.*hlgth;
resabs = resabs.*dhlgth;
resasc = resasc.*dhlgth;
abserr = abs((resk-resg).*hlgth);
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;
end
%DECK QMOMO
|
|
