GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dqk15i.f
Go to the documentation of this file.
1  SUBROUTINE dqk15i(F,BOUN,INF,A,B,RESULT,ABSERR,RESABS,RESASC,
2  1 IERR)
3 C***BEGIN PROLOGUE DQK15I
4 C***DATE WRITTEN 800101 (YYMMDD)
5 C***REVISION DATE 830518 (YYMMDD)
6 C***CATEGORY NO. H2A3A2,H2A4A2
7 C***KEYWORDS 15-POINT TRANSFORMED GAUSS-KRONROD RULES
8 C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
9 C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
10 C***PURPOSE THE ORIGINAL (INFINITE INTEGRATION RANGE IS MAPPED
11 C ONTO THE INTERVAL (0,1) AND (A,B) IS A PART OF (0,1).
12 C IT IS THE PURPOSE TO COMPUTE
13 C I = INTEGRAL OF TRANSFORMED INTEGRAND OVER (A,B),
14 C J = INTEGRAL OF ABS(TRANSFORMED INTEGRAND) OVER (A,B).
15 C***DESCRIPTION
16 C
17 C INTEGRATION RULE
18 C STANDARD FORTRAN SUBROUTINE
19 C DOUBLE PRECISION VERSION
20 C
21 C PARAMETERS
22 C ON ENTRY
23 C F - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND
24 C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
25 C DECLARED E X T E R N A L IN THE CALLING PROGRAM.
26 C
27 C BOUN - DOUBLE PRECISION
28 C FINITE BOUND OF ORIGINAL INTEGRATION
29 C RANGE (SET TO ZERO IF INF = +2)
30 C
31 C INF - INTEGER
32 C IF INF = -1, THE ORIGINAL INTERVAL IS
33 C (-INFINITY,BOUND),
34 C IF INF = +1, THE ORIGINAL INTERVAL IS
35 C (BOUND,+INFINITY),
36 C IF INF = +2, THE ORIGINAL INTERVAL IS
37 C (-INFINITY,+INFINITY) AND
38 C THE INTEGRAL IS COMPUTED AS THE SUM OF TWO
39 C INTEGRALS, ONE OVER (-INFINITY,0) AND ONE OVER
40 C (0,+INFINITY).
41 C
42 C A - DOUBLE PRECISION
43 C LOWER LIMIT FOR INTEGRATION OVER SUBRANGE
44 C OF (0,1)
45 C
46 C B - DOUBLE PRECISION
47 C UPPER LIMIT FOR INTEGRATION OVER SUBRANGE
48 C OF (0,1)
49 C
50 C ON RETURN
51 C RESULT - DOUBLE PRECISION
52 C APPROXIMATION TO THE INTEGRAL I
53 C RESULT IS COMPUTED BY APPLYING THE 15-POINT
54 C KRONROD RULE(RESK) OBTAINED BY OPTIMAL ADDITION
55 C OF ABSCISSAE TO THE 7-POINT GAUSS RULE(RESG).
56 C
57 C ABSERR - DOUBLE PRECISION
58 C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
59 C WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
60 C
61 C RESABS - DOUBLE PRECISION
62 C APPROXIMATION TO THE INTEGRAL J
63 C
64 C RESASC - DOUBLE PRECISION
65 C APPROXIMATION TO THE INTEGRAL OF
66 C ABS((TRANSFORMED INTEGRAND)-I/(B-A)) OVER (A,B)
67 C
68 C***REFERENCES (NONE)
69 C***ROUTINES CALLED D1MACH
70 C***END PROLOGUE DQK15I
71 C
72  DOUBLE PRECISION A,ABSC,ABSC1,ABSC2,ABSERR,B,BOUN,CENTR,DABS,DINF,
73  * dmax1,dmin1,d1mach,epmach,fc,fsum,fval1,fval2,fv1,fv2,hlgth,
74  * resabs,resasc,resg,resk,reskh,result,tabsc1,tabsc2,uflow,wg,wgk,
75  * xgk,fvalt
76  INTEGER INF,J
77  EXTERNAL f
78 C
79  dimension fv1(7),fv2(7),xgk(8),wgk(8),wg(8)
80 C
81 C THE ABSCISSAE AND WEIGHTS ARE SUPPLIED FOR THE INTERVAL
82 C (-1,1). BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND
83 C THEIR CORRESPONDING WEIGHTS ARE GIVEN.
84 C
85 C XGK - ABSCISSAE OF THE 15-POINT KRONROD RULE
86 C XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT
87 C GAUSS RULE
88 C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY
89 C ADDED TO THE 7-POINT GAUSS RULE
90 C
91 C WGK - WEIGHTS OF THE 15-POINT KRONROD RULE
92 C
93 C WG - WEIGHTS OF THE 7-POINT GAUSS RULE, CORRESPONDING
94 C TO THE ABSCISSAE XGK(2), XGK(4), ...
95 C WG(1), WG(3), ... ARE SET TO ZERO.
96 C
97  DATA wg(1) / 0.0d0 /
98  DATA wg(2) / 0.1294849661 6886969327 0611432679 082d0 /
99  DATA wg(3) / 0.0d0 /
100  DATA wg(4) / 0.2797053914 8927666790 1467771423 780d0 /
101  DATA wg(5) / 0.0d0 /
102  DATA wg(6) / 0.3818300505 0511894495 0369775488 975d0 /
103  DATA wg(7) / 0.0d0 /
104  DATA wg(8) / 0.4179591836 7346938775 5102040816 327d0 /
105 C
106  DATA xgk(1) / 0.9914553711 2081263920 6854697526 329d0 /
107  DATA xgk(2) / 0.9491079123 4275852452 6189684047 851d0 /
108  DATA xgk(3) / 0.8648644233 5976907278 9712788640 926d0 /
109  DATA xgk(4) / 0.7415311855 9939443986 3864773280 788d0 /
110  DATA xgk(5) / 0.5860872354 6769113029 4144838258 730d0 /
111  DATA xgk(6) / 0.4058451513 7739716690 6606412076 961d0 /
112  DATA xgk(7) / 0.2077849550 0789846760 0689403773 245d0 /
113  DATA xgk(8) / 0.0000000000 0000000000 0000000000 000d0 /
114 C
115  DATA wgk(1) / 0.0229353220 1052922496 3732008058 970d0 /
116  DATA wgk(2) / 0.0630920926 2997855329 0700663189 204d0 /
117  DATA wgk(3) / 0.1047900103 2225018383 9876322541 518d0 /
118  DATA wgk(4) / 0.1406532597 1552591874 5189590510 238d0 /
119  DATA wgk(5) / 0.1690047266 3926790282 6583426598 550d0 /
120  DATA wgk(6) / 0.1903505780 6478540991 3256402421 014d0 /
121  DATA wgk(7) / 0.2044329400 7529889241 4161999234 649d0 /
122  DATA wgk(8) / 0.2094821410 8472782801 2999174891 714d0 /
123 C
124 C
125 C LIST OF MAJOR VARIABLES
126 C -----------------------
127 C
128 C CENTR - MID POINT OF THE INTERVAL
129 C HLGTH - HALF-LENGTH OF THE INTERVAL
130 C ABSC* - ABSCISSA
131 C TABSC* - TRANSFORMED ABSCISSA
132 C FVAL* - FUNCTION VALUE
133 C RESG - RESULT OF THE 7-POINT GAUSS FORMULA
134 C RESK - RESULT OF THE 15-POINT KRONROD FORMULA
135 C RESKH - APPROXIMATION TO THE MEAN VALUE OF THE TRANSFORMED
136 C INTEGRAND OVER (A,B), I.E. TO I/(B-A)
137 C
138 C MACHINE DEPENDENT CONSTANTS
139 C ---------------------------
140 C
141 C EPMACH IS THE LARGEST RELATIVE SPACING.
142 C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
143 C
144 C***FIRST EXECUTABLE STATEMENT DQK15I
145  epmach = d1mach(4)
146  uflow = d1mach(1)
147  dinf = min0(1,inf)
148 C
149  centr = 0.5d+00*(a+b)
150  hlgth = 0.5d+00*(b-a)
151  tabsc1 = boun+dinf*(0.1d+01-centr)/centr
152  ierr = 0
153  CALL f(tabsc1,ierr,fval1)
154  IF (ierr .LT. 0) RETURN
155  IF(inf.EQ.2) THEN
156  CALL f(-tabsc1,ierr,fvalt)
157  IF (ierr .LT. 0) RETURN
158  fval1 = fval1+fvalt
159  ENDIF
160  fc = (fval1/centr)/centr
161 C
162 C COMPUTE THE 15-POINT KRONROD APPROXIMATION TO
163 C THE INTEGRAL, AND ESTIMATE THE ERROR.
164 C
165  resg = wg(8)*fc
166  resk = wgk(8)*fc
167  resabs = dabs(resk)
168  DO 10 j=1,7
169  absc = hlgth*xgk(j)
170  absc1 = centr-absc
171  absc2 = centr+absc
172  tabsc1 = boun+dinf*(0.1d+01-absc1)/absc1
173  tabsc2 = boun+dinf*(0.1d+01-absc2)/absc2
174  CALL f(tabsc1,ierr,fval1)
175  IF (ierr .LT. 0) RETURN
176  CALL f(tabsc2,ierr,fval2)
177  IF (ierr .LT. 0) RETURN
178  IF(inf.EQ.2) THEN
179  CALL f(-tabsc1,ierr,fvalt)
180  IF (ierr .LT. 0) RETURN
181  fval1 = fval1+fvalt
182  ENDIF
183  IF(inf.EQ.2) THEN
184  CALL f(-tabsc2,ierr,fvalt)
185  IF (ierr .LT. 0) RETURN
186  fval2 = fval2+fvalt
187  ENDIF
188  fval1 = (fval1/absc1)/absc1
189  fval2 = (fval2/absc2)/absc2
190  fv1(j) = fval1
191  fv2(j) = fval2
192  fsum = fval1+fval2
193  resg = resg+wg(j)*fsum
194  resk = resk+wgk(j)*fsum
195  resabs = resabs+wgk(j)*(dabs(fval1)+dabs(fval2))
196  10 CONTINUE
197  reskh = resk*0.5d+00
198  resasc = wgk(8)*dabs(fc-reskh)
199  DO 20 j=1,7
200  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
201  20 CONTINUE
202  result = resk*hlgth
203  resasc = resasc*hlgth
204  resabs = resabs*hlgth
205  abserr = dabs((resk-resg)*hlgth)
206  IF(resasc.NE.0.0d+00.AND.abserr.NE.0.d0) abserr = resasc*
207  * dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
208  IF(resabs.GT.uflow/(0.5d+02*epmach)) abserr = dmax1
209  * ((epmach*0.5d+02)*resabs,abserr)
210  RETURN
211  END
d1mach
double precision function d1mach(i)
Definition: d1mach.f:20
dimension
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)