| [f,a,b,result,abserr,resabs,resasc]=qk15(f,a,b,result,abserr,resabs,resasc); |
function [f,a,b,result,abserr,resabs,resasc]=qk15(f,a,b,result,abserr,resabs,resasc);
%***BEGIN PROLOGUE QK15
%***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 (QK15-S, DQK15-D)
%***KEYWORDS 15-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 rules
% 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 15-POINT
% KRONROD RULE (RESK) obtained by optimal addition
% of abscissae to the 7-POINT GAUSS RULE(RESG).
%
% ABSERR - Real
% Estimate of the modulus of the absolute error,
% which should not exceed ABS(I-RESULT)
%
% RESABS - Real
% Approximation to the integral J
%
% RESASC - Real
% Approximation to the integral of ABS(F-I/(B-A))
% over (A,B)
%
%***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 QK15
%
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,7); end;
if isempty(fv2), fv2=zeros(1,7); 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,4); end;
if isempty(wgk), wgk=zeros(1,8); end;
if isempty(xgk), xgk=zeros(1,8); 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 15-POINT KRONROD RULE
% XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT
% GAUSS RULE
% XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY
% ADDED TO THE 7-POINT GAUSS RULE
%
% WGK - WEIGHTS OF THE 15-POINT KRONROD RULE
%
% WG - WEIGHTS OF THE 7-POINT GAUSS RULE
%
if firstCall, xgk(1) =[0.9914553711208126e+00]; end;
if firstCall, xgk(2) =[0.9491079123427585e+00]; end;
if firstCall, xgk(3) =[0.8648644233597691e+00]; end;
if firstCall, xgk(4) =[0.7415311855993944e+00]; end;
if firstCall, xgk(5) =[0.5860872354676911e+00]; end;
if firstCall, xgk(6) =[0.4058451513773972e+00]; end;
if firstCall, xgk(7)=[0.2077849550078985e+00]; end;
if firstCall, xgk(8)=[0.0e+00]; end;
if firstCall, wgk(1) =[0.2293532201052922e-01]; end;
if firstCall, wgk(2) =[0.6309209262997855e-01]; end;
if firstCall, wgk(3) =[0.1047900103222502e+00]; end;
if firstCall, wgk(4) =[0.1406532597155259e+00]; end;
if firstCall, wgk(5) =[0.1690047266392679e+00]; end;
if firstCall, wgk(6) =[0.1903505780647854e+00]; end;
if firstCall, wgk(7)=[0.2044329400752989e+00]; end;
if firstCall, wgk(8)=[0.2094821410847278e+00]; end;
if firstCall, wg(1) =[0.1294849661688697e+00]; end;
if firstCall, wg(2) =[0.2797053914892767e+00]; end;
if firstCall, wg(3) =[0.3818300505051189e+00]; end;
if firstCall, wg(4)=[0.4179591836734694e+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 7-POINT GAUSS FORMULA
% RESK - RESULT OF THE 15-POINT KRONROD FORMULA
% 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 QK15
[epmach ]=r1mach(4);
[uflow ]=r1mach(1);
%
centr = 0.5e+00.*(a+b);
hlgth = 0.5e+00.*(b-a);
dhlgth = abs(hlgth);
%
% COMPUTE THE 15-POINT KRONROD APPROXIMATION TO
% THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
%
fc = f(centr);
resg = fc.*wg(4);
resk = fc.*wgk(8);
resabs = abs(resk);
for j = 1 : 3;
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(3+1);
for j = 1 : 4;
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(4+1);
reskh = resk.*0.5e+00;
resasc = wgk(8).*abs(fc-reskh);
for j = 1 : 7;
resasc = resasc + wgk(j).*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh));
end; j = fix(7+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 QK15I
|
|
